Filtering

Initialize NeuroAnalyzer

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

A filter ($h$) is a system that transforms an input signal ($x$) into an output signal ($y$). Filters are fundamental tools in signal processing, including EEG analysis, where they are used to:

Remove noise or artifacts.
Isolate specific frequency bands (e.g., alpha, beta, gamma).
Enhance signal-to-noise ratio (SNR).

Mathematically:

\[ x \xrightarrow{h} y \]

Temporal Filtering

Definition: A form of signal processing applied to time-series data that attenuates (reduces) signals at specific frequencies while allowing signals at other frequencies to pass unaffected.

Mathematically:

\[ x = x_1 + x_2 + x_3 \]

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

This property ensures that filtering is linear and additive, making it predictable and reliable.

Types of Temporal Filters

Type	Description	Example Use Case
Low-Pass Filter	Allows frequencies below a cutoff to pass; attenuates higher frequencies.	Removing high-frequency noise (e.g., muscle artifacts) in EEG.
High-Pass Filter	Allows frequencies above a cutoff to pass; attenuates lower frequencies.	Removing slow drifts or DC offsets in EEG.
Band-Pass Filter	Allows frequencies within a specific range to pass; attenuates frequencies outside the range.	Isolating alpha (8–12 Hz) or beta (13–30 Hz) bands in EEG.
Band-Stop Filter	Attenuates frequencies within a specific range; allows frequencies outside the range to pass.	Removing power-line noise (50/60 Hz) in EEG.

Edge Artifacts

Cause: Filters assume the signal is stationary (unchanging over time). At the edges of the signal, this assumption breaks down, leading to distortions.
Impact:
- Ringing artifacts: Oscillations at the edges of the signal.
- Phase shifts: Delays in the filtered signal.
Solution:
- Apply filters before epoching (segmenting into trials).
- Use tapers (e.g., Hanning, Hamming, or Tukey windows) to smooth the edges of the signal before filtering.

Tapering

Definition: Multiplying the signal by a window function to reduce discontinuities at the edges.
Common Tapers:
- Hanning Window: Smooths the edges with a cosine taper.
- Hamming Window: Similar to Hanning but with a different taper shape.
- Tukey Window: Combines a flat top with cosine tapers at the edges.
Effect: Reduces edge artifacts but may slightly widen the filter’s transition band.

Superposition Principle in Detail

The superposition principle states that the filtered output of a sum of signals is equal to the sum of the filtered outputs of the individual signals. This is a fundamental property of linear time-invariant (LTI) systems.

Mathematical Formulation

Input Signal:

\[ x(t) = x_1(t) + x_2(t) + \ldots + x_n(t) \]

Filtered Output:

\[ y(t) = h(x(t)) = h(x_1(t)) + h(x_2(t)) + \ldots + h(x_n(t)) \]

Example:

Input: $x(t) = \sin(2\pi \cdot 10t) + \sin(2\pi \cdot 50t)$ (10 Hz + 50 Hz sine waves).
Filter: Band-pass filter centered at 10 Hz.
Output: $y(t) = h((210t)) + h((250t)) = (210t) + 0 (the 50 Hz component is attenuated).

Practical Implications

Linearity:

Filters are linear systems, meaning they obey the superposition principle.
This allows for modular analysis (e.g., filtering individual frequency bands separately).

Time-Invariance:

Filters are time-invariant, meaning their behavior does not change over time.
This ensures that the filtering process is consistent across the entire signal.

Applications:

EEG Analysis: Isolate specific frequency bands (e.g., alpha, beta) for cognitive or clinical studies.
Noise Reduction: Remove artifacts (e.g., muscle activity, power-line noise) from signals.
Signal Enhancement: Improve the signal-to-noise ratio (SNR) for better interpretation.

Common Filtering Techniques

Technique	Description	Use Case
FIR (Finite Impulse Response)	Uses a finite number of samples to compute the output; no feedback.	Real-time filtering, high precision.
IIR (Infinite Impulse Response)	Uses feedback to compute the output; can achieve sharp cutoffs with fewer coefficients.	Low-latency filtering, efficiency.
Butterworth Filter	Maximally flat frequency response in the passband; no ripple.	EEG, audio processing.
Chebyshev Filter	Steeper roll-off than Butterworth but introduces ripple in the passband or stopband.	Applications requiring sharp cutoffs.
Elliptic Filter	Steepest roll-off but introduces ripple in both passband and stopband.	High-performance filtering.

Filter Design Considerations

Cutoff Frequency:

Choose a cutoff frequency that aligns with your analysis goals (e.g., 8 Hz for alpha band in EEG).

Filter Order:

Higher-order filters have sharper cutoffs but may introduce phase distortions or ringing artifacts.
Lower-order filters are smoother but may not adequately attenuate unwanted frequencies.

Phase Response:

Zero-phase filtering: Applies the filter forward and backward to eliminate phase shifts (common in EEG analysis).
Linear-phase filtering: Introduces a constant delay but preserves phase relationships.

Transition Band:

The range of frequencies between the pass-band and stop-band.
A narrower transition band requires a higher-order filter.

Filter Types

Finite impulse response (FIR) filters
Infinite impulse response (IIR) filters
Frequency domain filtering
Adaptive filters: Wiener and Kalman filters

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

FIR Filters

Finite Impulse Response (FIR) filters have a finite duration of impulse response. This means the filter’s output depends only on the current and past inputs (no feedback).

Linear Combination of Inputs: The output of an FIR filter is a linear combination of the current and past input samples (no feedback or past outputs are used).

Mathematically, the output $y[n]$ is a weighted sum of the current and past inputs $x[n], x[n-1], \ldots, x[n-N]$:

\[ y[n] = \sum_{k=0}^N h[k] \cdot x[n-k] \]

where:
$y[n]$ is the output signal at time $n$,
$h[k]$ are the filter coefficients (impulse response),
$x[n−k]$ is the input signal at time $n−k$,
$N$ is the order of the filter (number of coefficients).

Key Characteristics

Feature	Description
Impulse Response	Finite duration (tails off to zero).
Stability	Always stable (no feedback).
Phase Response	Can be designed to have linear phase (no phase distortion).
Computational Cost	Higher computational cost due to more coefficients for sharper cutoffs.
Design Flexibility	Highly flexible; can achieve arbitrary frequency responses.	a

Advantages

No feedback: Stable and easy to implement.
Linear phase: Preserves the shape of the signal (important for EEG and audio).
Robustness: Less sensitive to finite precision arithmetic.

Disadvantages

Higher computational cost: Requires more coefficients for sharp cutoffs.
Longer delay: Introduces a group delay equal to half the filter length.

Applications

EEG/MEG preprocessing: Removing artifacts (e.g., power-line noise, muscle activity).
Real-time systems: Due to stability and linear phase.

Design Methods

Window Method: Multiply the ideal impulse response by a window (e.g., Hamming, Hanning).
Frequency Sampling: Design the filter in the frequency domain.
Least Squares: Minimize the error between the desired and actual frequency response.

No Feedback: FIR filters do not use feedback (no poles in the transfer function). This makes them always stable. There is no risk of instability or oscillations (unlike IIR filters).

Impulse Response:

The impulse response of an FIR filter is the sequence of filter coefficients $h[k]$. Since the filter output is a weighted sum of a finite number of past inputs, the impulse response settles to zero after $N$ samples.

Mathematically, if the input is an impulse $\delta[n]$, the output is:

\[ y[n] = h[n] \quad \text{for} \quad n = 0, 1, \ldots, N-1 \]

and $y[n] = 0$ for $n \geq N$.

Finite Duration:

The impulse response lasts exactly $N$ samples (where $N$ is the filter order). This is why FIR filters are called “finite impulse response”.

Linear Phase:

FIR filters can be designed to have linear phase, meaning all frequency components of the signal are delayed by the same amount.

Linear phase ensures that the shape of the signal is preserved (no distortion).

Achieved by designing symmetric coefficient filters.

Symmetric Coefficients:

If the FIR filter coefficients are symmetric (i.e., $h[k] = h[N-1-k]$), the filter has linear phase. Symmetry simplifies the design and reduces computational complexity.

Example: Symmetric FIR Filter

Coefficients: $h = [1, 2, 3, 2, 1]$
Symmetric around the center coefficient.

Advantages of Symmetric FIR Filters

Linear phase: Preserves signal shape.
Reduced computational cost: Can be implemented with fewer multiplications.

Constant Group Delay:

The group delay (the delay of the envelope of a signal) is constant for all frequencies in a linear-phase FIR filter.

This is critical for applications like EEG and audio processing, where phase integrity is important.

Power Spectrum of the Kernel:

The frequency response of an FIR filter is the Fourier transform of its impulse response $h[k]$.

For a band-pass FIR filter, the frequency response resembles a rectangular or plateau-shaped pass-band (with smooth transitions to the stop-band).

Example: A complex Morlet wavelet kernel is often used for band-pass filtering in EEG analysis.

Higher Order (Longer Kernel):

FIR filters typically require a larger number of coefficients ($N=20 to 2000$) to achieve sharp cutoffs.

This makes them computationally expensive compared to IIR filters for the same cutoff sharpness.

The convolution operation (applying the filter) involves $N \times M$, where $M$ is the length of the input signal.

For real-time applications, optimized implementations (e.g., using FFT-based convolution) are often used.

Example: FIR Filter Order

Filter Type	Typical Order ($N$)	Applications
Low-pass	50–500	Removing high-frequency noise.
High-pass	50–500	Removing DC offsets.
Band-pass	100–2000	Isolating alpha/beta bands in EEG.
Band-stop	50–500	Removing power-line noise (50/60 Hz).

IIR Filters

Infinite Impulse Response (IIR) filters are a powerful class of digital filters that use feedback to achieve sharp frequency responses with fewer coefficients than FIR filters. However, this feedback introduces stability challenges and non-linear phase behavior.

Feedforward and Feedback Data Flow: The output of an IIR filter is a linear combination of:

Current and past inputs (feedforward path).
Past outputs (feedback path).

Mathematically:

\[ y[n] = \sum_{k=0}^{N} b_k \cdot x[n-k] - \sum_{k=1}^{M} a_k \cdot y[n-k] \]

where:
$y[n]$ is the output signal at time $n$,
$x[n−k]$ is the input signal at time $n-k$,
$b_k$ is the feedforward coefficients (zeros),
$a_k$ si the feedback coefficients (poles),
$N$ is the order of the feedforward path,
$M$ is the order of the feedback path.

Feedback Loop:

The output signal is fed back into the filter, creating a closed-loop system.
This feedback can lead to instability if not carefully designed.

Impulse Response:

The impulse response of an IIR filter does not settle to zero after a finite number of samples. Instead, it decays exponentially over time (if stable) or may oscillate indefinitely (if unstable).

Mathematically, if the input is an impulse $\delta[n]$, the output is:

\[ y[n] = h[n] \quad \text{for} \quad n \geq 0 \]

where $h[n]$ is the impulse response, which may never fully reach zero.

Persistent Output:

Even after the input signal goes to zero, the output may persist indefinitely due to feedback.
This is why IIR filters are called “infinite impulse response” filters.

Dependence on Signal:

IIR filters may become unstable depending on the input signal and the choice of coefficients.
Instability occurs if the feedback loop amplifies errors or noise instead of damping them.

Conditions for Stability:

All poles (roots of the denominator polynomial) must lie inside the unit circle in the z-plane.

Mathematically, for a transfer function:

\[ H(z) = \frac{B(z)}{A(z)} = \frac{b_0 + b_1 z^{-1} + \ldots + b_N z^{-N}}{1 + a_1 z^{-1} + \ldots + a_M z^{-M}} \]

the poles (roots of $A(z)$) must satisfy $\vert z \vert < 1$.

Unstable IIR filters can produce oscillations, divergence, or noise buildup.

Lower Order (Shorter Kernel):

IIR filters typically require fewer coefficients ($N=4 to 20$) to achieve sharp cutoffs compared to FIR filters.
This makes them faster to compute and suitable for real-time applications.
The convolution operation for IIR filters involves fewer multiplications than FIR filters of the same order.

Example: IIR Filter Order

Filter Type	Typical Order ($N$)	Applications
Low-pass	4–10	Removing high-frequency noise.
High-pass	4–10	Removing DC offsets.
Band-pass	6–20	Isolating alpha/beta bands in EEG.
Band-stop	4–10	Removing power-line noise (50/60 Hz).

Nonlinear Phase and Group Delay:

IIR filters do not have linear phase by default. This means:

Different frequency components are delayed by different amounts.
The group delay (delay of the envelope of a signal) is frequency-dependent.

Impact on Signal:

Nonlinear phase can distort the shape of the signal, which is problematic for applications like EEG and audio.

Poles and Zeros:

IIR filters are characterized by poles (feedback coefficients) and zeros (feedforward coefficients) in the z-plane.
The location of poles determines stability:
- Poles inside the unit circle: Stable.
- Poles outside the unit circle: Unstable (oscillations or divergence).

For a transfer function:

\[ H(z) = \frac{(z - z_1)(z - z_2)}{(z - p_1)(z - p_2)} \]

where:
$z_1, z_2$ are the zeros (feedforward),
$p_1, p_2$ are the poles (feedback).

If $\vert p_1 \vert > 1$ or $\vert p_2 \vert >1$, the filter is unstable.

Feedback Loops:

Feedback loops can amplify noise or errors due to arithmetic round-off, leading to noise buildup.
Example: In real-time systems, round-off errors in feedback can accumulate over time.

Steep Frequency Response:

IIR filters can achieve sharper cutoffs than FIR filters of the same order.
This is due to the feedback mechanism, which enhances the filter’s selectivity.

Lower Side Lobes in Stopband:

IIR filters typically have lower side lobes in the stopband, meaning they better suppress unwanted frequencies.

Ringing:

IIR filters may produce ringing artifacts (oscillations) when subjected to sudden changes (e.g., glitches, transients).
This is due to the feedback loop, which can cause the filter to “ring” before settling.

Noise Buildup

Noise terms created by arithmetic round-off errors can be fed back into the filter and amplified.
This is particularly problematic in fixed-point arithmetic or real-time systems.
Mitigation:
- Use higher-precision arithmetic (e.g., floating-point instead of fixed-point).
- Design filters with sufficient stability margins.

Easier to Design:

IIR filters are easier to design than FIR filters for achieving sharp cutoffs with fewer coefficients.
Common design methods include Butterworth, Chebyshev, and Elliptic filters.

Faster to Compute:

Due to their lower order, IIR filters are faster to compute and suitable for real-time applications.

IIR Filter Design Methods

Method	Description	Applications
Butterworth	Maximally flat frequency response in the passband.	EEG, audio preprocessing.
Chebyshev	Steeper roll-off but introduces ripple in the passband or stopband.	Narrowband filtering.
Elliptic	Steepest roll-off but introduces ripple in both passband and stopband.	High-performance filtering.

t = linspace(-1, 1, 1000)
s = generate_sinc(t; f = 10)
s ./= maximum(s)
w = generate_window(:hann, 1000)
NeuroAnalyzer.plot(t, s)
NeuroAnalyzer.plot(t, w)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/kJl0u/src/scenes.jl:264

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/kJl0u/src/scenes.jl:264

Causal and Non-Causal Filters

Filters can be classified as causal or non-causal based on their dependence on past, present, and future input samples. This classification is crucial for real-time signal processing and offline analysis.

Causal Filters

A causal filter is a linear and time-invariant system whose output at any time $t$ depends only on the current and past inputs (not on future inputs).

Mathematically, for a causal filter:

\[ y(t) = \int_{-\infty}^t h(t - \tau) x(\tau) d\tau \]

or in discrete time:

\[ y[n] = \sum_{k=0}^{N} h[k] \cdot x[n-k] \]

where:
$y(t)$ or $y[n]$ is the output signal,
$h(t)$ or $h[k]$ is the impulse response (causal if $h(t) = 0$ for $t < 0$),
$x(\tau)$ or $x[n-k]$ is the input signal.

Key Characteristics

Feature	Description
Impulse Response	$h(t) = 0$ for $t < 0$ (no values before $t = 0$).
Output Dependency	Depends only on past and present inputs (no future inputs).
Real-Time Use	Suitable for real-time systems (e.g., EEG, audio processing).
Group Delay	Introduces a constant delay (group delay) due to the convolution process.

Non-Causal Filters

A non-causal filter is a filter whose output depends on future inputs (in addition to past and present inputs).

Mathematically, for a non-causal filter:

\[ y(t) = \int_{-\infty}^{\infty} h(t - \tau) x(\tau) d\tau \]

or in discrete time:

\[ y[n] = \sum_{k=-M}^{N} h[k] \cdot x[n-k] \]

where: $h[k]$ can be non-zero for negative $k$ (future inputs).

Key Characteristics

Feature	Description
Impulse Response	$h(t) \neq 0$ for $t < 0$ (includes future inputs).
Output Dependency	Depends on past, present, and future inputs.
Zero Delay	Can be designed to have zero group delay (no phase distortion).
Real-Time Use	Not suitable for real-time systems (requires future data).

Advantages

Zero Delay: Non-causal filters can process signals without introducing delay or phase distortion.
Phase Integrity: Preserves the shape of the signal (no group delay).

Disadvantages

Not Real-Time: Requires future data, making it unsuitable for real-time applications.
Offline Processing: Typically used for post-processing of recorded signals.

Example: Non-Causal Moving Average Filter

A moving average filter is a classic example of a non-causal filter.
For a symmetric window (e.g., $[-2, -1, 0, 1, 2]$), the output at time $n$ depends on the current and future samples.
The group delay is zero because the filter is symmetric.

Comparison: Causal vs. Non-Causal Filters

Feature	Causal Filter	Non-Causal Filter
Definition	Output depends only on past/present inputs.	Output depends on past/present/future inputs.
Impulse Response	$h(t) = 0$ for $t < 0$.	$h(t) \neq 0$ for $t < 0$.
Group Delay	Introduces a constant delay.	Can have zero delay.
Real-Time Use	Suitable for real-time systems.	Not suitable for real-time systems.
Phase Distortion	May introduce phase shifts.	Preserves phase (no group delay).
Examples	Low-pass FIR/IIR filters, Butterworth filters.	Moving average, zero-phase filters.
Applications	EEG preprocessing, audio processing, real-time signal enhancement.	Offline signal analysis, post-processing.

Implications for Signal Processing

Real-Time Systems

Causal filters are essential for real-time applications, such as:
- EEG/MEG preprocessing (removing artifacts in real-time).
- Audio processing (real-time equalization).
- Biomedical signal monitoring (e.g., ECG, respiration).
Non-causal filters cannot be used in real-time because they require future data.

Offline Processing

Non-causal filters are ideal for offline analysis, where the entire signal is available for processing. Examples include:
- Post-processing EEG data (e.g., zero-phase filtering).
- Image processing (e.g., blurring, edge detection).
- Scientific data analysis (e.g., smoothing, trend removal).

Zero-Phase Filtering

Zero-phase filters are a hybrid approach:
1. Apply a causal filter forward.
2. Apply the same filter backward to eliminate phase distortion.
3. The result is a non-causal zero-phase filter with no group delay.
Example: In EEG analysis, zero-phase filters are used to remove line noise (50/60 Hz) without introducing delays.

Group Delay

Definition: The time delay introduced by the filter, which is the difference between the input and output phase at a specific frequency.

Cause: In a causal filter, the convolution of the input signal with the impulse response only includes part of the signal initially (until the full impulse response is applied).

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

Calculation:

For a linear-phase FIR filter, the group delay is:

\[ \tau_g = \frac{N-1}{2} \]

where $N$ is the filter order (length of the filter window, in samples).

Impact: The output signal is delayed by the group delay, which can be problematic for applications requiring zero delay.

Example: Causal FIR Filter

A low-pass FIR filter with symmetric coefficients has a constant group delay.
The output is a delayed version of the input signal.

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).

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 a convenient measure of the linearity of the phase with respect to frequency in a modulation system.

Introduction to the Z-Transform

The Z-transform is a powerful mathematical tool used in digital signal processing (DSP) to analyze and design discrete-time systems, such as digital filters. It generalizes the Discrete-Time Fourier Transform (DTFT) by incorporating exponential growth or decay into the signal analysis. This makes it particularly useful for studying the stability, frequency response, and design of Infinite Impulse Response (IIR) filters, analogies to the Laplace transform for continuous-time systems.

The Z-transform converts a discrete-time signal $x[n]$ (where $n$ is an integer representing the sample index) into a complex frequency domain representation $X(z)$. It is defined as:

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

where:
$x[n]$ is the discrete-time input signal,
$z$ is the complex variable in the Z-domain, often written as $z = re^{i\omega}$ (polar form), where:
$r$ is the magnitude (growth/decay factor),
$\omega$ is the phase angle (frequency in radians/sample).

Purpose

The Z-transform extends the Fourier transform to include exponential growth/decay in the signal, making it ideal for analyzing IIR filters and causal systems.
It provides a unified framework for studying the stability, frequency response, and design of digital filters.
In the Z-domain, systems are represented as rational functions (ratios of polynomials), which can be analyzed using poles and zeros.

The Discrete-Time Fourier Transform (DTFT) is a special case of the Z-transform where the magnitude of $z$ is 1 (i.e., $z = e^{i\omega}$). It represents the frequency response of a signal:

\[ X(e^{i\omega}) = \sum_{n=-\infty}^{\infty} x[n] \cdot e^{-i\omega n} \]

Interpretation: The DTFT treats the signal as a sum of complex exponentials (sines and cosines) of different frequencies.

Limitation: The DTFT assumes the signal is non-growing (i.e., $r = 1$). It cannot represent systems with exponential growth or decay.

Z-Transform The Z-transform generalizes the DTFT by allowing $z$ to vary in magnitude ($r$) and phase ($\omega$):

\[ X(z) = \sum_{n=-\infty}^{\infty} x[n] \cdot z^{-n} = \sum_{n=-\infty}^{\infty} x[n] \cdot (re^{i\omega})^{-n} \]

Advantage: By varying $r$, the Z-transform can represent systems with exponential growth or decay (e.g., IIR filters, unstable systems).

Linear Time-Invariant (LTI) Systems

The Z-transform is particularly useful for analyzing Linear Time-Invariant (LTI) systems, which are the foundation of digital filters.

For an LTI system with impulse response $h[n]$ and input $x[n]$, the output $y[n]$ is the convolution of $x[n]$ and $h[n]$:

\[ y[n] = x[n] * h[n] = \sum_{k=-\infty}^{\infty} x[k] \cdot h[n-k] \]

In the Z-domain, convolution becomes multiplication:

\[ Y(z) = X(z) \cdot H(z) \]

where:
$X(z)$ is the Z-transform of the input signal,
$H(z)$ is the Z-transform of the system’s impulse response (called the transfer function),
$Y(z)$ is the Z-transform of the output signal.

Transfer Function

The transfer function $H(z)$ describes how the system transforms input signals into output signals in the Z-domain.

It is typically a rational function (a ratio of polynomials):

\[ H(z) = \frac{B(z)}{A(z)} = \frac{b_0 + b_1 z^{-1} + \cdots + b_N z^{-N}}{1 + a_1 z^{-1} + \cdots + a_M z^{-M}}\]

where:
$B(z)$ is the numerator polynomial (feedforward coefficients),
$A(z)$ is the denominator polynomial (feedback coefficients).

Poles: Roots of $A(z)$ (feedback).
Zeros: Roots of $B(z)$ (feedforward).

Designing Filters

To design a filter, we define a desired frequency response (e.g., low-pass, band-pass) and find $H(z)$ that approximates this response.

The inverse Z-transform can then be used to obtain the impulse response $h[n]$ from $H(z)$.

Poles

Definition: Poles are the roots of the denominator polynomial A(z)A(z)A(z) in H(z)H(z)H(z).
Interpretation: A pole at z=pz = pz=p means the system has a natural mode that grows or decays exponentially with time.

Stability:

If all poles lie inside the unit circle ($\vert p \vert < 1$), the system is stable (output decays to zero for bounded input).
If any pole lies outside the unit circle ($\vert p \vert > 1$), the system is unstable (output grows without bound).
Poles on the unit circle ($\vert p \vert = 1$) indicate marginal stability (e.g., sustained oscillations).

Zeros

Definition: Zeros are the roots of the numerator polynomial $B(z)$ in $H(z)$. Interpretation: A zero at $z = q$ means the system attenuates specific frequency components.

Poles and zeros are plotted on the complex plane (Z-plane).

The unit circle ($\vert z \vert = 1$) represents the frequency axis for stable systems.

Creating IIR filter:

fs = 100
f = filter_create(fprototype = :butterworth,
                  ftype = :bp,
                  order = 8,
                  cutoff = (10, 20),
                  fs = fs)

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

Gain:

f.k

2.3959644103776224e-5

Drawing pole-zero map:

plot_polezero(f.p, f.z)

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.

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.

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,
                    bw = 1,
                    fs = 256)
t = linspace(-1, 1, length(flt))
NeuroAnalyzer.plot(t, flt)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/kJl0u/src/scenes.jl:264

And its frequency response:

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

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/kJl0u/src/scenes.jl:264

Detailed frequency response:

plot_filter(fprototype = :fir,
            ftype = :lp,
            cutoff = 10,
            bw = 1,
            order = 91,
            fs = 256,
            gui = false)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/kJl0u/src/scenes.jl:264

Same FIR band pass filter, but using Bartlett window:

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

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/kJl0u/src/scenes.jl:264

And its frequency response:

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

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/kJl0u/src/scenes.jl:264

Other impulse responses:

Low pass FIR filter:

flt = filter_create(fprototype = :fir, ftype = :lp, cutoff = 10, order = 91, bw = 0.5, fs = 256)
t = linspace(-1, 1, length(flt))
NeuroAnalyzer.plot(t, flt)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/kJl0u/src/scenes.jl:264

High pass FIR filter:

flt = filter_create(fprototype = :fir, ftype = :hp, cutoff = 10, order = 91, bw = 0.5, fs = 256)
NeuroAnalyzer.plot(t, flt)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/kJl0u/src/scenes.jl:264

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:

accurate for FIR filters only
for very large N use block convolution
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:

filter preparation (filter_create())
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)

Tip: bw is transition band width in Hz and a is attenuation in dB.

Tip: This estimation is just the entry point for further investigation. Filter response should be analyzed using plot_filter().

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)

NeuroAnalyzer.NEURO(NeuroAnalyzer.HEADER(Dict{Symbol, Any}(:weight => -1, :id => "", :middle_name => "", :height => -1, :head_circumference => -1, :handedness => "", :last_name => "528004  SIT 52, 20220831-122227-{d589f756-53fc-4f1b-915d-6e3b8c1560ad}", :first_name => ""), Dict{Symbol, Any}(:epoch_id => "", :channel_type => ["eeg", "eeg", "eeg", "eeg", "eeg", "eeg", "eeg", "eeg", "eeg", "eeg"  …  "eeg", "eeg", "ref", "ref", "eeg", "eeg", "eeg", "eog", "eog", "ecg"], :label => ["Fp1", "Fp2", "F3", "F4", "C3", "C4", "P3", "P4", "O1", "O2"  …  "T5", "T6", "A1", "A2", "Fz", "Cz", "Pz", "EOG1", "EOG2", "ECG"], :prefiltering => ["HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz"  …  "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz"], :gain => [0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218  …  0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218], :data_type => "eeg", :recording_notes => "", :recording_date => "31.08.22", :sampling_rate => 256, :file_type => "EDF"…), Dict(:name => "", :design => "", :notes => "")), ["filter_apply(OBJ, ch=[1], dir=twopass)", "filter_apply(OBJ, ch=[1], dir=twopass)", "filter_apply(OBJ, ch=[1], dir=twopass)", "filter_apply(OBJ, ch=[2], dir=twopass)", "filter_apply(OBJ, ch=[2], dir=twopass)", "filter_apply(OBJ, ch=[2], dir=twopass)", "filter_apply(OBJ, ch=[3], dir=twopass)", "filter_apply(OBJ, ch=[3], dir=twopass)", "filter_apply(OBJ, ch=[3], dir=twopass)", "filter_apply(OBJ, ch=[4], dir=twopass)"  …  "add_marker(OBJ, id=NA, start=62.0, len=1.0, value=DELETED, ch=0)", "trim(OBJ, seg=(62.0, 67.0), keep=false", "add_marker(OBJ, id=NA, start=74.0, len=1.0, value=DELETED, ch=0)", "trim(OBJ, seg=(74.0, 76.0), keep=false", "add_marker(OBJ, id=NA, start=119.0, len=1.0, value=DELETED, ch=0)", "trim(OBJ, seg=(119.0, 122.0), keep=false", "add_marker(OBJ, id=NA, start=348.0, len=1.0, value=DELETED, ch=0)", "trim(OBJ, seg=(348.0, 352.0), keep=false", "trim(OBJ, seg=(0, 1000), keep=true", "filter_apply(obj; ch=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24], dir=twopass)"], 0×5 DataFrame
 Row │ id      start    length   value   channel 
     │ String  Float64  Float64  String  Int64   
─────┴───────────────────────────────────────────, 23×9 DataFrame
 Row │ label   loc_radius  loc_theta  loc_x    loc_y    loc_z    loc_radius_sp ⋯
     │ String  Float64     Float64    Float64  Float64  Float64  Float64       ⋯
─────┼──────────────────────────────────────────────────────────────────────────
   1 │ Fp1           1.0       108.0    -0.31     0.95    -0.03            1.0 ⋯
   2 │ Fp2           1.0        72.0     0.31     0.95    -0.03            1.0
   3 │ F3            0.65      129.0    -0.55     0.67     0.5             1.0
   4 │ F4            0.65       51.0     0.55     0.67     0.5             1.0
   5 │ C3            0.51      180.0    -0.72     0.0      0.7             1.0 ⋯
   6 │ C4            0.51        0.0     0.72     0.0      0.7             1.0
   7 │ P3            0.65      231.0    -0.55    -0.67     0.5             1.0
   8 │ P4            0.65      309.0     0.55    -0.67     0.5             1.0
   9 │ O1            1.0       252.0    -0.31    -0.95    -0.03            1.0 ⋯
  10 │ O2            1.0       288.0     0.31    -0.95    -0.03            1.0
  11 │ F7            1.0       144.0    -0.81     0.59    -0.03            1.0
  ⋮  │   ⋮         ⋮           ⋮         ⋮        ⋮        ⋮           ⋮       ⋱
  14 │ T4            1.0         0.0     1.0      0.0     -0.03            1.0
  15 │ T5            1.0       216.0    -0.81    -0.59    -0.03            1.0 ⋯
  16 │ T6            1.0       324.0     0.81    -0.59    -0.03            1.0
  17 │ A1            1.0       192.0    -0.92    -0.23    -0.55            1.1
  18 │ A2            1.0       -12.0     0.92    -0.23    -0.55            1.1
  19 │ Fz            0.51       90.0     0.0      0.72     0.7             1.0 ⋯
  20 │ Cz            0.0         0.0     0.0      0.0      1.0             1.0
  21 │ Pz            0.51      270.0     0.0     -0.72     0.7             1.0
  22 │ EOG1          1.01      149.0    -0.87     0.51    -0.37            1.0
  23 │ EOG2          1.01       31.0     0.87     0.51    -0.37            1.0 ⋯
                                                    3 columns and 2 rows omitted, [0.0, 0.0039, 0.0078, 0.0117, 0.0156, 0.0195, 0.0234, 0.0273, 0.0312, 0.0352  …  999.9648, 999.9688, 999.9727, 999.9766, 999.9805, 999.9844, 999.9883, 999.9922, 999.9961, 1000.0], [0.0, 0.0039, 0.0078, 0.0117, 0.0156, 0.0195, 0.0234, 0.0273, 0.0312, 0.0352  …  999.9648, 999.9688, 999.9727, 999.9766, 999.9805, 999.9844, 999.9883, 999.9922, 999.9961, 1000.0], [-1.3373115403947367e-7 0.965874181893055 … 3.3420881950453407 -9.78082278502157e-8; 1.668413851518835e-8 -1.0656693698351298 … 2.209998047433227 -1.0494827423768172e-8; … ; -3.5604833747981957e-8 -1.1228098809810545 … -3.4418165483537275 -4.160514648354763e-7; -2.0092275221372802e-7 -2.767021563420194 … 1.3115379696702059 -1.3845311555016337e-6;;;])

Tip: 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)

Tip: bw is transition band width in Hz and a is attenuation in dB.

Tip: This estimation is just the entry point for further investigation. Filter response should be analyzed using plot_filter().

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)

Tip: 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))
NeuroAnalyzer.plot(t, flt)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/kJl0u/src/scenes.jl:264

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)

Tip: bw is transition band width in Hz and a is attenuation in dB.

Tip: This estimation is just the entry point for further investigation. Filter response should be analyzed using plot_filter().

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))
NeuroAnalyzer.plot(t, flt)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/kJl0u/src/scenes.jl:264

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)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/kJl0u/src/scenes.jl:264

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)

Tip: This estimation is just the entry point for further investigation. Filter response should be analyzed using plot_filter().

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);

Tip: 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)

NeuroAnalyzer.plot(t, y1)
NeuroAnalyzer.plot(t, y2)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/kJl0u/src/scenes.jl:264

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/kJl0u/src/scenes.jl:264

$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

NeuroAnalyzer.plot(t, y1)
NeuroAnalyzer.plot(t, y2)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/kJl0u/src/scenes.jl:264

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/kJl0u/src/scenes.jl:264

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]);

Tip: 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

NeuroAnalyzer.plot(t, y1)
NeuroAnalyzer.plot(t, y2)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/kJl0u/src/scenes.jl:264

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/kJl0u/src/scenes.jl:264

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
NeuroAnalyzer.plot(t, s)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/kJl0u/src/scenes.jl:264

Band pass Gaussian filter at 2 Hz:

s_new = filter_g(s; fs = fs, f = 2)
NeuroAnalyzer.plot(t, s_new)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/kJl0u/src/scenes.jl:264

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,
                     bw = 0.5,
                     preview = true)

Tip: When preview=true, signal will not be filtered.

Another way is to use plot_filter():

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

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/kJl0u/src/scenes.jl:264

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)

Feature	Description
Impulse Response	\(h(t) = 0\) for \(t < 0\) (no values before \(t = 0\)).
Output Dependency	Depends only on past and present inputs (no future inputs).
Real-Time Use	Suitable for real-time systems (e.g., EEG, audio processing).
Group Delay	Introduces a constant delay (group delay) due to the convolution process.

Feature	Description
Impulse Response	\(h(t) \neq 0\) for \(t < 0\) (includes future inputs).
Output Dependency	Depends on past, present, and future inputs.
Zero Delay	Can be designed to have zero group delay (no phase distortion).
Real-Time Use	Not suitable for real-time systems (requires future data).