incorrect error message in raw.notch_filter

[jasmainak]

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..match_module_callback at 0x7fd4f91fc820> Traceback (most recent call last): File "/autofs/space/meghnn_001/users/mjas/anaconda3/envs/mne/lib/python3.8/site-packages/threadpoolctl.py", line 400, in match_module_callback self._make_module_from_path(filepath) File "/autofs/space/meghnn_001/users/mjas/anaconda3/envs/mne/lib/python3.8/site-packages/threadpoolctl.py", line 515, in _make_module_from_path module = module_class(filepath, prefix, user_api, internal_api) File "/autofs/space/meghnn_001/users/mjas/anaconda3/envs/mne/lib/python3.8/site-packages/threadpoolctl.py", line 606, in init self.version = self.get_version() File "/autofs/space/meghnn_001/users/mjas/anaconda3/envs/mne/lib/python3.8/site-packages/threadpoolctl.py", line 646, in get_version config = get_config().split() AttributeError: 'NoneType' object has no attribute 'split' mne: 1.4.dev0 numpy: 1.23.0 {OpenBLAS 0.3.12 with 1 thread} scipy: 1.5.3 matplotlib: 3.3.3 {backend=Qt5Agg}

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

[sena-neuro]

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

[larsoner]

@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?

[sena-neuro]

Sure!