diff --git a/base/dsp.jl b/base/dsp.jl index e9a761f0bc41d..095f0f4fa9370 100644 --- a/base/dsp.jl +++ b/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 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 diff --git a/doc/stdlib/base.rst b/doc/stdlib/base.rst index bc64fda3dc5ee..77369545f7833 100644 --- a/doc/stdlib/base.rst +++ b/doc/stdlib/base.rst @@ -4664,15 +4664,15 @@ calling functions from `FFTW `_. 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) diff --git a/test/dsp.jl b/test/dsp.jl index f48367940ce16..b87894bc30e20 100644 --- a/test/dsp.jl +++ b/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.]