RFC: filt improvements

NEWS.md

@@ -317,6 +317,7 @@ Library improvements
 
   * New macro `@evalpoly` for efficient inline evaluation of polynomials ([#7146]).
 
+  * The signal filtering function `filt` now accepts an optional initial filter state vector. A new in-place function `filt!` is also exported. ([#7513])
 
 Build improvements
 ------------------
@@ -870,3 +871,4 @@ Too numerous to mention.
 [#7146]: https://github.com/JuliaLang/julia/issues/7146
 [#7373]: https://github.com/JuliaLang/julia/issues/7373
 [#7435]: https://github.com/JuliaLang/julia/issues/7435
+[#7513]: https://github.com/JuliaLang/julia/issues/7513

base/dsp.jl

@@ -3,71 +3,73 @@ module DSP
 importall Base.FFTW
 import Base.FFTW.normalization
 
-export FFTW, filt, deconv, conv, conv2, xcorr, fftshift, ifftshift,
+export FFTW, filt, filt!, deconv, conv, conv2, xcorr, fftshift, ifftshift,
        dct, idct, dct!, idct!, plan_dct, plan_idct, plan_dct!, plan_idct!,
        # the rest are defined imported from FFTW:
        fft, bfft, ifft, rfft, brfft, irfft,
        plan_fft, plan_bfft, plan_ifft, plan_rfft, plan_brfft, plan_irfft,
        fft!, bfft!, ifft!, plan_fft!, plan_bfft!, plan_ifft!
 
-function filt{T<:Number}(b::Union(AbstractVector{T}, T),a::Union(AbstractVector{T}, T),x::AbstractVector{T})
-    if isempty(b); error("b must be non-empty"); end
-    if isempty(a); error("a must be non-empty"); end
-    if a[1]==0; error("a[1] must be nonzero"); end
+function filt{T<:Number}(b::Union(AbstractVector{T}, T), a::Union(AbstractVector{T}, T),
+                         x::AbstractVector{T}; si::AbstractVector{T}=zeros(T, max(length(a), length(b))-1))
+    filt!(Array(T, size(x)), b, a, x, si)
+end
+
+# in-place filtering: returns results in the out argument, which may shadow x
+# (and does so by default)
+function filt!{T<:Number}(b::Union(AbstractVector{T}, T), a::Union(AbstractVector{T}, T), x::AbstractVector{T};
+                          si::AbstractVector{T}=zeros(T, max(length(a), length(b))-1), out::AbstractVector{T}=x)
+    filt!(out, b, a, x, si)
+end
+function filt!{T<:Number}(out::AbstractVector{T}, b::Union(AbstractVector{T}, T), a::Union(AbstractVector{T}, T),
+                          x::AbstractVector{T}, si::AbstractVector{T})
+    isempty(b) && error("b must be non-empty")
+    isempty(a) && error("a must be non-empty")
+    a[1] == 0  && error("a[1] must be nonzero")
+    size(x) != size(out) && error("out size must match x")
 
     as = length(a)
     bs = length(b)
     sz = max(as, bs)
-
-    if sz == 1
-        return (b[1]/a[1]).*x
-    end
-
-    if bs<sz
-        newb = zeros(T,sz)
-        newb[1:bs] = b
-        b = newb
-    end
-
+    silen = sz - 1
     xs = size(x,1)
-    y = Array(T, xs)
-    silen = sz-1
-    si = zeros(T, silen)
 
+    xs == 0 && return out
+    sz == 1 && return scale!(out, x, b[1]/a[1]) # Simple scaling without memory
+
+    # Filter coefficient normalization
     if a[1] != 1
         norml = a[1]
         a ./= norml
         b ./= norml
     end
+    # Pad the coefficients with zeros if needed
+    bs<sz   && (b = copy!(zeros(T,sz), b))
+    1<as<sz && (a = copy!(zeros(T,sz), a))
 
+    si = copy!(Array(T, silen), si)
     if as > 1
-        if as<sz
-            newa = zeros(T,sz)
-            newa[1:as] = a
-            a = newa
-        end
-
-        @inbounds begin
-            for i=1:xs
-                y[i] = si[1] + b[1]*x[i]
-                for j=1:(silen-1)
-                    si[j] = si[j+1] + b[j+1]*x[i] - a[j+1]*y[i]
-                end
-                si[silen] = b[silen+1]*x[i] - a[silen+1]*y[i]
+        @inbounds for i=1:xs
+            xi = x[i]
+            val = si[1] + b[1]*xi
+            for j=1:(silen-1)
+                si[j] = si[j+1] + b[j+1]*xi - a[j+1]*val
             end
+            si[silen] = b[silen+1]*xi - a[silen+1]*val
+            out[i] = val
         end
     else
-        @inbounds begin
-            for i=1:xs
-                y[i] = si[1] + b[1]*x[i]
-                for j=1:(silen-1)
-                    si[j] = si[j+1] + b[j+1]*x[i]
-                end
-                si[silen] = b[silen+1]*x[i]
+        @inbounds for i=1:xs
+            xi = x[i]
+            val = si[1] + b[1]*xi
+            for j=1:(silen-1)
+                si[j] = si[j+1] + b[j+1]*xi
             end
+            si[silen] = b[silen+1]*xi
+            out[i] = val
         end
     end
-    return y
+    return out
 end
 
 function deconv{T}(b::StridedVector{T}, a::StridedVector{T})

base/exports.jl

@@ -936,6 +936,7 @@ export
     fft,
     fftshift,
     filt,
+    filt!,
     idct!,
     idct,
     ifft!,

base/linalg/generic.jl

@@ -11,14 +11,18 @@ end
 
 scale{R<:Real}(s::Complex, X::AbstractArray{R}) = scale(X, s)
 
-function generic_scale!(X::AbstractArray, s::Number)
+generic_scale!(X::AbstractArray, s::Number) = generic_scale!(X, X, s)
+function generic_scale!(C::AbstractArray, X::AbstractArray, s::Number)
+    length(C) == length(X) || error("C must be the same length as X")
     for i = 1:length(X)
-        @inbounds X[i] *= s
+        @inbounds C[i] = X[i]*s
     end
-    X
+    C
 end
-scale!(X::AbstractArray, s::Number) = generic_scale!(X, s)
-scale!(s::Number, X::AbstractArray) = generic_scale!(X, s)
+scale!(C::AbstractArray, s::Number, X::AbstractArray) = generic_scale!(C, X, s)
+scale!(C::AbstractArray, X::AbstractArray, s::Number) = generic_scale!(C, X, s)
+scale!(X::AbstractArray, s::Number) = generic_scale!(X, X, s)
+scale!(s::Number, X::AbstractArray) = generic_scale!(X, X, s)
 
 cross(a::AbstractVector, b::AbstractVector) = [a[2]*b[3]-a[3]*b[2], a[3]*b[1]-a[1]*b[3], a[1]*b[2]-a[2]*b[1]]
 

doc/stdlib/base.rst

@@ -4663,9 +4663,15 @@ calling functions from `FFTW <http://www.fftw.org>`_.
 
    Undoes the effect of ``fftshift``.
 
-.. function:: filt(b,a,x)
+.. function:: filt(b, a, x; si=zeros(max(length(a),length(b))-1))
 
-   Apply filter described by vectors ``a`` and ``b`` to vector ``x``.
+   Apply filter described by vectors ``a`` and ``b`` to vector ``x``, with an
+   optional initial filter state keyword argument ``si`` (defaults to zeros).
+
+.. function:: filt!(b, a, x; si=zeros(max(length(a),length(b))-1), out=x)
+
+   Same as :func:`filt`, but stores the result in the ``out`` keyword argument,
+   which may alias the input ``x`` to modify it in-place (it does so by default).
 
 .. function:: deconv(b,a)
 

test/dsp.jl

@@ -3,6 +3,15 @@ b = [1., 2., 3., 4.]
 x = [1., 1., 0., 1., 1., 0., 0., 0.]
 @test filt(b, 1., x)  == [1., 3., 5., 8., 7., 5., 7., 4.]
 @test filt(b, [1., -0.5], x)  == [1., 3.5, 6.75, 11.375, 12.6875, 11.34375, 12.671875, 10.3359375]
+# With ranges
+@test filt(b, 1., 1.0:10.0) == [1., 4., 10., 20., 30., 40., 50., 60., 70., 80.]
+@test filt(1.:4., 1., 1.0:10.0) == [1., 4., 10., 20., 30., 40., 50., 60., 70., 80.]
+# With initial conditions: a lowpass 5-pole butterworth filter with W_n = 0.25,
+# and a stable initial filter condition matched to the initial value
+b = [0.003279216306360201,0.016396081531801006,0.03279216306360201,0.03279216306360201,0.016396081531801006,0.003279216306360201]
+a = [1.0,-2.4744161749781606,2.8110063119115782,-1.703772240915465,0.5444326948885326,-0.07231566910295834]
+init = [0.9967207836936347,-1.4940914728163142,1.2841226760316475,-0.4524417279474106,0.07559488540931815]
+@test_approx_eq filt(b, a, ones(10), si=init) ones(10) # Shouldn't affect DC offset
 
 # Convolution
 a = [1., 2., 1., 2.]