Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Merged
merged 3 commits into from
Jul 11, 2014
Merged

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
20430b3
Enable filtering (filt) across columns of arrays
mbauman Jul 10, 2014
128c4a6
Relax filt signatures to allow heterogeneous types
mbauman Jul 10, 2014
06803c9
API change to fix the type instability
mbauman Jul 11, 2014
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
65 changes: 36 additions & 29 deletions base/dsp.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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!,
Expand All @@ -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")
Expand All @@ -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
Expand All @@ -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]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a case where N == 1 && col != 1?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, I see si could be a vector but x and out could be matrices.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yup, exactly. It allows the default case of zeros to just allocate a vector instead of an array, and that branch should be statically inlined.

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
Expand Down
10 changes: 5 additions & 5 deletions doc/stdlib/base.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
12 changes: 10 additions & 2 deletions test/dsp.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.]
Expand Down