Preprocessing with Finite Impulse Response (FIR) Filters

Author

Dr. Ole Bialas

Download Materials

import numpy as np
import matplotlib.pyplot as plt
import mne

EEG and MEG recordings are broadband signals that mix neural activity across a wide range of frequencies with noise and artifacts of both biological and environmental origin. Filtering is the primary tool for separating these components: by selectively passing or suppressing specific frequency ranges, we can remove slow baseline drifts, attenuate powerline interference, or isolate a particular neural oscillation of interest such as alpha or theta activity.

However, filtering is not a neutral operation. Every filter makes a trade-off between frequency selectivity and time-domain fidelity. A filter that is very sharp in frequency must be long in time, and a long filter can smear, delay, or distort the very signals it is meant to clean. Understanding these trade-offs is essential for correctly interpreting filtered neural data — in particular, it is easy to mistake a filtering artifact for a genuine neural effect.

In this session, we’ll explore common applications of filters in preprocessing and visualize their effects in the time and frequency domain. The first two sections use magnetometer recordings from the MNE sample dataset. The later sections use synthetic signals to make the artifacts clearly visible under controlled conditions.

Setup

Utility Functions

def plot_filtered(before, after, tmin, tmax):
    plt.plot(before.times, before.get_data()[0], label="original")
    plt.plot(after.times, after.get_data()[0], label="after filtering")
    plt.xlim(tmin, tmax)
    plt.xlabel("Time [s]")
    plt.ylabel("Amplitude [a.u.]")
    plt.legend()
    plt.show()

def make_oscillation(dur, freq, t_start, t_stop, sfreq=500, plot=False):
    n_samples = int(sfreq * dur)
    t = np.linspace(0, dur, n_samples)
    signal = np.zeros(n_samples)
    burst_mask = (t >= t_start) & (t <= t_stop)
    signal[burst_mask] = np.sin(2 * np.pi * freq * t[burst_mask])
    info = mne.create_info( ch_names=['ch1'], sfreq=sfreq, ch_types=['eeg'])
    raw = mne.io.RawArray(signal[np.newaxis, :], info)
    if plot:
        plt.plot(raw.times, raw.get_data()[0])
        plt.xlabel("Time [s]")
        plt.ylabel(r"Amplitude $\mu$V")
        plt.show()
    return raw

def make_impulse(dur, t_impulse, width, amplitude=1.0, sfreq=500, plot=False):
    n_samples = int(sfreq * dur)
    signal = np.zeros(n_samples)
    impulse_start = int(t_impulse * sfreq)
    impulse_stop = int((t_impulse+width) * sfreq)
    signal[impulse_start:impulse_stop] = amplitude
    info = mne.create_info(ch_names=['ch1'], sfreq=sfreq, ch_types=['eeg'])
    raw = mne.io.RawArray(signal[np.newaxis, :], info)
    if plot:
        plt.plot(raw.times, raw.get_data()[0])
        plt.xlabel("Time [s]")
        plt.ylabel(r"Amplitude $\mu$V")
        plt.title("Impulse Signal")
        plt.show()
    return raw

def make_step(dur, t_step, amplitude=1.0, sfreq=500, plot=False):
    n_samples = int(sfreq * dur)
    t = np.linspace(0, dur, n_samples)
    signal = np.zeros(n_samples)
    signal[t >= t_step] = amplitude
    info = mne.create_info(ch_names=['ch1'], sfreq=sfreq, ch_types=['eeg'])
    raw = mne.io.RawArray(signal[np.newaxis, :], info)
    if plot:
        plt.plot(raw.times, raw.get_data()[0])
        plt.xlabel("Time [s]")
        plt.ylabel(r"Amplitude $\mu$V")
        plt.title("Step Function")
        plt.show()
    return raw

class utils:
    plot_filtered = plot_filtered
    make_oscillation = make_oscillation
    make_impulse = make_impulse
    make_step = make_step

Data Preparation

sample_data_folder = mne.datasets.sample.data_path()
sample_data_raw_file = sample_data_folder / "MEG" / "sample" / "sample_audvis_raw.fif"
raw = mne.io.read_raw_fif(sample_data_raw_file)
raw.crop(0, 60).pick(picks=["mag"]).load_data();
raw.del_proj("all");

Opening raw data file /home/atleeri/mne_data/MNE-sample-data/MEG/sample/sample_audvis_raw.fif...
    Read a total of 3 projection items:
        PCA-v1 (1 x 102)  idle
        PCA-v2 (1 x 102)  idle
        PCA-v3 (1 x 102)  idle
    Range : 25800 ... 192599 =     42.956 ...   320.670 secs
Ready.
Reading 0 ... 36037  =      0.000 ...    60.000 secs...

Section 1: Removing Drift and Powerline Artifacts

Background

EEG and MEG recordings are contaminated by different artifacts. There are two ubiquitous noise sources that can be addressed with filters:

Drift is a slow, low-frequency fluctuation in the baseline signal. It can arise from electrode impedance changes, slow movements of the subject, or gradual shifts in amplifier offset. Because neural signals of interest rarely occur below ~0.1–1 Hz, a highpass filter can suppress drift without removing meaningful neural activity.
Powerline noise is a narrowband artifact at the frequency of the AC power supply (50 Hz in Europe, 60 Hz in North America) and its integer multiples (harmonics). It is picked up by the recording equipment from the electrical environment. A lowpass filter or notch filter can remove the powerline noise from the recording.

Exercises

In the following exercises, you are going to filter raw MEG using lowpass, highpass, bandpass (combination of low- and highpass) and notch filters. To visualize the effect of filtering, you are going to plot the filtered signal and its power spectral density (PSD). Here are some useful code examples:

Code	Description
`raw.plot(duration=10)`	Plot a `10` second segment of the data stored in `raw`
`raw2 = raw.copy()`	Create a new copy of the `raw` object named `raw2`
`raw.filter(l_freq=1, h_freq=None)`	Apply a 1 Hz highpass filter to `raw`
`raw.filter(l_freq=None, h_freq=40)`	Apply a 40 Hz lowpass filter to `raw`
`raw.filter(l_freq=1, h_freq=40)`	Apply a 1 Hz to 40 Hz bandpass filter to `raw`
`raw.notch_filter(freqs=(50, 100))`	Apply a notch filter at `50` and `100` Hz to `raw`
`raw.notch_filter(freqs=(50, 100), notch_widths=2)`	Set the width of the filter’s notches to `2` Hz
`spectrum = raw.compute_psd()`	Compute the power spectral density (PSD) of `raw`
`spectrum.plot(average=True)`	Plot the PSD, averaged across all sensors
`spectrum.plot(average=True, xscale="log")`	Plot the average PSD with a logarithmically scaled x-axis

Run the cell below to plot one minute of raw MEG data. Several channels show drifts — slow fluctuations in the measured signals that can be caused for example by changes in temperature or environmental magnetic fields. These drift artifacts can be removed using a highpass filter.

raw.plot(duration=60);

Example: Create a new copy of raw called raw_hp. Then, apply a 0.1 Hz highpass filter to raw_hp and plot the filtered data. As you can see, the filter did not remove the drifts completely.

raw_hp = raw.copy()
raw_hp.filter(l_freq=0.1, h_freq=None)
raw_hp.plot(duration=60);

Filtering raw data in 1 contiguous segment
Setting up high-pass filter at 0.1 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 0.10
- Lower transition bandwidth: 0.10 Hz (-6 dB cutoff frequency: 0.05 Hz)
- Filter length: 19821 samples (33.001 s)

Exercise: Create a new copy of raw called raw_hp. Then, apply a 1 Hz highpass filter to raw_hp and plot the filtered data.

Solution

raw_hp = raw.copy()
raw_hp.filter(l_freq=1, h_freq=None)
raw_hp.plot(duration=60);

Filtering raw data in 1 contiguous segment
Setting up high-pass filter at 1 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 1.00
- Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)
- Filter length: 1983 samples (3.302 s)

Exercise: Create a new copy of raw called raw_hp. Then, apply a 0.3 Hz highpass filter to raw_hp and plot the filtered data.

Solution

raw_hp = raw.copy()
raw_hp.filter(l_freq=0.3, h_freq=None)
raw_hp.plot(duration=60);

Filtering raw data in 1 contiguous segment
Setting up high-pass filter at 0.3 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 0.30
- Lower transition bandwidth: 0.30 Hz (-6 dB cutoff frequency: 0.15 Hz)
- Filter length: 6607 samples (11.000 s)

Example: Compute the power spectral density (PSD) which shows the frequency content of the raw signal in decibels per Hz. Observe the sharp spikes at 60 Hz and its multiples. This is an artifact produced by the alternating current of the power grid.

spectrum = raw.compute_psd()
spectrum.plot()

Effective window size : 3.410 (s)
Plotting power spectral density (dB=True).

Exercise: Create a new copy of raw called raw_lp. Then, apply a 40 Hz lowpass filter to raw_lp and compute and plot the filtered signal’s PSD.

Solution

raw_lp = raw.copy()
raw_lp.filter(l_freq=None, h_freq=40)
spectrum = raw_lp.compute_psd()
spectrum.plot()

Filtering raw data in 1 contiguous segment
Setting up low-pass filter at 40 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal lowpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Upper passband edge: 40.00 Hz
- Upper transition bandwidth: 10.00 Hz (-6 dB cutoff frequency: 45.00 Hz)
- Filter length: 199 samples (0.331 s)

Effective window size : 3.410 (s)
Plotting power spectral density (dB=True).

Exercise: Create a new copy of raw called raw_lp. Then, apply a 10 Hz lowpass filter to raw_lp and compute and plot the filtered signal’s PSD. Also plot the filtered signal using the .plot() method.

Tip: Use xscale="log" to get a better view of the lower end of the PSD.

Solution

raw_lp = raw.copy()
raw_lp.filter(l_freq=None, h_freq=10)
spectrum = raw_lp.compute_psd()
spectrum.plot(xscale="log")

Filtering raw data in 1 contiguous segment
Setting up low-pass filter at 10 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal lowpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Upper passband edge: 10.00 Hz
- Upper transition bandwidth: 2.50 Hz (-6 dB cutoff frequency: 11.25 Hz)
- Filter length: 793 samples (1.320 s)

Effective window size : 3.410 (s)
Plotting power spectral density (dB=True).

Solution

raw_lp.plot();

Example: Create a new copy of raw called raw_notch. Then, apply a notch filter at 180 and 240 Hz to raw_notch and compute and plot the filtered signal’s PSD.

raw_notch = raw.copy()
raw_notch.notch_filter(freqs=(180, 240))
spectrum = raw_notch.compute_psd()
spectrum.plot(average=True)

Filtering raw data in 1 contiguous segment
Setting up band-stop filter

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandstop filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower transition bandwidth: 0.50 Hz
- Upper transition bandwidth: 0.50 Hz
- Filter length: 3965 samples (6.602 s)

Effective window size : 3.410 (s)
Plotting power spectral density (dB=True).

Exercise: Create a new copy of raw called raw_notch. Then, apply a notch filter at 60, 120, 180 and 240 Hz to raw_notch and compute and plot the filtered signal’s PSD.

Solution

raw_notch = raw.copy()
raw_notch.notch_filter(freqs=(60, 120, 180, 240))
spectrum = raw_notch.compute_psd()
spectrum.plot(average=True)

Filtering raw data in 1 contiguous segment
Setting up band-stop filter

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandstop filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower transition bandwidth: 0.50 Hz
- Upper transition bandwidth: 0.50 Hz
- Filter length: 3965 samples (6.602 s)

Effective window size : 3.410 (s)
Plotting power spectral density (dB=True).

Exercise: Repeat the notch filter from the previous exercises but set notch_widths=2.

Solution

raw_notch = raw.copy()
raw_notch.notch_filter(freqs=(60, 120, 180, 240), notch_widths=2)
spectrum = raw_notch.compute_psd()
spectrum.plot(average=True)

Filtering raw data in 1 contiguous segment
Setting up band-stop filter

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandstop filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower transition bandwidth: 0.50 Hz
- Upper transition bandwidth: 0.50 Hz
- Filter length: 3965 samples (6.602 s)

Effective window size : 3.410 (s)
Plotting power spectral density (dB=True).

Example: Create a new copy of raw called raw_bp. Then apply a bandpass filter with a low cutoff of 0.1 Hz and a high cutoff of 50 Hz to raw_bp and compute and plot the filtered signal’s PSD.

raw_bp = raw.copy()
raw_bp.filter(l_freq=0.1, h_freq=50)
spectrum = raw_bp.compute_psd()
spectrum.plot()

Filtering raw data in 1 contiguous segment
Setting up band-pass filter from 0.1 - 50 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 0.10
- Lower transition bandwidth: 0.10 Hz (-6 dB cutoff frequency: 0.05 Hz)
- Upper passband edge: 50.00 Hz
- Upper transition bandwidth: 12.50 Hz (-6 dB cutoff frequency: 56.25 Hz)
- Filter length: 19821 samples (33.001 s)

Effective window size : 3.410 (s)
Plotting power spectral density (dB=True).

Exercise: Create a new copy of raw called raw_bp. Then apply a bandpass filter with a low cutoff of 0.3 Hz and a high cutoff of 40 Hz to raw_bp and compute and plot the filtered signal’s PSD.

Solution

raw_bp = raw.copy()
raw_bp.filter(0.3, 40)
spectrum = raw_bp.compute_psd()
spectrum.plot()

Filtering raw data in 1 contiguous segment
Setting up band-pass filter from 0.3 - 40 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 0.30
- Lower transition bandwidth: 0.30 Hz (-6 dB cutoff frequency: 0.15 Hz)
- Upper passband edge: 40.00 Hz
- Upper transition bandwidth: 10.00 Hz (-6 dB cutoff frequency: 45.00 Hz)
- Filter length: 6607 samples (11.000 s)

Effective window size : 3.410 (s)
Plotting power spectral density (dB=True).

Exercise: Create a new copy of raw called raw_bp. Then apply a bandpass filter that isolates activity in the theta band (4-8 Hz) and compute and plot the filtered signal’s PSD (TIP: use a logarithmically scaled x-axes). Also plot the filtered signal using the .plot() method.

Solution

raw_bp = raw.copy()
raw_bp.filter(4, 8)
spectrum = raw_bp.compute_psd()
spectrum.plot(xscale="log")

Filtering raw data in 1 contiguous segment
Setting up band-pass filter from 4 - 8 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 4.00
- Lower transition bandwidth: 2.00 Hz (-6 dB cutoff frequency: 3.00 Hz)
- Upper passband edge: 8.00 Hz
- Upper transition bandwidth: 2.00 Hz (-6 dB cutoff frequency: 9.00 Hz)
- Filter length: 993 samples (1.653 s)

Effective window size : 3.410 (s)
Plotting power spectral density (dB=True).

Solution

raw_bp.plot();

Exercise: Create a new copy of raw called raw_bp. Then apply a bandpass filter that isolates activity in the alpha band (8-13 Hz) and compute and plot the filtered signal’s PSD (TIP: use a logarithmically scaled x-axes). Also plot the filtered signal using the .plot() method.

Solution

raw_bp = raw.copy()
raw_bp.filter(8, 13)
spectrum = raw_bp.compute_psd()
spectrum.plot(xscale="log")

Filtering raw data in 1 contiguous segment
Setting up band-pass filter from 8 - 13 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 8.00
- Lower transition bandwidth: 2.00 Hz (-6 dB cutoff frequency: 7.00 Hz)
- Upper passband edge: 13.00 Hz
- Upper transition bandwidth: 3.25 Hz (-6 dB cutoff frequency: 14.62 Hz)
- Filter length: 993 samples (1.653 s)

Effective window size : 3.410 (s)
Plotting power spectral density (dB=True).

Solution

raw_bp.plot();

Section 2: Filter properties

Background

The magnitude response of a FIR (Finite Impulse Response) filter consists of a passband (frequencies that pass through unchanged) and a stopband (frequencies that are suppressed). The frequency range over which the filter transitions from passband to stopband is called the transition bandwidth. A narrower transition band produces a sharper, more frequency-selective filter — but this comes at a cost, because FIR filter length and transition bandwidth are inversely related:

\text{filter length} \propto \frac{1}{\text{transition bandwidth}}

A narrower transition band therefore requires a longer filter, which means:

Better frequency selectivity — steeper roll-off between passband and stopband
More time-domain spread — the longer impulse response can smear or distort brief neural events in time
Greater data requirements — if the required filter length exceeds the length of the recording, MNE will warn you and truncate the filter, producing a degraded frequency response

Additionally, because the filter is finite, the magnitude response never achieves perfect, flat attenuation in the stopband. Instead it oscillates slightly around the target level — an effect known as stopband ripple.

Exercises

In the following exercises you are going to create filters with different transition bandwidths and visualize their impulse and magnitude response to see the relationship between filter length and steepness.

Code	Description
`filt = mne.filter.create_filter(data, sfreq, l_freq, h_freq)`	Create a filter for the given `data` and sampling frequency (`sfreq`) with the cutoffs `l_freq` and `h_freq` and store the filter coefficients in `filt`
`filt = mne.filter.create_filter(data, sfreq, l_freq, h_freq, h_trans_bandwidth=5)`	Create a filter and set the width of the upper transition band to `5` Hz
`filt = mne.filter.create_filter(data, sfreq, l_freq, h_freq, l_trans_bandwidth=0.5)`	Create a filter and set the width of the lower transition band to `0.5` Hz
`mne.viz.plot_filter(filt, sfreq, plot=("time", "magnitude"))`	Plot the impulse response (`"time"`) and magnitude response (`"magnitude"`) of `filt`
`mne.viz.plot_filter(filt, sfreq, plot=("time", "magnitude"), flim=(0.01, 5))`	Same as above, but limit the frequency axis of the magnitude plot to `0.01`–`5` Hz

Example: Make a 40 Hz lowpass filter for raw. Then, plot its "time" and "magnitude" response. The top panel shows the filter’s impulse response — its coefficients over time, which are convolved with the signal during filtering. The bottom panel shows the filter’s magnitude response, which describes how much each frequency is attenuated (this filter passes everything below 40 Hz and suppresses everything above it).

filt = mne.filter.create_filter(
    raw.get_data(), raw.info["sfreq"], l_freq=None, h_freq=40
    )
mne.viz.plot_filter( filt, raw.info["sfreq"], plot=("time", "magnitude"));

Setting up low-pass filter at 40 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal lowpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Upper passband edge: 40.00 Hz
- Upper transition bandwidth: 10.00 Hz (-6 dB cutoff frequency: 45.00 Hz)
- Filter length: 199 samples (0.331 s)

Exercise: Make a 40 Hz lowpass filter for raw and set h_trans_bandwidth=1 (the default is 10 Hz). How does reducing the transition band width affect the "time" and "magnitude" response of the filter?

Solution

filt = mne.filter.create_filter(
    raw.get_data(), raw.info["sfreq"], l_freq=None, h_freq=40, h_trans_bandwidth=1,
    )
mne.viz.plot_filter( filt, raw.info["sfreq"], plot=("time", "magnitude"));

Setting up low-pass filter at 40 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal lowpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Upper passband edge: 40.00 Hz
- Upper transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 40.50 Hz)
- Filter length: 1983 samples (3.302 s)

Exercise: Make a 40 Hz lowpass filter for raw and set h_trans_bandwidth=20. How does increasing the transition band width affect the "time" and "magnitude" response of the filter?

filt = mne.filter.create_filter(
    raw.get_data(), raw.info["sfreq"], l_freq=None, h_freq=40, h_trans_bandwidth=20,
    )
mne.viz.plot_filter( filt, raw.info["sfreq"], plot=("time", "magnitude"));

Setting up low-pass filter at 40 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal lowpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Upper passband edge: 40.00 Hz
- Upper transition bandwidth: 20.00 Hz (-6 dB cutoff frequency: 50.00 Hz)
- Filter length: 101 samples (0.168 s)

Exercise: Make a 1 Hz highpass filter for raw and plot its "time" and "magnitude" response. Use the flim argument to restrict the plot to frequencies between 0.01 and 5 Hz.

Solution

filt = mne.filter.create_filter(
    raw.get_data(), raw.info["sfreq"], l_freq=1, h_freq=None)
mne.viz.plot_filter( filt, raw.info["sfreq"], flim=(0.01, 5), plot=("time", "magnitude"));

Setting up high-pass filter at 1 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 1.00
- Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)
- Filter length: 1983 samples (3.302 s)

Exercise: Make a 1 Hz highpass filter for raw and set l_trans_bandwidth=0.5 (the default is 1 Hz). Then, plot its "time" and "magnitude" response between 0.01 and 5 Hz. How does reducing the transition band width affect the filter’s time and magnitude response?

Solution

filt = mne.filter.create_filter(
    raw.get_data(), raw.info["sfreq"], l_freq=1, h_freq=None, l_trans_bandwidth=0.5)
mne.viz.plot_filter( filt, raw.info["sfreq"], flim=(0.001, 5), plot=("time", "magnitude"));

Setting up high-pass filter at 1 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 1.00
- Lower transition bandwidth: 0.50 Hz (-6 dB cutoff frequency: 0.75 Hz)
- Filter length: 3965 samples (6.602 s)

Exercise: Make a 1 Hz highpass filter for raw and set l_trans_bandwidth=0.01 (the default is 1 Hz). Then, plot its "time" and "magnitude" response between 0.01 and 5 Hz. What warning does create_filter show you, and why?

Solution

filt = mne.filter.create_filter(
    raw.get_data(), raw.info["sfreq"], l_freq=1, h_freq=None, l_trans_bandwidth=0.01)
mne.viz.plot_filter( filt, raw.info["sfreq"], flim=(0.001, 5), plot=("time", "magnitude"));

Setting up high-pass filter at 1 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 1.00
- Lower transition bandwidth: 0.01 Hz (-6 dB cutoff frequency: 0.99 Hz)
- Filter length: 198203 samples (330.000 s)

Section 3: Filtering Artifacts

Background

Filtering is not a lossless operation — even when correctly applied, a filter can introduce its own distortions into the signal. Understanding these filtering artifacts is essential for correctly interpreting EEG and MEG data.

The most common filtering artifact is ringing. When a signal contains a sharp transition — such as a sudden step change or a brief pulse — the filter’s finite impulse response cannot represent this discontinuity perfectly. Instead, the filtered signal oscillates around the transition point, producing wave-like ripples that can extend well before and after the event. These are the same oscillations that cause stopband ripple in the frequency domain (see the previous section).

The severity and duration of ringing depend on filter length: a longer filter (i.e., a narrower transition band) produces more pronounced and more temporally extended ringing. This creates a practical dilemma — the filter settings chosen for better frequency selectivity are exactly those that produce the worst time-domain distortions. For a detailed discussion on filtering artifacts and how to avoid them, see De Cheveigné & Nelken, 2019.

Exercises

In this section, you are going to apply filters to simulated data (impulses and step functions) to see the resulting artifacts and how they relate to the properties of the filter. While these artificial signals differ from real neural data, they provide us with a simple scenario for testing the side effects of filtering.

Code	Description
`raw = utils.make_step(dur, t_step, plot=True)`	Create a step signal of `dur` seconds that jumps from 0 to 1 at `t_step` seconds
`raw = utils.make_impulse(dur, t_impulse, width, plot=True)`	Create a rectangular pulse of `width` seconds starting at `t_impulse` seconds
`utils.plot_filtered(before, after, tmin, tmax)`	Plot the first channel of `before` and `after` overlaid, between `tmin` and `tmax` seconds
`raw_filt = raw.copy().filter(l_freq=None, h_freq=40)`	Create a filtered copy of `raw` with a `40` Hz lowpass filter
`raw_filt = raw.copy().filter(l_freq=1, h_freq=None)`	Create a filtered copy of `raw` with a `1` Hz highpass filter
`raw_filt = raw.copy().filter(l_freq=1, h_freq=None, l_trans_bandwidth=0.5)`	Same as above, but set the lower transition bandwidth to `0.5` Hz

Execute the cell below to generate a step signal that jumps from 0 to 1 at 15 seconds after signal onset. While this is an artificial signal, certain neural events can be modeled as step functions (see, for example Southwell et al., 2017).

raw = utils.make_step(dur=60, t_step=15, plot=True)

Creating RawArray with float64 data, n_channels=1, n_times=30000
    Range : 0 ... 29999 =      0.000 ...    59.998 secs
Ready.

Example: Apply a 40 Hz lowpass filter to a copy of raw and plot the signal before and after filtering between 14 and 16 seconds (i.e. 1 second before and after the step)

raw_filt = raw.copy().filter(l_freq=None, h_freq=40)
utils.plot_filtered(before=raw, after=raw_filt, tmin=14, tmax=16)

Filtering raw data in 1 contiguous segment
Setting up low-pass filter at 40 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal lowpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Upper passband edge: 40.00 Hz
- Upper transition bandwidth: 10.00 Hz (-6 dB cutoff frequency: 45.00 Hz)
- Filter length: 165 samples (0.330 s)

Exercise: Apply a 40 Hz lowpass filter with h_trans_bandwidth=1 to a copy of raw and plot the signal before and after filtering between 14 and 16 seconds.

Solution

raw_filt = raw.copy().filter(l_freq=None, h_freq=40, h_trans_bandwidth=1)
utils.plot_filtered(before=raw, after=raw_filt, tmin=14, tmax=16)

Filtering raw data in 1 contiguous segment
Setting up low-pass filter at 40 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal lowpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Upper passband edge: 40.00 Hz
- Upper transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 40.50 Hz)
- Filter length: 1651 samples (3.302 s)

Exercise: Apply a 1 Hz highpass filter to a copy of raw and plot the signal before and after filtering between 10 and 20 seconds.

Solution

raw_filt = raw.copy().filter(l_freq=1, h_freq=None)
utils.plot_filtered(before=raw, after=raw_filt, tmin=10, tmax=20)

Filtering raw data in 1 contiguous segment
Setting up high-pass filter at 1 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 1.00
- Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)
- Filter length: 1651 samples (3.302 s)

Exercise: Apply a 1 Hz highpass filter with l_trans_bandwidth=0.06 to a copy of raw and plot the signal before and after filtering between 10 and 20 seconds. How does reducing the transition band width affect the ringing artifact?

Solution

raw_filt = raw.copy().filter(l_freq=1, h_freq=None, l_trans_bandwidth=0.06)
utils.plot_filtered(before=raw, after=raw_filt, tmin=10, tmax=20)

Filtering raw data in 1 contiguous segment
Setting up high-pass filter at 1 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 1.00
- Lower transition bandwidth: 0.06 Hz (-6 dB cutoff frequency: 0.97 Hz)
- Filter length: 27501 samples (55.002 s)

Run the cell below to generate a 1 second long pulse at 20 seconds after signal onset. While this is an artificial signal, certain neural events can be very pulse-like such as the K-complexes and vertex waves that can be seen during sleep.

raw = utils.make_impulse(dur=60, t_impulse=20, width=1, plot=True)

Creating RawArray with float64 data, n_channels=1, n_times=30000
    Range : 0 ... 29999 =      0.000 ...    59.998 secs
Ready.

Exercise: Apply a 40 Hz lowpass filter and plot the signal before and after filtering between 19 and 22 seconds.

Solution

raw_filt = raw.copy().filter(l_freq=None, h_freq=40)
utils.plot_filtered(before=raw, after=raw_filt, tmin=19, tmax=22)

Filtering raw data in 1 contiguous segment
Setting up low-pass filter at 40 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal lowpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Upper passband edge: 40.00 Hz
- Upper transition bandwidth: 10.00 Hz (-6 dB cutoff frequency: 45.00 Hz)
- Filter length: 165 samples (0.330 s)

Exercise: Apply a 0.1 Hz highpass filter to a copy of raw and plot the signal before and after filtering between 0 and 60 seconds.

Solution

raw_filt = raw.copy().filter(l_freq=0.1, h_freq=None)
utils.plot_filtered(before=raw, after=raw_filt, tmin=0, tmax=60)

Filtering raw data in 1 contiguous segment
Setting up high-pass filter at 0.1 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 0.10
- Lower transition bandwidth: 0.10 Hz (-6 dB cutoff frequency: 0.05 Hz)
- Filter length: 16501 samples (33.002 s)

Exercise: Apply a 0.1 Hz highpass filter with l_trans_bandwidth=0.06 to a copy of raw and plot the signal before and after filtering between 0 and 60 seconds. How does reducing the transition band width affect the filtered signal?

Solution

raw_filt = raw.copy().filter(l_freq=0.1, h_freq=None, l_trans_bandwidth=0.06)
utils.plot_filtered(before=raw, after=raw_filt, tmin=0, tmax=60)

Filtering raw data in 1 contiguous segment
Setting up high-pass filter at 0.1 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 0.10
- Lower transition bandwidth: 0.06 Hz (-6 dB cutoff frequency: 0.07 Hz)
- Filter length: 27501 samples (55.002 s)

Section 4: Causal and Acausal Filtering

Background

In the previous sections, we saw that filtering artifacts can extend into the time range before the (simulated) neural event. This is because, by default, MNE uses acausal filters. In acausal filters, each sample of the output depends on past and future samples of the input. In contrast, with causal filters, the output depends only on past and current samples of the input.

(from De Cheveigné & Nelken, 2019)

This has important consequences for neural data analysis: acausal filters introduce no net phase delay, preserving the exact timing of all events, which is essential for ERP analysis. However, they produce ringing artifacts that can extend into the period before the event and may be mistakenly interpreted as an anticipatory neural signal (for an example of this see Zoefel & Heil, 2013).

To prevent such misinterpretations, it is useful to test the analysis using a causal filter to see how this changes the result. MNE provides two kinds of causal FIR filters:

Linear-phase filters delay every frequency in the passband by the same amount, which is proportional to the length of the filter. Because the delay is uniform, the shape of any signal in the passband is perfectly preserved; it just arrives later.
Minimum-phase filters introduce the same total delay as linear-phase filters but distribute it unequally: most of the delay is concentrated in the stopband frequencies. The advantage is that the delay of the passband frequencies is minimised; the downside is that the signal is distorted by the non-uniform delay.

Exercises

In this section, you’ll filter the same simulated signals as in the previous section using filters with different phase responses to see how this affects the temporal relationship between input and output. This requires setting the phase parameter of the .filter() method as shown below:

Code	Description
`raw.filter(l_freq, h_freq, phase='zero')`	Apply a zero-phase (acausal) filter — MNE default; symmetric ringing, no timing shift
`raw.filter(l_freq, h_freq, phase='linear')`	Apply a linear-phase (causal) filter; signal shape preserved, uniform timing shift
`raw.filter(l_freq, h_freq, phase='minimum')`	Apply a minimum-phase (causal) filter; minimum latency, non-uniform timing shift

Example: Simulate a step function, filter it with a "zero" phase 1 Hz highpass filter and plot the result. This is the same filter used in the previous section (phase="zero" is the default in MNE).

raw = utils.make_step(dur=60, t_step=15, plot=False)
raw_filt = raw.copy().filter(l_freq=1, h_freq=None, phase="zero")
utils.plot_filtered(before=raw, after=raw_filt, tmin=10, tmax=20)

Creating RawArray with float64 data, n_channels=1, n_times=30000
    Range : 0 ... 29999 =      0.000 ...    59.998 secs
Ready.
Filtering raw data in 1 contiguous segment
Setting up high-pass filter at 1 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 1.00
- Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)
- Filter length: 1651 samples (3.302 s)

Exercise: Filter the simulated step function with a "linear" phase 1 Hz highpass filter and plot the result.

Solution

raw_filt = raw.copy().filter(l_freq=1, h_freq=None, phase="linear")
utils.plot_filtered(before=raw, after=raw_filt, tmin=10, tmax=20)

Filtering raw data in 1 contiguous segment
Setting up high-pass filter at 1 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower transition bandwidth: 1.00 Hz
- Filter length: 1651 samples (3.302 s)

Exercise: Filter the simulated step function with a "minimum" phase 1 Hz highpass filter and plot the result.

Solution

raw_filt = raw.copy().filter(l_freq=1, h_freq=None, phase="minimum")
utils.plot_filtered(before=raw, after=raw_filt, tmin=10, tmax=20)

Filtering raw data in 1 contiguous segment
Setting up high-pass filter at 1 Hz

FIR filter parameters
---------------------
Designing a one-pass, non-linear phase, causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower transition bandwidth: 1.00 Hz
- Filter length: 1651 samples (3.302 s)

Example: Simulate a square pulse, filter it using a "zero" phase 40 Hz lowpass filter and plot the result.

raw = utils.make_impulse(dur=30, t_impulse=10, width=1, plot=False)
raw_filt = raw.copy().filter(l_freq=None, h_freq=40, phase="zero")
utils.plot_filtered(before=raw, after=raw_filt, tmin=9, tmax=12)

Creating RawArray with float64 data, n_channels=1, n_times=15000
    Range : 0 ... 14999 =      0.000 ...    29.998 secs
Ready.
Filtering raw data in 1 contiguous segment
Setting up low-pass filter at 40 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal lowpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Upper passband edge: 40.00 Hz
- Upper transition bandwidth: 10.00 Hz (-6 dB cutoff frequency: 45.00 Hz)
- Filter length: 165 samples (0.330 s)

Exercise: Filter the square pulse using a "linear" phase 40 Hz lowpass filter and plot the result.

Solution

raw_filt = raw.copy().filter(l_freq=None, h_freq=40, phase="linear")
utils.plot_filtered(before=raw, after=raw_filt, tmin=9, tmax=12)

Filtering raw data in 1 contiguous segment
Setting up low-pass filter at 40 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal lowpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Upper transition bandwidth: 10.00 Hz
- Filter length: 165 samples (0.330 s)

Exercise: Filter the square pulse using a "minimum" phase 40 Hz lowpass filter and plot the result.

Solution

raw_filt = raw.copy().filter(l_freq=None, h_freq=40, phase="minimum")
utils.plot_filtered(before=raw, after=raw_filt, tmin=9, tmax=12)

Filtering raw data in 1 contiguous segment
Setting up low-pass filter at 40 Hz

FIR filter parameters
---------------------
Designing a one-pass, non-linear phase, causal lowpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Upper transition bandwidth: 10.00 Hz
- Filter length: 165 samples (0.330 s)