Merge pull request #7560 from mbauman/filtarray

Filtering (filt) across columns of arrays and method signature tweaks

base/dsp.jl

@@ -2,6 +2,7 @@ module DSP
 
 importall Base.FFTW
 import Base.FFTW.normalization
+import Base.trailingsize
 
 export FFTW, filt, filt!, deconv, conv, conv2, xcorr, fftshift, ifftshift,
        dct, idct, dct!, idct!, plan_dct, plan_idct, plan_dct!, plan_idct!,
@@ -10,19 +11,17 @@ export FFTW, filt, filt!, deconv, conv, conv2, xcorr, fftshift, ifftshift,
        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}; si::AbstractVector{T}=zeros(T, max(length(a), length(b))-1))
-    filt!(Array(T, size(x)), b, a, x, si)
+_zerosi(b,a,T) = zeros(promote_type(eltype(b), eltype(a), T), max(length(a), length(b))-1)
+
+function filt{T,S}(b::Union(AbstractVector, Number), a::Union(AbstractVector, Number),
+                   x::AbstractArray{T}, si::AbstractArray{S}=_zerosi(b,a,T))
+    filt!(Array(promote_type(eltype(b), eltype(a), T, S), 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})
+function filt!{T,S,N}(out::AbstractArray, b::Union(AbstractVector, Number), a::Union(AbstractVector, Number),
+                      x::AbstractArray{T}, si::AbstractArray{S,N}=_zerosi(b,a,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")
@@ -33,6 +32,10 @@ function filt!{T<:Number}(out::AbstractVector{T}, b::Union(AbstractVector{T}, T)
     sz = max(as, bs)
     silen = sz - 1
     xs = size(x,1)
+    ncols = trailingsize(x,2)
+
+    size(si, 1) != silen && error("si must have max(length(a),length(b))-1 rows")
+    N > 1 && trailingsize(si,2) != ncols && error("si must either be a vector or have the same number of columns as x")
 
     xs == 0 && return out
     sz == 1 && return scale!(out, x, b[1]/a[1]) # Simple scaling without memory
@@ -44,29 +47,33 @@ function filt!{T<:Number}(out::AbstractVector{T}, b::Union(AbstractVector{T}, T)
         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))
+    bs<sz   && (b = copy!(zeros(eltype(b), sz), b))
+    1<as<sz && (a = copy!(zeros(eltype(a), sz), a))
 
-    si = copy!(Array(T, silen), si)
-    if as > 1
-        @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
+    initial_si = si
+    for col = 1:ncols
+        # Reset the filter state
+        si = initial_si[:, N > 1 ? col : 1]
+        if as > 1
+            @inbounds for i=1:xs
+                xi = x[i,col]
+                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,col] = val
             end
-            si[silen] = b[silen+1]*xi - a[silen+1]*val
-            out[i] = val
-        end
-    else
-        @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
+        else
+            @inbounds for i=1:xs
+                xi = x[i,col]
+                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,col] = val
             end
-            si[silen] = b[silen+1]*xi
-            out[i] = val
         end
     end
     return out

doc/stdlib/base.rst

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

test/dsp.jl

@@ -6,12 +6,20 @@ x = [1., 1., 0., 1., 1., 0., 0., 0.]
 # 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.]
+# Across an array is the same as channel-by-channel
+@test filt(b, 1., [x 1.0:8.0]) == [filt(b, 1., x) filt(b, 1., 1.0:8.0)]
+@test filt(b, [1., -0.5], [x 1.0:8.0]) == [filt(b, [1., -0.5], x) filt(b, [1., -0.5], 1.0:8.0)]
+si = zeros(3)
+@test filt(b, 1., [x 1.0:8.0], si) == [filt(b, 1., x, si) filt(b, 1., 1.0:8.0, si)]
+@test si == zeros(3) # Will likely fail if/when arrayviews are implemented
+si = [zeros(3) ones(3)]
+@test filt(b, 1., [x 1.0:8.0], si) == [filt(b, 1., x, zeros(3)) filt(b, 1., 1.0:8.0, ones(3))]
 # 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
+si = [0.9967207836936347,-1.4940914728163142,1.2841226760316475,-0.4524417279474106,0.07559488540931815]
+@test_approx_eq filt(b, a, ones(10), si) ones(10) # Shouldn't affect DC offset
 
 # Convolution
 a = [1., 2., 1., 2.]