filters.jl
This unit implements banks of band-pass filters interfacing the DSP.jl package. Those banks are used to create time-frequency analytic signal representations via the Hilbert transform.
At the bottom of this page you can find notes on DSP package useful functions.
FilterDesign
FilterDesign = Union{ZeroPoleGain, FIRWindow}
Those are possible filter designs implemented in the DSP package. For design methods see the DSP documentation.
FourierAnalysis.filterbank
— Functionfunction filterbank(x :: Vector{T},
sr :: Int,
bandwidth :: IntOrReal = 2;
filtkind :: FilterDesign = Butterworth(2),
fmin :: IntOrReal = bandwidth,
fmax :: IntOrReal = sr÷2,
⏩ :: Bool = true) where T<:Real
Pass signal vector x
throught a bank of band-pass filters, given sampling rate sr
and bandwidth
parameter.
The kind of filters is specified by optional keyword argument filtkind
, of type FilterDesign, using the DSP package. By default filtkind
is a forward-backward (linear phase response) Butterworth IIR filter of order 2 in each direction (hence, 4th order total). See notes on DSP package useful functions for tips on how to design other IIR and FIR filters.
Return 2-tuple (f, Y)
, where f
is a vector holding the center frequencies of the filter bank band-pass regions and Y
a matrix holding in the $i^{th}$ column the signal x
band-pass iltered by the $i^{th}$ band-pass filter. Hence, size(Y, 1)=length(x)
and size(Y, 2)=length(f)
.
The filter bank is designed by means of argument bandwidth
and optional keyword arguments fmin
and fmax
. These three arguments can be passed either as integers or real numbers. All band-pass regions have bandwidth equal to the bandwidth
argument and overlap with adjacent band-pass regions. By default the lower limit of the first band-pass region is set at bandwidth
Hz, and successive band-pass regions are defined up to, but excluding, the Nyquist frequency ($sr/2$). If fmin
is specified (in Hz), the center frequency of the first band-pass region is set as close as possible, but not below, fmin
. If fmax
is specified (in Hz), the center frequency of the last band-pass region is set as close as possible, but not above, fmax
.
Here are some examples of filter bank definitions given different arguments (sr
=16 in all examples).
bandwidth | fmin | fmax | center frequencies | band-pass regions |
---|---|---|---|---|
4 | - | - | 4 | 2-6 |
2 | - | - | 2, 3, 4, 5, 6 | 1-3, 2-4, 3-5, 4-6, 5-7 |
2 | 3 | 7 | 3, 4, 5, 6 | 2-4, 3-5, 4-6, 5-7 |
1 | 3 | 5 | 3, 3.5, 4, 4.5, 5 | 2.5-3.5, 3-4, 3.5-4.5, 4-5, 4.5-5.5 |
1.1 | 3 | 5 | 2.75, 3.3, 3.85, 4.4, 4.95 | 2.2-3.3, 2.75-8.85, 3.3-4.4, 3.85-4.95, 4.4-5.5 |
1.9 | 3 | 5 | 2.85, 3.8, 4.75 | 1.9-3.8, 2.85-4.75, 3.8-5.7 |
At least sr
samples should be trimmed at the beginning and end of the output signal Y
, as those samples are severely distorted by the filtering process.
If keyword optional argument ⏩ is true (default), the filtering is multi-threaded across band-pass filters. See Threads.
This function is called by the following functions operating on time-frequency reprsentations: TFanalyticsignal
, TFamplitude
, TFphase
, meanAmplitude
, concentration
, meanDirection
, comodulation
, coherence
.
Examples:
using FourierAnalysis, DSP, Plots
# generate a sinusoidal + noise
f, sr, t = 8, 128, 512
v=sinusoidal(1., f, sr, t, 0)
x=v+randn(t)
flabels, Y=filterbank(x, 128)
flabels, Y=filterbank(x, 128; fmin=4, fmax=32)
flabels, Y=filterbank(x, 128, 4; fmin=4, fmax=32)
flabels, Y=filterbank(x, 128, 4;
filtkind=Chebyshev2(8, 10),
fmin=4,
fmax=16)
# trick for plotting the signal filtered in the band-pass regions
for i=1:size(Y, 2) Y[:, i].+=convert(eltype(Y), i)*1.5 end
mylabels=Array{String}(undef, 1, length(flabels))
for i=1:length(flabels) mylabels[1, i]=string(flabels[i])*" Hz" end
plot(Y; c=:grey, labels=mylabels)
notes on DSP package useful functions
using DSP, Plots
### create a time-series with dominant frequency at 10Hz
f, sr, t = 10, 128, 512
v=sinusoidal(1., f, sr, t, 0)+randn(512)
plot(v)
### IIR filters
filterBP = digitalfilter(Bandpass(8/(sr/2), 12/(sr/2)), Chebyshev2(8, 10))
filterLP = digitalfilter(Lowpass(13/(sr/2)), Butterworth(4))
filterHP = digitalfilter(Highpass(7/(sr/2)), Butterworth(4))
# forward filter, unlinear phase response
z=filt(filterBP, v)
plot([v, z])
# forward-backward filter, linear phase response
z=filtfilt(filterBP, v)
plot([v, z])
### FIR filters
responsetype = Bandpass(8, 12; fs=sr)
filtkind = FIRWindow(hanning(64))
filterBP = digitalfilter(responsetype, filtkind)
# forward filter, unlinear phase response
z=filt(digitalfilter(responsetype, filtkind), v)
plot([v, z])
# forward-backward filter, linear phase response
z=filtfilt(digitalfilter(responsetype, filtkind), v)
plot([v, z])
### resampling
z1=resample(x, 0.5) # downsample by 2
z2=resample(x, 2) # upsample by 2
plot([x, z1, z2])
### estimate peak frequency using ESPRIT algorithm
peakv=esprit(v, 200, 2, sr)[1]
# this may fails.
# If the approximate location of the peak is sought,
first band-pass filter the data
filter = digitalfilter(Bandpass(7/(sr/2), 13/(sr/2)), Chebyshev2(8, 10))
z = filtfilt(filter, v)
peakz=esprit(z, 200, 2, sr)[1]
### mean power frequency
meanfreq(v, t)
### root mean square
rms(v)
### cross-correlation
y=sinusoidal(1., f, sr, t, 0)+randn(512)
plot(xcorr(v, y))
### convolution (using the FFT algorithm)
plot(conv(v, y))
### closest product of 2, 3, 5, and 7 greater than or equal to n.
# This is used to find a window length for which FFTW will be able
# to find an efficient FFT plan
nextfastfft(1031)
### estimate the delay by locating the peak of the cross-correlation.
# The output delay will be positive when v is delayed with respect y,
# negative if advanced, 0 otherwise.
finddelay(v, y)
finddelay([0, 0, 1, 2, 3], [1, 2, 3]) # return 2
finddelay([1, 2, 3], [0, 0, 1, 2, 3]) # return -2
### shift elements of signal v in time by a given amount s of samples
### and fill the spaces with zeros.
# For circular shifting, use circshift.
shiftsignal([1, 2, 3], 2) # return [0, 0, 1]
shiftsignal([1, 2, 3], -2) # return [3, 0, 0]
### use finddelay() and shiftsignal() to time align v to y.
# Also return the delay of v with respect to y.
alignsignals([0, 0, 1, 2, 3], [1, 2, 3]) # return ([1, 2, 3, 0, 0], 2)
alignsignals([1, 2, 3], [0, 0, 1, 2, 3]) # return ([0, 0, 1], -2)