#!/usr/bin/env python
# coding: utf-8

'''
Filtering in jumeg.

Comparison of the various filters available for use in mne/jumeg.
'''

import numpy as np
import matplotlib.pyplot as plt

import mne
from mne.datasets import sample

from jumeg.filter import jumeg_filter

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

# plot the raw psd
fig = 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=False, n_jobs=4, line_alpha=None,
                   spatial_colors=None, xscale='linear', verbose='ERROR');
fig.suptitle('Original sprectrum')

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

# apply FIR filter
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')

# apply IIR filter
filter_type = 'butter'
filt_method = 'fft'
iir_params={'ftype': filter_type, 'order': 4}

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

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

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

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

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

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

fir_filt.plot_psd(tmin=0.0, tmax=160., fmin=0, fmax=70., n_fft=None,
                  picks=picks, ax=ax2, color='red', 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');

iir_filt.plot_psd(tmin=0.0, tmax=160., fmin=0., fmax=70., n_fft=None,
                  picks=picks, ax=ax2, color='yellow', 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');
ax2.set_xlim(0., 75.);
ax2.set_title('MNE FIR (red)/ MNE IIR (yellow)');

fig.tight_layout()
plt.show()

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

@@ -942,156 +942,8 @@
 ')%0A%0A
-# filter data using mne filter_data%0A# (mne.filter.band_pass_filter, jumeg.jumeg_preprocessing.apply_filter%0A# all use the same underlying function)%0A%0A
 l_fr
@@ -1406,16 +1406,35 @@
 r': 4%7D%0A%0A
+# apply IIR filter%0A
 iir_filt
@@ -1783,24 +1783,8 @@
  mne
- (just to check)
 %0Afil
@@ -2986,16 +2986,38 @@
 meg MNE 
+and jumeg Butterworth 
 filtered
@@ -3020,24 +3020,25 @@
 ered raw psd
+s
 %0Aju_mne_filt
@@ -3384,54 +3384,8 @@
 );%0A%0A
-# plot the jumeg Butterworth filtered raw psd%0A
 ju_b
@@ -3806,16 +3806,63 @@
 en)');%0A%0A
+# plot the MNE FIR and IIR filterered raw psds%0A
 fir_filt

Add comments.