mne-python
mne-python copied to clipboard
mne-tools
incorrect error message in raw.notch_filter
Description of the problem
When I try to fiddle with filter_length MNE gives an incorrect error message
Steps to reproduce
import numpy as np
import mne
sample_data_folder = mne.datasets.sample.data_path()
sample_data_raw_file = (sample_data_folder / 'MEG' / 'sample' /
'sample_audvis_filt-0-40_raw.fif')
raw = mne.io.read_raw_fif(sample_data_raw_file, preload=True)
raw.crop(tmax=1.)
raw.notch_filter(60, trans_bandwidth=10., filter_length='330ms')
Link to data
No response
Expected results
It should say that the the transition bandwidth is 10 Hz or 5 Hz (if only considering upper/lower) but it reports 2.5 Hz.
Actual results
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
File /autofs/homes/001/mj513/Desktop/mne_bug.py:10
7 raw = mne.io.read_raw_fif(sample_data_raw_file, preload=True)
9 raw.crop(tmax=1.)
---> 10 raw.notch_filter(60, trans_bandwidth=10., filter_length='330ms')
File <decorator-gen-229>:12, in notch_filter(self, freqs, picks, filter_length, notch_widths, trans_bandwidth, n_jobs, method, iir_params, mt_bandwidth, p_value, phase, fir_window, fir_design, pad, skip_by_annotation, verbose)
File /autofs/space/meghnn_001/users/mjas/github_repos/mne-python/mne/io/base.py:1065, in BaseRaw.notch_filter(self, freqs, picks, filter_length, notch_widths, trans_bandwidth, n_jobs, method, iir_params, mt_bandwidth, p_value, phase, fir_window, fir_design, pad, skip_by_annotation, verbose)
1062 logger.info('Filtering raw data in %d contiguous segment%s'
1063 % (len(onsets), _pl(onsets)))
1064 for si, (start, stop) in enumerate(zip(onsets, ends)):
-> 1065 self._data = notch_filter(
1066 self._data[:, start:stop], fs, freqs,
1067 filter_length=filter_length, notch_widths=notch_widths,
1068 trans_bandwidth=trans_bandwidth, method=method,
1069 iir_params=iir_params, mt_bandwidth=mt_bandwidth,
1070 p_value=p_value, picks=picks, n_jobs=n_jobs, copy=False,
1071 phase=phase, fir_window=fir_window, fir_design=fir_design,
1072 pad=pad)
1073 return self
File <decorator-gen-145>:12, in notch_filter(x, Fs, freqs, filter_length, notch_widths, trans_bandwidth, method, iir_params, mt_bandwidth, p_value, picks, n_jobs, copy, phase, fir_window, fir_design, pad, verbose)
File /autofs/space/meghnn_001/users/mjas/github_repos/mne-python/mne/filter.py:1187, in notch_filter(x, Fs, freqs, filter_length, notch_widths, trans_bandwidth, method, iir_params, mt_bandwidth, p_value, picks, n_jobs, copy, phase, fir_window, fir_design, pad, verbose)
1183 lows = [freq - nw / 2.0 - tb_2
1184 for freq, nw in zip(freqs, notch_widths)]
1185 highs = [freq + nw / 2.0 + tb_2
1186 for freq, nw in zip(freqs, notch_widths)]
-> 1187 xf = filter_data(x, Fs, highs, lows, picks, filter_length, tb_2, tb_2,
1188 n_jobs, method, iir_params, copy, phase, fir_window,
1189 fir_design, pad=pad)
1190 elif method == 'spectrum_fit':
1191 xf = _mt_spectrum_proc(x, Fs, freqs, notch_widths, mt_bandwidth,
1192 p_value, picks, n_jobs, copy, filter_length)
File <decorator-gen-143>:12, in filter_data(data, sfreq, l_freq, h_freq, picks, filter_length, l_trans_bandwidth, h_trans_bandwidth, n_jobs, method, iir_params, copy, phase, fir_window, fir_design, pad, verbose)
File /autofs/space/meghnn_001/users/mjas/github_repos/mne-python/mne/filter.py:812, in filter_data(data, sfreq, l_freq, h_freq, picks, filter_length, l_trans_bandwidth, h_trans_bandwidth, n_jobs, method, iir_params, copy, phase, fir_window, fir_design, pad, verbose)
810 data = _check_filterable(data)
811 iir_params, method = _check_method(method, iir_params)
--> 812 filt = create_filter(
813 data, sfreq, l_freq, h_freq, filter_length, l_trans_bandwidth,
814 h_trans_bandwidth, method, iir_params, phase, fir_window, fir_design)
815 if method in ('fir', 'fft'):
816 data = _overlap_add_filter(data, filt, None, phase, picks, n_jobs,
817 copy, pad)
File <decorator-gen-144>:12, in create_filter(data, sfreq, l_freq, h_freq, filter_length, l_trans_bandwidth, h_trans_bandwidth, method, iir_params, phase, fir_window, fir_design, verbose)
File /autofs/space/meghnn_001/users/mjas/github_repos/mne-python/mne/filter.py:1062, in create_filter(data, sfreq, l_freq, h_freq, filter_length, l_trans_bandwidth, h_trans_bandwidth, method, iir_params, phase, fir_window, fir_design, verbose)
1059 raise ValueError('Stop bands are not sufficiently '
1060 'separated.')
1061 if method == 'fir':
-> 1062 out = _construct_fir_filter(sfreq, freq, gain, filter_length, phase,
1063 fir_window, fir_design)
1064 return out
File /autofs/space/meghnn_001/users/mjas/github_repos/mne-python/mne/filter.py:385, in _construct_fir_filter(sfreq, freq, gain, filter_length, phase, fir_window, fir_design)
383 h = minimum_phase(h)
384 else:
--> 385 h = fir_design(N, freq, gain, window=fir_window)
386 assert h.size == N
387 att_db, att_freq = _filter_attenuation(h, freq, gain)
File /autofs/space/meghnn_001/users/mjas/github_repos/mne-python/mne/filter.py:304, in _firwin_design(N, freq, gain, window, sfreq)
302 this_N += (1 - this_N % 2) # make it odd
303 if this_N > N:
--> 304 raise ValueError('The requested filter length %s is too short '
305 'for the requested %0.2f Hz transition band, '
306 'which requires %s samples'
307 % (N, transition * sfreq / 2., this_N))
308 # Construct a lowpass
309 this_h = firwin(this_N, (prev_freq + this_freq) / 2.,
310 window=window, pass_zero=True, fs=freq[-1] * 2)
ValueError: The requested filter length 51 is too short for the requested 2.50 Hz transition band, which requires 99 samples
Additional information
Platform: Linux-3.10.0-1160.76.1.el7.x86_64-x86_64-with-glibc2.10 Python: 3.8.6 | packaged by conda-forge | (default, Oct 7 2020, 19:08:05) [GCC 7.5.0] Executable: /autofs/space/meghnn_001/users/mjas/anaconda3/envs/mne/bin/python3.8 CPU: x86_64: 64 cores Memory: 125.4 GB
Exception ignored on calling ctypes callback function: <function _ThreadpoolInfo._find_modules_with_dl_iterate_phdr.
sklearn: 0.23.2 numba: Not found nibabel: 3.2.1 nilearn: 0.7.0 dipy: 1.3.0 openmeeg: Not found cupy: Not found pandas: 1.1.5 pyvista: 0.36.1 {OpenGL 4.5.0 NVIDIA 455.45.01 via Quadro P5000/PCIe/SSE2} pyvistaqt: 0.9.0 ipyvtklink: Not found vtk: 9.0.1 qtpy: 2.0.1 {PyQt5=5.12.9} ipympl: Not found /autofs/space/meghnn_001/users/mjas/anaconda3/envs/mne/lib/python3.8/site-packages/pyqtgraph/colors/palette.py:1: RuntimeWarning: PyQtGraph supports Qt version >= 5.15, but 5.12.9 detected. from ..Qt import QtGui pyqtgraph: 0.13.1 pooch: v1.5.2
mne_bids: Not found mne_nirs: Not found mne_features: Not found mne_qt_browser: 0.4.0 mne_connectivity: Not found mne_icalabel: Not found
Hi everyone,
I think the issue stems from the _firwin_design function, specifically these lines:
transition = (prev_freq - this_freq) / 2.0
this_N = int(round(_length_factors[window] / transition))
this_N += 1 - this_N % 2 # make it odd
if this_N > N:
raise ValueError(
f"The requested filter length {N} is too short for the requested "
f"{transition * sfreq / 2.0:0.2f} Hz transition band, which "
f"requires {this_N} samples"
)
prev_freq and this_freq are normalized by the Nyquist frequency, as seen in _construct_fir_filter() in filter.py at line 485:
# normalize frequencies
freq = np.array(freq) / (sfreq / 2.0)
The line transition = (prev_freq - this_freq) / 2.0 seems unusual. However, this_N is still calculated correctly. This is because in the equation transition width = 3.3/N (where 3.3 is for the Hamming window), the transition is actually normalized by the sampling frequency and not the Nyquist frequency [refer to the calculation at the end of page 182]. When the transition width is normalized with the Nyquist frequency, the equation should be transition width = 6.6/N [reference].
Therefore, prev_freq - this_freq is normalized by the Nyquist frequency, but the equation needs normalization with the sampling frequency, which is why we divide by 2.
We can modify the code as follows:
transition = prev_freq - this_freq # transition width normalized to Nyquist frequency
transition_norm_fs = transition / 2
this_N = int(round(_length_factors[window] / transition_norm_fs))
this_N += 1 - this_N % 2 # make it odd
if this_N > N:
raise ValueError(
f"The requested filter length {N} is too short for the requested "
f"{transition * sfreq / 2.0:0.2f} Hz transition band, which "
f"requires {this_N} samples"
)
or
transition = (prev_freq - this_freq) / 2.0 # normalization from Nyquist freq. to sampling freq.
this_N = int(round(_length_factors[window] / transition))
this_N += 1 - this_N % 2 # make it odd
if this_N > N:
raise ValueError(
f"The requested filter length {N} is too short for the requested "
f"{transition * sfreq:0.2f} Hz transition band, which "
f"requires {this_N} samples"
)
@sena-neuro that sounsd like it's probably correct (and the second one looks better), would you be up for making a PR to fix it?
