Frequency domain plot no different after MNE low-pass and high-pass filtering

[oskarhibbert (Oskar)]

Hi! After importing my EEG recording using:

raw = mne.io.read_raw_edf(filepath, preload=True).filter(l_freq=1, h_freq=40)
# ... lots of processing or 'raw'
# EEG = 
              FP1-F7         F7-T7  ...      FT9-FT10       FT10-T8
0       1.016440e-14 -3.388132e-15  ... -1.270549e-14 -2.752857e-15
1      -2.580525e+01 -1.245503e+01  ...  6.489457e+00  8.350538e+00
2      -4.289616e+01 -2.009946e+01  ...  1.073307e+01  1.373020e+01
3      -4.793277e+01 -2.150957e+01  ...  1.196083e+01  1.504048e+01
4      -4.432927e+01 -1.919635e+01  ...  1.115241e+01  1.352766e+01
             ...           ...  ...           ...           ...
921595  6.392553e+00 -2.020362e+01  ...  8.076414e+00  7.673939e+00
921596  5.983158e+00 -1.690656e+01  ...  2.825419e+00  7.266273e+00
921597  4.482021e+00 -1.176326e+01  ... -4.654217e-01  6.606688e+00
921598  2.356530e+00 -5.916955e+00  ... -1.160807e+00  4.144133e+00
921599  5.293956e-15  2.964615e-15  ... -2.329341e-15 -2.541099e-15

[921600 rows x 22 columns]

And applying SciPy’s Fourier transform functions to plot the frequency domain of my recording:

N = len(EEG[['FP1-F7']]) #921600
f = 256 #256 Hz

xf = fftfreq(N, 1/f)
yf = fft(EEG[['FP1-F7']])

plt.plot(xf, np.abs(yf))
plt.xlabel("Frequency")
plt.ylabel("Amplitude")
plt.title("Frequency domain (Fourier transform) of single signal")
plt.show()

I still receive a frequency domain plot ranging from -128 to +128 Hz with considerable amplitudes across all frequencies, e.g. a large spike at 60 Hz still (Imgur: The magic of the Internet).

Shouldn’t the high-pass, low-pass filtering filter out the frequencies outside the range specified (in my case: 1 to 40Hz) as per the filtering applied, more or less flattening the amplitudes outside of that range? I.e., in my plot, should the amplitudes before 0 and after 40 Hz be almost flat?

MNE version and OS information:

• MNE version: 1.2.1
• operating system: macOS 12.6

[mscheltienne (Mathieu Scheltienne)]

Hello,

It looks weird, but it’s difficult to tell what has been going wrong from those code snippets.
Can you do:

raw = mne.io.read_raw_edf(filepath, preload=True)
raw.plot_psd()
raw.filter(1., 40.)
raw.plot_psd()

And share both plots?

Mathieu

[oskarhibbert (Oskar)]

Hi Mathieu! Thank you sincerely for your reply. You’re right! The filter is applied as can be seen in the following two graphs (without and with filtering resp.):

Indeed after applying the filter I checked the EEG data-frame I am creating, and the values do change.

So I must be fundamentally mistaken how Fourier Transform works (in Python).

If I were to plot the frequency domain of one of the EEG channels after filtering, shouldn’t the plot reflect this change? I.e., shouldn’t the power in the red boxed be close to zero in the following plot?

My code:

raw = mne.io.read_raw_edf(filepath, verbose='Critical', preload=True)
raw.filter(1., 40.)
    
raw_ch_names = raw.info['ch_names']
raw_data = raw.get_data().T
eeg = pd.DataFrame(data=raw_data*10**6, columns=raw_ch_names)
eeg = eeg.loc[:, ~eeg.head(10).T.duplicated()] 

def fourier_transform(df):
    
    N = len(df)
    f = 256

    xf = fftfreq(N, 1/f)
    yf = fft(df)
    
    return xf, yf

xf, yf = fourier_transform(eeg[['FP1-F7']])

plt.plot(xf, np.abs(yf))
plt.xlabel("Frequency")
plt.ylabel("Amplitude")
plt.title("Frequency domain (Fourier transform) of single signal")
plt.show()

[mscheltienne (Mathieu Scheltienne)]

From your code, I don’t see anything that looks wrong. But I am not super fluent with pandas, I might be missing something.

Here is a short example of the EEG channels of a sample MEG recording you can play with. I’ve used the FFT from numpy directly:

import numpy as np
from matplotlib import pyplot as plt
from mne.io import read_raw_fif
from mne.datasets import sample


#%% load raw data and drop MEG channels
directory = sample.data_path() / "MEG" / "sample"
raw = read_raw_fif(directory / "sample_audvis_raw.fif", preload=True)
raw.pick_types(eeg=True, exclude="bads")

#%% high-pass to avoid the DC-offset
raw.filter(1., None)

#%% before filter
data = raw.get_data()
winsize = data.shape[-1]
frequencies = np.fft.rfftfreq(winsize, 1 / raw.info["sfreq"])
fftval = np.abs(np.fft.rfft(data, axis=1))
# average across channels
fftval = np.average(fftval, axis=0)
# plots
fig, ax = plt.subplots(2, 1)
raw.plot_psd(ax=ax[0])
ax[1].plot(frequencies, fftval)
ax[0].set_title("PSD")
ax[1].set_title("abs(FFT)")

#%% after filter
raw.filter(1., 40.)
data = raw.get_data()
winsize = data.shape[-1]
frequencies = np.fft.rfftfreq(winsize, 1 / raw.info["sfreq"])
fftval = np.abs(np.fft.rfft(data, axis=1))
# average across channels
fftval = np.average(fftval, axis=0)
# plots
fig, ax = plt.subplots(2, 1)
raw.plot_psd(ax=ax[0])
ax[1].plot(frequencies, fftval)
ax[0].set_title("PSD")
ax[1].set_title("abs(FFT)")

As you can see, the filter effect is very clear.

[oskarhibbert (Oskar)]

You’re right! This worked for me! Thank you so much.

Side-note, do you have any idea what those two large spikes might be in my data? Some noise presumably? I hate to take more of your time so please ignore if you don’t feel like / have time for responding!

[mscheltienne (Mathieu Scheltienne)]

No idea, but you should investigate. Try to look at raw.plot() and to find this activity in the temporal domain, try to figure out if it is a constant background noise or if it just happens at specific timings/intervals.

I’m going to guess an ICA should be able to isolate it well.

[drammock (Dan McCloy)]

@oskarhibbert

raw_ch_names = raw.info['ch_names']
raw_data = raw.get_data().T
eeg = pd.DataFrame(data=raw_data*10**6, columns=raw_ch_names)
eeg = eeg.loc[:, ~eeg.head(10).T.duplicated()] 

ASIDE: Just an FYI to make your life easier: there is a raw.to_data_frame() method that handles column names automatically, scales EEG from V to uV, lets you subselect specific channels or time ranges, etc. mne.io.Raw — MNE 1.6.0 documentation