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 — Function

function 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

Nota Bene

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)

source

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)