Skip to content

Commit

Permalink
Faster filt!
Browse files Browse the repository at this point in the history
The array pointer loads cannot be hoisted out of the loop if this code
is in the same function that does the checking and potential
modification of b and a. This change gives a 1.2-1.5x perf gain. It
also mitigates potential effects of type instability of b and a due to
coefficient normalization.
  • Loading branch information
simonster committed Jul 23, 2014
1 parent d4f8c31 commit 3cf4fca
Showing 1 changed file with 29 additions and 20 deletions.
49 changes: 29 additions & 20 deletions base/dsp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -31,13 +31,12 @@ function filt!{T,S,N}(out::AbstractArray, b::Union(AbstractVector, Number), a::U
bs = length(b)
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
size(x,1) == 0 && return out
sz == 1 && return scale!(out, x, b[1]/a[1]) # Simple scaling without memory

# Filter coefficient normalization
Expand All @@ -55,30 +54,40 @@ function filt!{T,S,N}(out::AbstractArray, b::Union(AbstractVector, Number), a::U
# 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
_filt_iir!(out, b, a, x, si, col)
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
_filt_fir!(out, b, x, si, col)
end
end
return out
end

function _filt_iir!(out, b, a, x, si, col)
silen = length(si)
@inbounds for i=1:size(x, 1)
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
end

function _filt_fir!(out, b, x, si, col)
silen = length(si)
@inbounds for i=1:size(x, 1)
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
end

function deconv{T}(b::StridedVector{T}, a::StridedVector{T})
lb = size(b,1)
la = size(a,1)
Expand Down

0 comments on commit 3cf4fca

Please sign in to comment.