Filtering in jumeg

Comparison of the various filters availables for use in mne/jumeg



In [1]:

    
%matplotlib inline



In [4]:

    
import numpy as np
import matplotlib.pyplot as plt

import mne
from mne.datasets import sample

data_path = sample.data_path()
raw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif'

# load the raw data
raw = mne.io.Raw(raw_fname, preload=True)
picks = mne.pick_types(raw.info, meg='mag', exclude='bads')
print('Sampling frequency is %f' % raw.info['sfreq'])









    



Opening raw data file /Users/psripad/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.
Current compensation grade : 0
Reading 0 ... 166799  =      0.000 ...   277.714 secs...
Sampling frequency is 600.614990



In [5]:

    
# plot the raw psd
raw.plot_psd(tmin=0.0, tmax=160., fmin=0, fmax=np.inf, n_fft=None,
             picks=picks, ax=None, color='black', area_mode='std',
             area_alpha=0.33, n_overlap=0, dB=True, average=True,
             show=True, n_jobs=4, line_alpha=None,
             spatial_colors=None, xscale='linear', verbose='ERROR');



In [6]:

    
# filter data using mne filter_data
# (mne.filter.band_pass_filter, jumeg.jumeg_preprocessing.apply_filter all use the same underlying function)

l_freq, h_freq = 1., 45.
fir_filt = raw.copy().filter(l_freq, h_freq, picks=picks, filter_length='auto',
                             l_trans_bandwidth='auto', h_trans_bandwidth='auto',
                             n_jobs=4, method='fir', iir_params=None, phase='zero',
                             fir_window='hamming', verbose='ERROR')

# plot the FIR filtered raw psd
fig = fir_filt.plot_psd(tmin=0.0, tmax=160., fmin=0, fmax=70., n_fft=None,
                  picks=picks, ax=None, color='black', area_mode='std',
                  area_alpha=0.33, n_overlap=0, dB=True, average=True,
                  show=False, n_jobs=4, line_alpha=None,
                  spatial_colors=None, xscale='linear', verbose='ERROR');
fig.gca().set_title('MNE FIR filter');



In [7]:

    
filter_type = 'butter'
filt_method = 'fft'
iir_params={'ftype': filter_type, 'order': 4}

# apply IIR filter
iir_filt = raw.copy().filter(l_freq, h_freq, picks=picks, filter_length='auto',
                             l_trans_bandwidth='auto', h_trans_bandwidth='auto',
                             n_jobs=4, method='iir', iir_params=iir_params, phase='zero',
                             fir_window='hamming', verbose='ERROR')

# plot the IIR filtered raw psd
fig = iir_filt.plot_psd(tmin=0.0, tmax=160., fmin=l_freq, fmax=h_freq, n_fft=None,
                  picks=picks, ax=None, color='black', area_mode='std',
                  area_alpha=0.33, n_overlap=0, dB=True, average=True,
                  show=False, n_jobs=4, line_alpha=None,
                  spatial_colors=None, xscale='linear', verbose='ERROR');
fig.gca().set_title('MNE BW IIR filter');



In [8]:

    
from jumeg.filter import jumeg_filter



In [10]:

    
# apply the jumeg filter using mne (just to check)
filt_ju_mne = jumeg_filter(filter_method='mne', filter_type='bp', fcut1=l_freq, fcut2=h_freq,
                           remove_dcoffset=True, sampling_frequency=raw.info['sfreq'],
                           filter_window='hamming', notch=np.array([50., 60.]), notch_width=1.0, order=4, njobs=4,
                           mne_filter_method='fft',mne_filter_length='10s',trans_bandwith=0.5)
ju_mne_filt = raw.copy()  # make a copy
ju_mne_filt._data = filt_ju_mne.apply_filter(ju_mne_filt._data, picks)

# plot the jumeg MNE filtered raw psd
fig = ju_mne_filt.plot_psd(tmin=0.0, tmax=160., fmin=0, fmax=70., n_fft=None,
                     picks=picks, ax=None, color='black', area_mode='std',
                     area_alpha=0.33, n_overlap=0, dB=True, average=True,
                     show=False, n_jobs=4, line_alpha=None,
                     spatial_colors=None, xscale='linear', verbose=None);
fig.gca().set_title('jumeg filter using MNE FIR/fft');









    



Effective window size : 3.410 (s)






    



[Parallel(n_jobs=4)]: Using backend LokyBackend with 4 concurrent workers.
[Parallel(n_jobs=4)]: Done   2 out of   4 | elapsed:    0.4s remaining:    0.4s
[Parallel(n_jobs=4)]: Done   4 out of   4 | elapsed:    0.7s remaining:    0.0s
[Parallel(n_jobs=4)]: Done   4 out of   4 | elapsed:    0.7s finished



In [16]:

    
# apply the jumeg filter using bw
filt_ju_bw = jumeg_filter(filter_method="bw", filter_type='bp', fcut1=l_freq, fcut2=h_freq,
                          remove_dcoffset=True, sampling_frequency=raw.info['sfreq'],
                          filter_window='hamming', notch=np.array([50., 60.]), notch_width=1.0, order=4, njobs=4,
                          mne_filter_method='fft',mne_filter_length='10s',trans_bandwith=0.5)
filt_ju_bw.verbose = False
ju_bw_filt = raw.copy()  # make a copy
filt_ju_bw.apply_filter(ju_bw_filt._data, picks)

# plot the jumeg MNE filtered raw psd
fig = ju_bw_filt.plot_psd(tmin=0.0, tmax=160., fmin=0.1, fmax=70., n_fft=None,
                    picks=picks, ax=None, color='black', area_mode='std',
                    area_alpha=0.33, n_overlap=0, dB=True, average=True,
                    show=False, n_jobs=4, line_alpha=None,
                    spatial_colors=None, xscale='linear', verbose=None);
fig.gca().set_title('Jumeg BW filter');









    



Effective window size : 3.410 (s)






    



[Parallel(n_jobs=4)]: Using backend LokyBackend with 4 concurrent workers.
[Parallel(n_jobs=4)]: Done   2 out of   4 | elapsed:    1.2s remaining:    1.2s
[Parallel(n_jobs=4)]: Done   4 out of   4 | elapsed:    2.1s remaining:    0.0s
[Parallel(n_jobs=4)]: Done   4 out of   4 | elapsed:    2.1s finished



In [12]:

    
print(filt_ju_bw.filter_name_postfix)









    



fibp1-45n2



In [17]:

    
## plot the window sinc filter results also
# from jumeg.filter.jumeg_filter_ws import JuMEG_Filter_Ws
# filt_ju_ws = jumeg_filter(filter_method='ws', filter_type='bp', fcut1=l_freq, fcut2=h_freq,
#                           remove_dcoffset=True, sampling_frequency=raw.info['sfreq'],
#                           filter_window='hamming', notch=np.array([50., 60.]), notch_width=1.0, order=4, njobs=4,
#                           mne_filter_method='fft',mne_filter_length='10s',trans_bandwith=0.5)

# ju_ws_filt = raw.copy()  # make a copy
# ju_ws_filt.notch_filter(np.array([50., 60.]), picks=picks, filter_length='auto',
#                         notch_widths=None, trans_bandwidth=1.0, n_jobs=1, method='fir',
#                         iir_params=None, mt_bandwidth=None, p_value=0.05, phase='zero',
#                         fir_window='hamming', verbose=None)

# filt_ju_ws.apply_filter(ju_ws_filt._data, picks)

## plot the jumeg Window Sinc filtered raw psd
# fig = ju_ws_filt.plot_psd(tmin=0.0, tmax=160., fmin=0, fmax=70., n_fft=None,
#                     picks=picks, ax=None, color='black', area_mode='std',
#                     area_alpha=0.33, n_overlap=0, dB=True, average=True,
#                     show=False, n_jobs=4, line_alpha=None,
#                     spatial_colors=None, xscale='linear', verbose=None);
# fig.gca().set_title('Jumeg Window Sinc filter');









    



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

Effective window size : 3.410 (s)






    



[Parallel(n_jobs=4)]: Using backend LokyBackend with 4 concurrent workers.
[Parallel(n_jobs=4)]: Done   2 out of   4 | elapsed:    0.5s remaining:    0.5s
[Parallel(n_jobs=4)]: Done   4 out of   4 | elapsed:    0.7s remaining:    0.0s
[Parallel(n_jobs=4)]: Done   4 out of   4 | elapsed:    0.7s finished



In [27]:

    
plt.close()

fig, ax1 = plt.subplots(1, 1)

# plot the jumeg MNE filtered raw psd
ju_mne_filt.plot_psd(tmin=0.0, tmax=160., fmin=0, fmax=70., n_fft=None,
                     picks=picks, ax=ax1, color='blue', area_mode='std',
                     area_alpha=0.33, n_overlap=0, dB=True, average=True,
                     show=False, n_jobs=4, line_alpha=None,
                     spatial_colors=None, xscale='linear', verbose=None);

# plot the jumeg Butterworth filtered raw psd
ju_bw_filt.plot_psd(tmin=0.0, tmax=160., fmin=0, fmax=70., n_fft=None,
                    picks=picks, ax=ax1, color='green', area_mode='std',
                    area_alpha=0.33, n_overlap=0, dB=True, average=True,
                    show=False, n_jobs=4, line_alpha=None,
                    spatial_colors=None, xscale='linear', verbose=None);

# # plot the jumeg Windowsinc filtered raw psd
# fig = ju_ws_filt.plot_psd(tmin=0.0, tmax=160., fmin=0, fmax=70., n_fft=None,
#                     picks=picks, ax=ax1, color='red', area_mode='std',
#                     area_alpha=0.33, n_overlap=0, dB=True, average=True,
#                     show=True, n_jobs=4, line_alpha=None,
#                     spatial_colors=None, xscale='linear', verbose=None);

ax1.set_xlim(0., 75.);
ax1.set_title('MNE (blue)/jumeg BW (green)');









    



Effective window size : 3.410 (s)






    



[Parallel(n_jobs=4)]: Using backend LokyBackend with 4 concurrent workers.
[Parallel(n_jobs=4)]: Done   2 out of   4 | elapsed:    1.1s remaining:    1.1s






    



Effective window size : 3.410 (s)






    



[Parallel(n_jobs=4)]: Done   4 out of   4 | elapsed:    2.1s remaining:    0.0s
[Parallel(n_jobs=4)]: Done   4 out of   4 | elapsed:    2.1s finished
[Parallel(n_jobs=4)]: Using backend LokyBackend with 4 concurrent workers.
[Parallel(n_jobs=4)]: Done   2 out of   4 | elapsed:    1.5s remaining:    1.5s
[Parallel(n_jobs=4)]: Done   4 out of   4 | elapsed:    3.3s remaining:    0.0s
[Parallel(n_jobs=4)]: Done   4 out of   4 | elapsed:    3.3s finished

Various available filters are applied on raw MEG data of sampling frequency 1017.25 above.

The mne filter uses filter_data routine to perform band pass filtering and the notch_filter routine to perform notch filtering.
MNE filter with FIR filtering performs better than using an IIR Butterworth filter with oder 4.
Jumeg filter module has three filters implemented
- MNE filter (which wraps around MNE modules) (can be combined with notch filters)
- jumeg Window Sinc filter (WS) (notch filter combination not implemented)
- jumeg Butterworth filter (BW) (can be combined with notch filters)
Of all the three jumeg implementations above, mne FIR filter shows the best performance.
The WS implementation performs better than the jumeg BW (with order 4).

Recommendation: Presently, it is best to use the MNE FIR filter for our filtering requirements. The jumeg filter module (which wraps around mne filter) may be used when the MNE band_pass_filter needs to be combined with the notch_filter. In cases where the notch frequencies are removed using the noise reducer, the mne FIR filter can be directly used.