NeuroAnalyzer tutorials: Process: filtering

Initialize NeuroAnalyzer 
using DSP
using NeuroAnalyzer
using Plots
eeg = load("files/eeg.hdf");

Filter as a system

Filter (h): a system that changes input signal into output signal.

Temporal filtering: a form of signal processing applied to data time series that attenuates signals at specific frequencies while allowing signal at other frequencies to pass unaffected.

Temporal filtering can introduce distortions of the data at edges of the time series, known as edge artifacts (therefore filters should be applied before epoching). Tapering the signal before filtering may reduce edge artifacts.

Superposition: if the input signal is the sum of some signals, then filtered output signal is the sum of filtered individual signal making the input signal. \[ x = x_1 + x-2 + x_3 \]

\[ h(x) = y \]

\[ y = h(x_1) + h(x_2) + h(x_3) \]

Superposition: signals \(x_1\), \(x_2\), \(x_n\) sums as \(x\), filtered via \(h\) gives \[ x_1 + x_2 + \ldots + x_n \rightarrow h \rightarrow y_1 + y_2 + \ldots + y_n \]

\[ x_1 + x_2 + \ldots + x_n = x \]

\[ y_1 + y_2 + \ldots + y_n = y \]

\[ x \rightarrow h \rightarrow y \]

Filter types

  1. finite impulse response (FIR)
  2. infinite impulse response (IIR)
  3. frequency domain filtering
  4. adaptive filters: Wiener and Kalman filters

Adaptive filter: changes its properties over time in response to signal.

FIR filters properties:

  • filter output is a linear combination only of the input signal at the current time and past times
  • its output necessarily becomes zero within a finite amount of time after the input signal goes to zero
  • impulse response is of finite duration, because it settles to zero in finite time
  • can have linear phase response
  • constant group delay for all frequencies
  • requires higher order (longer kernel), usually 20-2000 = slower to compute
  • no feedback (no poles, just zeros) = always stable
  • BP filter kernel looks like like a complex Morlet wavelet
  • frequency representation (power spectrum) of the kernel is almost rectangular, plateau shaped

FIR filter express each output sample as a weighted sum of the last N input samples, where N is the order of the filter.

If the FIR coefficients are symmetrical (often the case), then such a filter is linear phase.

IIR filters properties:

  • filter output is a linear combination of both the input signal at the current time and past times (feed-forward data flow) and the output signal at past times (feedback data flow)
  • its output may persist indefinitely in the absence of any further input, because the output signal itself is fed back into the filter (therefore it may become unstable)
  • there is an infinite number of impulse response samples, therefore may continue to respond indefinitely (usually decaying)
  • never fully “settles down” after turning off the input, gets quieter and quieter though
  • sometimes unstable, depends on the signal
  • lower order (shorter kernel), usually 4-20 = faster to compute
  • nonlinear phase delay
  • group delay is frequency dependent
  • have feedback loops that may become unstable and oscillate (have poles and zeros)
  • have lower side lobes in the stop-band
  • more ringing on glitches
  • may cause noise buildup, because noise terms created by arithmetic round-off errors are fed back into the filter and amplified
  • IIR filter is steeper than FIR of the same order
  • easier to design, faster to compute (may be used in real-time applications)

Length of the FIR filter response is finite because the response of an ideal filter is multiplied by a window in the time domain, so left and rights parts are zeroed.

t = linspace(-1, 1, 1000)
s = generate_sinc(t, f=10)
s ./= maximum(s)
w = generate_window(:hann, 1000)
plot_signal(t, s)
Plots.plot!(t, w, legend=false)

Causal filter: linear and time-invariant system.

Causal means that impulse response has no values before time \(t = 0\), meaning that the filter output depends only on past and present inputs.

Non-causal filter: a filter whose output depends only on future inputs is anti-causal.

Non-causal filter may be designed to have zero delay (not introduce group delay/phase distortions).

Example: moving average filter.

Real-time filter must be causal.

Causal filter introduces group delay. Initially the convolution of the signal with the impulse response only includes part of the signal. The number of samples at the begging of the signal (when convolution with the impulse response starts) until it includes a full signal is group delay.

For a symmetric FRI filter its: \[ \tau_g = \frac{M - 1}{2} \] \(M\): length of the filter window, in samples

Group delays is a negative derivative of phase, so constant group delays means the filter is linear phase filter.

Group/phase delay: functions that describe the delay times experienced by a signal’s various sinusoidal frequency components as they pass through a filter. These delays are sometimes frequency dependent (different sinusoid frequency components experience different time delays). Delays result in distortions.

Phase delay is in units of time and equals the negative of the phase shift at each frequency divided by the value of that frequency.

Linear phase filter: has a flat phase delay property, meaning that all frequencies are delayed in the same way. Non-flat phase delay causes time delay differences among the signal’s various sinusoidal frequency components leading to distortions.

Group delay is the negative derivative of phase shift with respect to frequency.

Group delay is a convenient measure of the linearity of the phase with respect to frequency in a modulation system.

A stable filter assures that every limited input signal produces a limited filter response. FIR filters are always stable.

Introduction to Z transform

Z transform is used to design complex filters, e.g. for audio signal parametric equalization.

Z transform is applied on a discrete signal and transforms it into Z domain (time \(\rightarrow\) Z domain): \[ X(z) = \sum_{n=-\infty}^{\infty}x[n] \times z^{-n} \] \(n\): sample number

\(z\): variable z

In discrete-time Fourier transform (time \(\rightarrow\) frequency): \[ X(k) = \sum_{n=-\infty}^{\infty}x[n] \times e^{-1i k n} \]

We add combinations of signal and complex sine (cos + sine) waves of various frequencies to find the best correlation between them. When the signal and CSW is not similar they produce product with positive and negative values, which when summed gives low correlation coefficient. CSW is used since signal and real-valued cosine wave could be not in phase and still, while having similar frequency, produce low score. When CSW is used, we get two components: real (correlation with cos) and imaginary (correlation with sine) and from the Euler’s formula: \(e^{1i \times \theta} = \cos(\theta) + 1i \times \sin(\theta)\).

Fourier transform gives information about the signal. However, for filters we need information about the system that produces a signal and how to design the system to get an output we need.

Z transform: \[ X(z) = \sum_{n=-\infty}^{\infty}x[n] \times e^{-1i \omega n} \times r^{-n} \]

This way it gives a correlation to sinusoids (\(e^{-1i \omega n}\)) and exponentials (\(r^{-n}\)) (\(n>1\): growth, \(n<1\): decay). \[ e^{-1i \omega n} \times r^{-n} = (re^{i \omega})^{-n} \]

\[ z = r \times e^{i \omega} \]

and thus \[ X(z) = \sum_{n=-\infty}^{\infty}x[n] \times z^{-n} \]

By knowing the frequency and exponential components of the output system of the response, given a known input, then we know exactly the system the generated that output.

When doing Z transform, we are looking for exponentials and check various values of \(r\) (probing exponential). At some \(r\) value, the probing signal growth at the same rate that the time signal decays. At the value the product is infinite.

Pole: threshold \(r\) value for which the summation reaches infinity. If poles are outside the unit circle, the filter becomes unstable.

Below that value exponential is raising faster than the system is decaying, thus the product is \(> \infty\). The pole is the moment when if jumps from convergence to infinity.

Zero: \(r\) value for which the summation is 0.

If pole is at the value \(r < 1\), then the system is stable (energy of \(x[n]\) exponentially decays over time).

For a finite signal, the summation will never be infinite. Therefore, the Z transform assumes that the signal is infinite. If it is not, we assume what the infinite rest of the data would be.

Z transform can be applied on a model, not on a real signal since it required infinite input. Therefore, there is no fft() equivalent for Z transform. It can be used symbolically: \[ Z(\sin(x)) = \frac{z \sin(1)}{z^2 - 2 \cos(1)z + 1} \]

Z transform provides a tool for analyzing stability issues of digital IIR filters (it is analogous to the Laplace transform, which is used to design and analyze analog IIR filters).

Filter: \[ y(n) = s(n) \times h(n) \] \(\times\): convolution

\(s(n)\): signal

\(h(n)\): impulse response

\[ Y(z) = S(z) \times H(z) \] \(S\): Z transform of \(s\) \(H\): Z transform of \(h\) \(Y\): Z transform of \(y\)

Designing a filter is all about finding an appropriate h(i), which, under the Z transform, reduces to finding an appropriate H(z). Once H(z) is found, h(i) can be obtained from the inverse Z transform.

For continuous signal the Laplace transform converts it into continuous S domain.

Convert from Z to S domain: \[ z = e^{st} \]

Discrete Z domain is the result of converting a discrete signal using the Z transform:

Pole-zero map is plotted in the Z domain:

For IIR filter:

fs = 100
s = rand(1000)
f = filter_create(fprototype=:butterworth, ftype=:bp, order=8, cutoff=(10, 20), fs=fs)
ZeroPoleGain{:z, ComplexF64, ComplexF64, Float64}(ComplexF64[1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, -1.0 + 0.0im, -1.0 + 0.0im, -1.0 + 0.0im, -1.0 + 0.0im, -1.0 + 0.0im, -1.0 + 0.0im, -1.0 + 0.0im, -1.0 + 0.0im], ComplexF64[0.772730222120692 - 0.5642337926102768im, 0.29386652993967477 + 0.8838716409781953im, 0.772730222120692 + 0.5642337926102768im, 0.29386652993967477 - 0.8838716409781953im, 0.695601902379151 - 0.5300448742165345im, 0.2993479526642546 + 0.7577390080587313im, 0.695601902379151 + 0.5300448742165345im, 0.2993479526642546 - 0.7577390080587313im, 0.6099410063393964 - 0.5137953559736712im, 0.3439507174518397 + 0.6493627863090522im, 0.6099410063393964 + 0.5137953559736712im, 0.3439507174518397 - 0.6493627863090522im, 0.5150237094980352 - 0.5227595531392271im, 0.4201703270003071 + 0.5677142673254709im, 0.5150237094980352 + 0.5227595531392271im, 0.4201703270003071 - 0.5677142673254709im], 2.3959644103776224e-5)

Poles:

f.p
16-element Vector{ComplexF64}:
   0.772730222120692 - 0.5642337926102768im
 0.29386652993967477 + 0.8838716409781953im
   0.772730222120692 + 0.5642337926102768im
 0.29386652993967477 - 0.8838716409781953im
   0.695601902379151 - 0.5300448742165345im
  0.2993479526642546 + 0.7577390080587313im
   0.695601902379151 + 0.5300448742165345im
  0.2993479526642546 - 0.7577390080587313im
  0.6099410063393964 - 0.5137953559736712im
  0.3439507174518397 + 0.6493627863090522im
  0.6099410063393964 + 0.5137953559736712im
  0.3439507174518397 - 0.6493627863090522im
  0.5150237094980352 - 0.5227595531392271im
  0.4201703270003071 + 0.5677142673254709im
  0.5150237094980352 + 0.5227595531392271im
  0.4201703270003071 - 0.5677142673254709im

Zeros:

f.z
16-element Vector{ComplexF64}:
  1.0 + 0.0im
  1.0 + 0.0im
  1.0 + 0.0im
  1.0 + 0.0im
  1.0 + 0.0im
  1.0 + 0.0im
  1.0 + 0.0im
  1.0 + 0.0im
 -1.0 + 0.0im
 -1.0 + 0.0im
 -1.0 + 0.0im
 -1.0 + 0.0im
 -1.0 + 0.0im
 -1.0 + 0.0im
 -1.0 + 0.0im
 -1.0 + 0.0im

Draw pole-zero map:

plot_polezero(f.p, f.z)

Gain:

f.k
2.3959644103776224e-5

FIR filters contain as many poles as they have zeros, but all of the poles are located at the origin, \(z=0\). Because all of the poles are located inside the unit circle, the FIR filter is always stable.

Filter design

Filter filters a signal at a cutoff frequency.

A low pass (LP) filter is used to cut unwanted high-frequency signals, above the cutoff.

A high pass (HP) filter passes high frequencies, cuts below the cutoff.

HP is all pass filter minus LP.

A band pass (BP) filter passes a limited range of frequencies within the cutoff.

A band stop (BS) filter passes frequencies above and below the cutoff range. A very narrow band stop filter is known as a notch filter.

Window design method: required frequency response of a filter response is converted into a filter kernel.

Ideal frequency response:

Next, an impulse response (e.g. sinc) is generated from the ideal filter.

The infinite impulse is multiplied by a window function. Multiplying by window function in the time domain results in the frequency response of the IIR being convolved with the Fourier transform of the window function.

The result of the frequency domain convolution is that the edges of the rectangle are tapered, and ripples appear in the pass-band and stop-band. Working backward, one can specify the slope (or width) of the tapered region (transition band) and the height of the ripples, and thereby derive the frequency-domain parameters of an appropriate window function.

Left: Hamming window, right: rectangular window.

The trick is to select the window type and filter length that will give a filter with the correct rate of roll-off and level of attenuation in the stop band.

To avoid ripple artifacts (at the cutoff frequency) in the time domain the filter should be designed to have a smooth transition in the frequency domain.

Impulse response: filter output when the input signal is the Kronecker delta function (single impulse at time 0 and zeros everywhere else).

Filter -3 dB = attenuation by 30% = 70% of the signal is passed through.

Cutoff: half-amplitude frequency = frequency that is suppressed in 50%, e.g. 25 Hz LP filter: passes 10% of the signal at 50 Hz.

Order: roll off of a digital filter = number of points that are averaged by the filter at any one time; higher order means steeper roll-off.

For FIR filters: filter order (taps) is the number of samples used in calculation; more samples means more accuracy (better frequency response), but is slower.

Pass-band: frequency at which the output is preserved.

Stop-band: frequency at which the output is severely attenuated.

The designed filter is translated into filter kernel, that gives more or less imperfect frequency response:

Frequency response of a filter is the amplitude spectrum of the filter kernel.

Frequency response of a low pass filter:

-60 dB is the filter attenuation.

In NeuroAnalyzer transition band is calculated as:

  • for LP filter: f_pass = cutoff[1] - (bw / 2), f_stop = cutoff[1] + (bw / 2)
  • for HP filter: f_pass = cutoff[1] + (bw / 2), f_stop = cutoff[1] - (bw / 2)
  • for BP filter: f1_stop = cutoff[1] - (bw / 2), f1_pass = cutoff[1] + (bw / 2), f2_pass = cutoff[2] - (bw / 2), f2_stop = cutoff[2] + (bw / 2)
  • for BS filter: f1_pass = cutoff[1] - (bw / 2), f1_stop = cutoff[1] + (bw / 2), f2_stop = cutoff[2] - (bw / 2), f2_pass = cutoff[2] + (bw / 2)

where bw is the transition band width in Hz and cutoff specifies cutoff frequency for LP and HP or frequencies for BP and BS.

Wider plateau means lower frequency precision and higher temporal precision.

Wider transition zone means less pass band ripple = cleaner, better filter kernel.

Narrower transition zone mean steeper filter, but also less stop band attenuation (higher amplitude of ripples in the stop band).

Left: 31 taps, right: 11 taps.

Example: impulse response of FIR LP filter, cut off at 10 Hz, default Hamming window:

flt = filter_create(fprototype=:fir, ftype=:lp, cutoff=10, order=91, fs=256)
t = linspace(-1, 1, length(flt))
plot_signal(t, flt, ylims=(-0.11, 0.11))

And its frequency response:

f, _ = freqs(flt, 256)
_, _, p, _ = ftransform(flt)
plot_psd(f, p, frq_lim=(0, 40))

Detailed frequency response:

plot_filter_response(fprototype=:fir, ftype=:lp, cutoff=10, order=91, fs=256)

Same FIR band pass filter, but using Bartlett window:

flt = filter_create(fprototype=:fir, ftype=:lp, cutoff=10, order=91, fs=256, w=DSP.bartlett(91))
t = linspace(-1, 1, length(flt))
plot_signal(t, flt, ylims=(-0.11, 0.11))

And its frequency response:

f, _ = freqs(flt, 256)
_, _, p, _ = ftransform(flt)
plot_psd(f, p, frq_lim=(0, 40))

Detailed frequency response:

plot_filter_response(fprototype=:fir, ftype=:lp, cutoff=10, order=91, fs=256, w=DSP.bartlett(91))

Other impulse responses:

Low pass FIR filter:

flt = filter_create(fprototype=:fir, ftype=:lp, cutoff=10, order=91, fs=256)
t = linspace(-1, 1, length(flt))
plot_signal(t, flt, ylims=(-0.11, 0.11))

High pass FIR filter:

flt = filter_create(fprototype=:fir, ftype=:hp, cutoff=10, order=91, fs=256)
Plots.plot!(t, flt, ylims=(-0.11, 0.11), label=false)

Therefore, when designing a filter, check its amplitude spectrum over the ideal design.

When FIR filter is applied on the signal, its frequency response modifies individual frequencies of the input signal:

Bode plot: graphical representation that shows how the gain (magnitude) and phase of a system respond over a range of frequencies.

Applying a filter introduces a phase shift in the data. The amount of the phase shift is determined by the length of the filter kernel. These phase shifts are an artifact and should be corrected.

Phase shift correction:

  • filter the data
  • flip the time series backwards
  • apply the same filter again
  • flip the data forward again

Phase distortions introduced by filtering are reintroduced in reversed-time, thus reverse-distorting the distortions.

FFT filtering

General rule: signal → FFT → signal .* kernel → IFFT

Constraints:

  1. accurate for FIR filters only
  2. for very large N use block convolution
  3. zero-pad signal to avoid time-domain aliasing

Perform filtering in short windows (1 second) (= block convolution).

Use windowing to remove edge artifacts.

To avoid loss of information due to windowing use overlapping (by 25-75%) epochs.

HP kernel (at 7 Hz): 1 / (1 + exp(-frq + 7)), where frq = linspace(0, fs, n) and n = length(signal)

LP kernel (at 7 hz): 1 - (1 / (1 + exp(-frq + 7)))

BP kernel: plateau-shaped, e.g. rectangular

Order: number of time points of kernel (length(kernel) = order + 1)

For EEG signals, a useful rule of thumb is that order should be 4 to 5 times the number of time points in one cycle of the lowest analyzed frequency. For example, in gamma rhythm (35-45 Hz) one cycle is about 25 ms and the order of the filter would be be 100 ms. If the sampling rate is 500 Hz, then the filter order would be 50 data points.

Practical aspects of EEG filtering

HP 0.05-1.0 Hz: removes slow artifacts, reduces offset/slow trends (demeaning/detrending)

LP 35.0-70.0 Hz: removes fast artifacts

Notch (BS) 50.0 Hz: removes AC noise

HP filter should always be applied to continuous, non-epoched data. Filtering at 0.5 Hz gives edge artifacts lasting up to 6 seconds.

HP makes analysis of very slow oscillations difficult. Solution: use DC-coupled amplifiers. For very fast (>100 Hz) oscillations MEG is better.

For analyzing a time-domain signal in a specified frequency range: FIR = IIR.

For analyzing phase values of the filtered signal: FIR > IIR.

NeuroAnalyzer filters

The following filters (specified using the fprototype parameter) are available:

FIR filters:

  • window method (:fir)
  • weighted least-squares design (:firls)
  • Remez (:remez)

IIR filters:

  • Butterworth (:butterworth)
  • Chebyshev1 (:chebyshev1)
  • Chebyshev2 (:chebyshev2)
  • Elliptic (:elliptic)
  • second-order IIR notch filter (:iirnotch)

All major filter types (specified using the ftype parameter) are available: low pass (:lp), high pass (:hp), band pass (:bp) and band stop (:bs).

Filter cutoff must be specified in Hz (cutoff parameter). For band pass and band stop filters, a pair of frequencies must be provided.

Filter order is specified usign the order parameter: order is the filter steepness for :butterworth, :chebyshev1, :chebyshev2, :elliptic or number of taps for :fir, :firls and :remez filters

Filtering includes two steps:

  1. filter preparation (filter_create())
  2. filtering (filter_apply())

The filter() function combines both steps.

Filtering direction is set using the dir parameter:

  • :twopass: two passes, the resulting signal has zero phase distortion, the effective filter order is doubled
  • :onepass: single pass
  • :reverse: one pass, reverse direction

Also, smoothing (LP) filters are available:

  • moving average
  • moving median
  • polynomial
  • Savitzky-Golay

FFT filtering is performed using the Gaussian filter.

FIR filters

Order is equal to the number of filter taps.

Filter order can be estimated using F. Harris formula: \[ N = \frac{fs}{bw} \frac{A}{22} \] \(fs\): sampling rate

\(bw\): transition band width

\(A\): attenuation, in dB

n = fir_order_bw(eeg, bw=0.2, a=50)
2909

(!) bw is transition band width in Hz and a is attenuation in dB.

(!) This estimation is just the entry point for further investigation. Filter response should be analyzed using plot_filter_response() or ifilter().

Default window is Hamming. Its length is equal to the filter order.

Custom windows are specified using w parameter. If both window and order are specified, the window length must be equal to the filter order.

Order or custom window length must be odd for HP, BP and BS filters.

FIR HP filter at 0.5 Hz, transition band width of 0.2 Hz and order of 2901, with Hamming window:

NeuroAnalyzer.filter(eeg,
                     ch="all",
                     fprototype=:fir,
                     ftype=:hp,
                     cutoff=0.5,
                     order=2901,
                     bw=0.2);

Same filter with custom Hanning window:

NeuroAnalyzer.filter(eeg,
                     ch="all",
                     fprototype=:fir,
                     ftype=:hp,
                     cutoff=0.5,
                     w=DSP.hanning(2901),
                     bw=0.2);

(!) Certain windows (e.g. hanning()) may require DSP.jl to be loaded. NeuroAnalyzer allows generating additional windows, see generate_window().

FIR weighted least-squares design

Transition bandwidth must be provided (bw) in Hz. Based on bw value, f_pass and f_stop are calculated.

Order is equal to the number of filter taps.

Filter weights may be provided (w) (4 number for LP and HP, 6 numbers for BP and BS filters). This values modifies the impact of the filter on individual filter’s frequency response points.

Filter order can be estimated using Harris formula:

n = fir_order_bw(eeg, bw=1.5, a=50)
388

(!) bw is transition band width in Hz and a is attenuation in dB.

(!) This estimation is just the entry point for further investigation. Filter response should be analyzed using plot_filter_response() or ifilter().

Sometimes, FIR filter order needs to be calculated for a specific frequencies studied. A useful rule of thumb is to use 4 to 5 cycles of the lowest frequency studied. For example, when studying the gamma rhythm in the range of 35-45 Hz, one slowest cycle is about 0.028 s (25 ms) (\(1 / 35\)). The order would be from 100 (\(4 \times 25\)) to 125 (\(5 \times 25\)) ms. We need the order to be in samples, so for the sampling rate of 500 this would be 50 t2s(0.1, 500) to 63 t2s(0.125, 500) samples. This can be also calculated using:

fir_order_f(eeg, f=35)
(32, 40)

(!) Our sampling rate is 256 here.

FIR LP filter at 40.0 Hz, transition band width of 1.5 Hz and order of 200:

NeuroAnalyzer.filter(eeg,
                     ch="all",
                     fprototype=:firls,
                     ftype=:lp,
                     cutoff=40,
                     order=200,
                     bw=1.5);

Filter kernel:

flt = filter_create(fprototype=:firls, ftype=:lp, cutoff=40, order=200, fs=256, bw=1.5)
t = linspace(-1, 1, length(flt))
plot_signal(t, flt, ylims=(-0.5, 0.5))

FIR Remez filter

Calculate the minimax optimal filter using the Remez exchange algorithm.

Transition bandwidth must be provided (bw) in Hz. Based on bw value, f_pass and f_stop are calculated.

Order is equal to the number of filter taps.

Filter weights may be provided (w) (4 number for LP and HP, 6 numbers for BP and BS filters). This values modifies the impact of the filter on individual filter’s frequency response points.

Filter order can be estimated using harris formula:

n = fir_order_bw(eeg, bw=1.5, a=50)
388

(!) bw is transition band width in Hz and a is attenuation in dB.

(!) This estimation is just the entry point for further investigation. Filter response should be analyzed using plot_filter_response() or ifilter().

FIR BP filter at 8.0 to 14.0 Hz, transition band width of 0.5 Hz and order of 200:

NeuroAnalyzer.filter(eeg,
                     ch="all",
                     fprototype=:remez,
                     ftype=:bp,
                     order=200,
                     cutoff=(8, 14),
                     bw=0.5);

Filter kernel:

flt = filter_create(fprototype=:remez, ftype=:bp, cutoff=(8, 14), order=200, fs=256, bw=0.5)
t = linspace(-1, 1, length(flt))
plot_signal(t, flt, ylims=(-0.5, 0.5))

IIR filters

Filter response type:

cutoff = 40
ftype = DSP.Lowpass(cutoff)
Lowpass{Float64}(40.0)

Filter prototype:

order = 8
fprototype = DSP.Butterworth(order)
ZeroPoleGain{:s, ComplexF64, ComplexF64, Float64}(ComplexF64[], ComplexF64[-0.19509032201612828 + 0.9807852804032304im, -0.19509032201612828 - 0.9807852804032304im, -0.5555702330196022 + 0.8314696123025452im, -0.5555702330196022 - 0.8314696123025452im, -0.8314696123025452 + 0.5555702330196022im, -0.8314696123025452 - 0.5555702330196022im, -0.9807852804032304 + 0.19509032201612828im, -0.9807852804032304 - 0.19509032201612828im], 1.0)
plot_polezero(fprototype.p, fprototype.z)

Available filter prototypes are :butterworth, :elliptic, :chebyshev1 and :chebyshev2:

Chebyshev filter has sharper transition from the bandpass to the bandstop than Butterworth.

IIR filter order can be estimated using:

n = iir_order(eeg, fprototype=:butterworth, ftype=:lp, cutoff=35, bw=1.5)
53

(!) This estimation is just the entry point for further investigation. Filter response should be analyzed using plot_filter_response() or ifilter().

IIR Butterworth LP filter at 45 Hz and order of 53:

NeuroAnalyzer.filter(eeg,
                     ch="all",
                     fprototype=:butterworth,
                     ftype=:lp,
                     cutoff=35.0,
                     order=53);

(!) For :butterworth, :elliptic, :chebyshev1 and :chebyshev2 we do not specify transition band width.

Second-order IIR notch filter:

NeuroAnalyzer.filter(eeg,
                     ch="all",
                     fprototype=:iirnotch,
                     cutoff=50,
                     order=50,
                     bw=2);

Moving average filter

Moving average filter (MAF): simplest low pass filter.

MAF averages d samples before and after the current sample: f = mean(x(n - d) : x(n + d)).

y1 = rand(100)
t = range(0, 1, 100)

d = 3
x = ones(d) ./ d
y2 = DSP.filt(x, y1)

plot_signal(t, y1)
Plots.plot!(t, y2, label=false)

\(d\): must be chosen empirically, based on results

MAF works well when the noise is random and with small amplitude relative to the signal.

MAF order: window length (number of averaged points) (\(2 \times d + 1\)).

Zeros of MAF (frequencies at which amplitude becomes 0) and phase changes are at:

sampling rate / order
2 × sampling rate / order
3 × sampling rate / order
4 × sampling rate / order

For example for order 8 and 1000 Hz: \[ 1000/8 = 125\ \mathrm{Hz} \]

\[ 1000/4 = 250\ \mathrm{Hz} \]

\[ 3 \times 1000 / 8 = 375\ \mathrm{Hz} \]

\[ 1000/2 = 500\ \mathrm{Hz} \]

y1 = rand(100)
t = range(0, 1, 100)

y2 = filter_mavg(y1, k=2) # order is 2 * k + 1

plot_signal(t, y1)
Plots.plot!(t, y2)

If window length is \(2 \times d + 1\) then for cutoff frequency \(f\), d is \(\sqrt{0.196202 + f^2} / f\).

If window length is \(2 \times d + 1\) then for d, the approximate cutoff frequency: \(0.442947 / \sqrt{(2 * d + 1)^2 - 1} * fs\).

NeuroAnalyzer.filter_mavg(eeg,
                          ch="all",
                          k=2,
                          ww=[1, 2, 3, 2, 1]);

(!) The weighting window parameter ww=[1, 2, 3, 2, 1] will produce a triangular-window moving average.

Moving median filter

Moving median filter (MMF): similar to MAF: f = median(x(n - d) : x(n + d)).

MMF is particularly useful if there are a few noise spikes of very large amplitude.

MMF may be applied several times, each time using a different d parameter, this is useful for time series with considerable noise.

After MMF is computed several times, the resulting median-filtered signal can be averaged together.

y1 = rand(100)
t = range(0, 1, 100)

y2 = filter_mmed(y1, k=2) # order is 2 * k + 1

plot_signal(t, y1)
Plots.plot!(t, y2)
NeuroAnalyzer.filter_mmed(eeg,
                          ch="all");

Polynomial filter

Polynomial fitting filter

Polynomial order should not be set higher than necessary.

Models with a large number of parameters require a lot of data for reliable parameter estimation.

Polynomial order must roughly match the changes in the signal; under-fitting and over-fitting will produce suboptimal results by being insensitive to the signal changes or by being too sensitive to noise.

Polynomial fitting can be a useful de-noising strategy when the signal and noise can be distinguished by rapidly vs. slowly changing time series features, even if they do not have well-defined frequency characteristics (e.g. slow signal, fast noise or vice versa).

Filter using a polynomial of the 8th order:

NeuroAnalyzer.filter_poly(eeg,
                          ch="all",
                          order=8);

Savitzky-Golay filter

Savitzky–Golay filter: digital filter that can be applied to a set of digital data points for the purpose of smoothing the data, that is, to increase the precision of the data without distorting the signal tendency.

This is achieved, in a process known as convolution, by fitting successive sub-sets of adjacent data points with a low-degree polynomial by the method of linear least squares.

NeuroAnalyzer.filter_sg(eeg,
                        ch="all",
                        order=6);

Gaussian filter

Gaussian filters extracts a specific frequency component of FTT:

and computes IFFT of this part of FTT.

fs = 256
t = collect(0:(1 / fs):10)
a = [5, 5]
f = [2, 5]
s = zeros(length(t))
for idx in 1:length(a)
    s = s .+ generate_sine(f[idx], t, a[idx])
end
plot_signal(t, s)

Band pass Gaussian filter at 2 Hz:

s_new = filter_g(s, fs=fs, f=2)
plot_signal(t, s_new)

Plot filter response

Filter response can be plotted using the preview parameter:

NeuroAnalyzer.filter(eeg,
                     ch="eeg",
                     fprototype=:fir,
                     ftype=:lp,
                     cutoff=25,
                     order=15,
                     preview=true)

(!) When preview=true, signal will not be filtered.

Another way is to use plot_filter_response():

plot_filter_response(fs=sr(eeg),
                     fprototype=:butterworth,
                     ftype=:bs,
                     cutoff=(45, 55),
                     order=8)

Filter designer

You may also use the interactive filter designer to design the filter:

flt = ifilter(eeg)

and later apply the designed filter:

filter_apply(eeg,
             ch="all",
             flt=flt)