htwangtw
commented
1 year ago The TR used in the script is too small to do meaningful filtering with those parameters. It will not create cosine drift term in the first place, so the result you got for the cosine method was only detrended. If we change the TR to 1.49, we get pretty consistent results:
Note that the consine method doesn't perform low pass filtering.
We definitely want to raise a warning if no cosine drift term is created. As for what should be the expected behaviour for the rest, I would suggest to have someone better at signal processing to weigh in
bthirion
tsalo
ymzayek
Is there an existing issue for this?
Operating system
Operating system version
Ubuntu 20.04
Python version
nilearn version
0.10.1.dev
Expected behavior
There seems to be unexpected behavior of filtering with
signal.clean. This was pointed out by a colleague (@demianw) who implemented his own filtering close to cosine. I generated some data and plots to reproduce the problem (see the code in the last section to reproduce the plots). Here I present 4 plots. The first is the original uncleaned signal. The second is an implementation by Demian and the third is nilearn cosine filtering. It is expected that the second and third plot look similarly but it looks like there is no filtering being done in the third. The last plot is basically what nilearn does for butterworth but reproduced with the underlying scipy functions and it leaves no signal. @htwangtw any ideas?Current behavior & error messages
Steps and code to reproduce bug
Code to reproduce plots
```python # Paste your code here import numpy as np import matplotlib.pyplot as plt from scipy import signal as sp_signal, fft from nilearn import signal as nisignal # create noisy data with linear trend rng = np.random.RandomState(0) npoints = 1000 noise = rng.standard_normal(npoints) x = 3 + 2*np.linspace(0, 1, npoints) + noise low_pass = 0.1 high_pass = 0.01 t_r = 1/26 # filtering implementation by Demian s = sp_signal.detrend(x) fs = fft.fft(s) t_r = 1/26 freqs = fft.fftfreq(s.size, d=t_r) freqs_pass = (np.abs(freqs) <= low_pass) * (np.abs(freqs) >= high_pass) fs_filtered = fs * freqs_pass s_filtered = fft.ifft(fs_filtered) x_filtered = np.real(s_filtered) # nilearn cosine filtering cleaned_signal = nisignal.clean( x[:, None], detrend=True, filter="cosine", low_pass=low_pass, # not used for cosine high_pass=high_pass, t_r=t_r, ) # scipy butterworth bandpass filtering (basically what nilearn does I think) order = 5 # same as nilearn default sampling_rate = 1/t_r # same as defined in nilearn b, a = sp_signal.butter(order, [high_pass, low_pass], 'band', fs=sampling_rate) filteredBandPass = sp_signal.filtfilt(b, a, x) # plot original and cleaned signals fig, axs = plt.subplots(4) axs[0].set_title("Original signal") axs[0].plot(x) axs[1].set_title("Filtering by Demian") axs[1].plot(x_filtered) axs[2].set_title("Nilearn Cosine") axs[2].plot(cleaned_signal) axs[3].set_title("Scipy butterworth with nilearn defaults") axs[3].plot(filteredBandPass) plt.show() ```