Skip to content

Commit bcca68b

Browse files
simonstersimonster
committed
Merge pull request #100 from JuliaDSP/sjk/stateful-filters
Add stateful filters
2 parents a009d8c + 8ddce88 commit bcca68b

File tree

6 files changed

+184
-39
lines changed

6 files changed

+184
-39
lines changed

‎NEWS.md‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,8 @@
1+
# DSP v0.0.8 Release Notes
2+
3+
- The `DF2TFilter` object provides a filter that preserves state
4+
between invocations ([#100](https://github.com/JuliaDSP/DSP.jl/pull/100)).
5+
16
# DSP v0.0.7 Release Notes
27

38
- Filter coefficient types have been renamed to distinguish them from implementations ([#96](https://github.com/JuliaDSP/DSP.jl/pull/96)):

‎doc/filters.rst‎

Lines changed: 41 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,19 @@
11
:mod:`Filters` - filter design and filtering
22
============================================
33

4-
Linear time-invariant filter representations
5-
--------------------------------------------
4+
DSP.jl differentiates between filter coefficients and filters. Filter
5+
coefficient objects specify the response of the filter in one of
6+
several standard forms. Filter objects carry the state of the filter
7+
together with filter coefficients in an implementable form
8+
(``PolynomialRatio``, ``Biquad``, or ``SecondOrderSections``).
9+
When invoked on a filter coefficient object, ``filt`` preserves state
10+
between invocations.
611

7-
DSP.jl supports common filter representations. Filters can be converted
8-
from one type to another using ``convert``.
12+
Linear time-invariant filter coefficient objects
13+
------------------------------------------------
14+
15+
DSP.jl supports common filter representations. Filter coefficients can
16+
be converted from one type to another using ``convert``.
917

1018
.. function:: ZeroPoleGain(z, p, k)
1119

@@ -45,28 +53,45 @@ from one type to another using ``convert``.
4553
sections and gain. ``biquads`` must be specified as a vector of
4654
``Biquads``.
4755

56+
57+
Filter objects
58+
----------------------------------
59+
60+
.. function:: DF2TFilter(coef[, si])
61+
62+
Construct a stateful direct form II transposed filter with
63+
coefficients ``coef``. ``si`` is an optional array representing the
64+
initial filter state (defaults to zeros). If ``f`` is a
65+
``PolynomialRatio``, ``Biquad``, or ``SecondOrderSections``,
66+
filtering is implemented directly. If ``f`` is a ``ZeroPoleGain``
67+
object, it is first converted to a ``SecondOrderSections`` object.
68+
69+
4870
Filter application
4971
------------------
5072

5173
.. function:: filt(f, x[, si])
5274

53-
Apply filter ``f`` along the first dimension of array ``x``. ``si``
54-
is an optional array representing the initial filter state
55-
(defaults to zeros). If ``f`` is a ``PolynomialRatio``, ``Biquad``,
56-
or ``SecondOrderSections``, filtering is implemented directly. If ``f`` is a
57-
``ZeroPoleGain``, it is first converted to an ``SecondOrderSections``.
75+
Apply filter or filter coefficients ``f`` along the first dimension
76+
of array ``x``. If ``f`` is a filter coefficient object, ``si``
77+
is an optional array representing the initial filter state (defaults
78+
to zeros). If ``f`` is a ``PolynomialRatio``, ``Biquad``, or
79+
``SecondOrderSections``, filtering is implemented directly. If
80+
``f`` is a ``ZeroPoleGain`` object, it is first converted to a
81+
``SecondOrderSections`` object.
5882

5983
.. function:: filt!(out, f, x[, si])
6084

6185
Same as :func:`filt()` but writes the result into the ``out``
6286
argument, which may alias the input ``x`` to modify it in-place.
6387

64-
.. function:: filtfilt(f, x)
88+
.. function:: filtfilt(coef, x)
6589

66-
Filter ``x`` in the forward and reverse directions. The initial
67-
state of the filter is computed so that its response to a step
68-
function is steady state. Before filtering, the data is
69-
extrapolated at both ends with an odd-symmetric extension of length
90+
Filter ``x`` in the forward and reverse directions using filter
91+
coefficients ``f``. The initial state of the filter is computed so
92+
that its response to a step function is steady state. Before
93+
filtering, the data is extrapolated at both ends with an
94+
odd-symmetric extension of length
7095
``3*(max(length(b), length(a))-1)``.
7196

7297
Because ``filtfilt`` applies the given filter twice, the effective
@@ -84,6 +109,7 @@ Filter application
84109
choosing the optimal algorithm based on the lengths of ``b`` and
85110
``x``.
86111

112+
87113
Filter design
88114
-------------
89115

@@ -95,6 +121,7 @@ Filter design
95121

96122
Construct a digital filter.
97123

124+
98125
Filter response types
99126
---------------------
100127

‎src/Filters/Filters.jl‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,11 @@
11
module Filters
22
using Polynomials, Compat, ..Util
33

4-
include("types.jl")
5-
export Filter, ZeroPoleGain, PolynomialRatio, Biquad, SecondOrderSections, coefa, coefb
4+
include("coefficients.jl")
5+
export FilterCoefficients, ZeroPoleGain, PolynomialRatio, Biquad, SecondOrderSections, coefa, coefb
66

77
include("filt.jl")
8-
export filtfilt, fftfilt, firfilt
8+
export DF2TFilter, filtfilt, fftfilt, firfilt
99

1010
include("design.jl")
1111
export FilterType, Butterworth, Chebyshev1, Chebyshev2, Elliptic,

src/Filters/types.jl renamed to src/Filters/coefficients.jl

File renamed without changes.

‎src/Filters/filt.jl‎

Lines changed: 120 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,26 @@ Base.filt(f::PolynomialRatio, x, si=_zerosi(f, x)) = filt(coefb(f), coefa(f), x,
1616
_zerosi{T,G,S}(f::SecondOrderSections{T,G}, x::AbstractArray{S}) =
1717
zeros(promote_type(T, G, S), 2, length(f.biquads))
1818

19+
# filt! algorithm (no checking, returns si)
20+
function _filt!{S,N}(out::AbstractArray, si::AbstractArray{S,N}, f::SecondOrderSections,
21+
x::AbstractArray, col::Int)
22+
g = f.g
23+
biquads = f.biquads
24+
n = length(biquads)
25+
@inbounds for i = 1:size(x, 1)
26+
yi = x[i, col]
27+
for fi = 1:n
28+
biquad = biquads[fi]
29+
xi = yi
30+
yi = si[1, fi] + biquad.b0*xi
31+
si[1, fi] = si[2, fi] + biquad.b1*xi - biquad.a1*yi
32+
si[2, fi] = biquad.b2*xi - biquad.a2*yi
33+
end
34+
out[i, col] = yi*g
35+
end
36+
si
37+
end
38+
1939
function Base.filt!{S,N}(out::AbstractArray, f::SecondOrderSections, x::AbstractArray,
2040
si::AbstractArray{S,N}=_zerosi(f, x))
2141
biquads = f.biquads
@@ -26,19 +46,10 @@ function Base.filt!{S,N}(out::AbstractArray, f::SecondOrderSections, x::Abstract
2646
error("si must be 2 x nbiquads or 2 x nbiquads x nsignals")
2747

2848
initial_si = si
49+
g = f.g
50+
n = length(biquads)
2951
for col = 1:ncols
30-
si = initial_si[:, :, N > 2 ? col : 1]
31-
@inbounds for i = 1:size(x, 1)
32-
yi = x[i, col]
33-
for fi = 1:length(biquads)
34-
biquad = biquads[fi]
35-
xi = yi
36-
yi = si[1, fi] + biquad.b0*xi
37-
si[1, fi] = si[2, fi] + biquad.b1*xi - biquad.a1*yi
38-
si[2, fi] = biquad.b2*xi - biquad.a2*yi
39-
end
40-
out[i, col] = yi*f.g
41-
end
52+
_filt!(out, initial_si[:, :, N > 2 ? col : 1], f, x, col)
4253
end
4354
out
4455
end
@@ -50,6 +61,20 @@ Base.filt{T,G,S<:Number}(f::SecondOrderSections{T,G}, x::AbstractArray{S}, si=_z
5061
_zerosi{T,S}(f::Biquad{T}, x::AbstractArray{S}) =
5162
zeros(promote_type(T, S), 2)
5263

64+
# filt! algorithm (no checking, returns si)
65+
function _filt!(out::AbstractArray, si1::Number, si2::Number, f::Biquad,
66+
x::AbstractArray, col::Int)
67+
@inbounds for i = 1:size(x, 1)
68+
xi = x[i, col]
69+
yi = si1 + f.b0*xi
70+
si1 = si2 + f.b1*xi - f.a1*yi
71+
si2 = f.b2*xi - f.a2*yi
72+
out[i, col] = yi
73+
end
74+
(si1, si2)
75+
end
76+
77+
# filt! variant that preserves si
5378
function Base.filt!{S,N}(out::AbstractArray, f::Biquad, x::AbstractArray,
5479
si::AbstractArray{S,N}=_zerosi(f, x))
5580
ncols = Base.trailingsize(x, 2)
@@ -60,15 +85,7 @@ function Base.filt!{S,N}(out::AbstractArray, f::Biquad, x::AbstractArray,
6085

6186
initial_si = si
6287
for col = 1:ncols
63-
si1 = initial_si[1, N > 1 ? col : 1]
64-
si2 = initial_si[2, N > 1 ? col : 1]
65-
@inbounds for i = 1:size(x, 1)
66-
xi = x[i, col]
67-
yi = si1 + f.b0*xi
68-
si1 = si2 + f.b1*xi - f.a1*yi
69-
si2 = f.b2*xi - f.a2*yi
70-
out[i, col] = yi
71-
end
88+
_filt!(out, initial_si[1, N > 1 ? col : 1], initial_si[2, N > 1 ? col : 1], f, x, col)
7289
end
7390
out
7491
end
@@ -80,6 +97,88 @@ Base.filt{T,S<:Number}(f::Biquad{T}, x::AbstractArray{S}, si=_zerosi(f, x)) =
8097
Base.filt(f::FilterCoefficients, x) = filt(convert(SecondOrderSections, f), x)
8198
Base.filt!(out, f::FilterCoefficients, x) = filt!(out, convert(SecondOrderSections, f), x)
8299

100+
#
101+
# Direct form II transposed filter with state
102+
#
103+
104+
immutable DF2TFilter{T<:FilterCoefficients,S<:Array}
105+
coef::T
106+
state::S
107+
108+
function DF2TFilter(coef::PolynomialRatio, state::Vector)
109+
length(state) == length(coef.a)-1 == length(coef.b)-1 ||
110+
throw(ArgumentError("length of state vector must match filter order"))
111+
new(coef, state)
112+
end
113+
function DF2TFilter(coef::SecondOrderSections, state::Matrix)
114+
(size(state, 1) == 2 && size(state, 2) == length(coef.biquads)) ||
115+
throw(ArgumentError("state must be 2 x nbiquads"))
116+
new(coef, state)
117+
end
118+
function DF2TFilter(coef::Biquad, state::Vector)
119+
length(state) == 2 || throw(ArgumentError("length of state must be 2"))
120+
new(coef, state)
121+
end
122+
end
123+
124+
## PolynomialRatio
125+
DF2TFilter{T,S}(coef::PolynomialRatio{T},
126+
state::Vector{S}=zeros(T, max(length(coef.a), length(coef.b))-1)) =
127+
DF2TFilter{PolynomialRatio{T}, Vector{S}}(coef, state)
128+
129+
function Base.filt!{T,S}(out::AbstractVector, f::DF2TFilter{PolynomialRatio{T},Vector{S}}, x::AbstractVector)
130+
length(x) != length(out) && throw(ArgumentError("out size must match x"))
131+
132+
si = f.state
133+
# Note: these are in the Polynomials.jl convention, which is
134+
# reversed WRT the usual DSP convention
135+
b = f.coef.b.a
136+
a = f.coef.a.a
137+
n = length(b)
138+
if n == 1
139+
scale!(out, x, b[1])
140+
else
141+
@inbounds for i=1:length(x)
142+
xi = x[i]
143+
val = si[1] + b[n]*xi
144+
for j=1:n-2
145+
si[j] = si[j+1] + b[n-j]*xi - a[n-j]*val
146+
end
147+
si[n-1] = b[1]*xi - a[1]*val
148+
out[i] = val
149+
end
150+
end
151+
out
152+
end
153+
154+
## SecondOrderSections
155+
DF2TFilter{T,G,S}(coef::SecondOrderSections{T,G},
156+
state::Matrix{S}=zeros(promote_type(T, G), 2, length(coef.biquads))) =
157+
DF2TFilter{SecondOrderSections{T,G}, Matrix{S}}(coef, state)
158+
159+
function Base.filt!{T,G,S}(out::AbstractVector, f::DF2TFilter{SecondOrderSections{T,G},Matrix{S}}, x::AbstractVector)
160+
length(x) != length(out) && throw(ArgumentError("out size must match x"))
161+
_filt!(out, f.state, f.coef, x, 1)
162+
out
163+
end
164+
165+
## Biquad
166+
DF2TFilter{T,S}(coef::Biquad{T}, state::Vector{S}=zeros(T, 2)) =
167+
DF2TFilter{Biquad{T}, Vector{S}}(coef, state)
168+
function Base.filt!{T,S}(out::AbstractVector, f::DF2TFilter{Biquad{T},Vector{S}}, x::AbstractVector)
169+
length(x) != length(out) && throw(ArgumentError("out size must match x"))
170+
si = f.state
171+
(si[1], si[2]) =_filt!(out, si[1], si[2], f.coef, x, 1)
172+
out
173+
end
174+
175+
# Variant that allocates the output
176+
Base.filt{T,S<:Array}(f::DF2TFilter{T,S}, x::AbstractVector) =
177+
filt!(Array(eltype(S), length(x)), f, x)
178+
179+
# Fall back to SecondOrderSections
180+
DF2TFilter(coef::FilterCoefficients) = DF2TFilter(convert(SecondOrderSections, coef))
181+
83182
#
84183
# filtfilt
85184
#

‎test/filt.jl‎

Lines changed: 15 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -21,22 +21,36 @@ for n = 1:6
2121
res = filt(sos, x)
2222

2323
# Test with filt with tf/sos
24-
@test_approx_eq res filt(tf, x)
24+
tfres = filt(tf, x)
25+
@test_approx_eq res tfres
2526
@test_approx_eq res filt!(similar(x), sos, x)
2627
@test_approx_eq res filt!(similar(x), tf, x)
2728

2829
# For <= 2 poles, test with biquads
2930
if n <= 2
3031
@test_approx_eq res filt(bq, x)
3132
@test_approx_eq res filt!(similar(x), bq, x)
33+
f = DF2TFilter(bq)
34+
@test tfres == [filt(f, x[1:50]); filt(f, x[51:end])]
3235
end
3336

3437
# Test that filt with zpk converts
3538
@test res == filt(zpk, x)
3639
@test res == filt!(similar(x), zpk, x)
40+
41+
# Test with DF2TFilter
42+
f = DF2TFilter(sos)
43+
@test res == [filt(f, x[1:50]); filt(f, x[51:end])]
44+
f = DF2TFilter(tf)
45+
@test tfres == [filt(f, x[1:50]); filt(f, x[51:end])]
46+
f = DF2TFilter(zpk)
47+
@test res == [filt(f, x[1:50]); filt(f, x[51:end])]
3748
end
3849
end
3950

51+
# Test simple scaling with DF2TFilter
52+
@test filt(DF2TFilter(PolynomialRatio([3.7], [4.2])), x) == scale(x, 3.7/4.2)
53+
4054
#
4155
# filtfilt
4256
#

0 commit comments

Comments
 (0)