1
$\begingroup$

I am analyzing functional MRI data with a sampling rate of 1 second (1 Hz). The frequency band that I am interested in is 0.01-0.2 Hz.

Regarding measurements, I am interested in computing (1) the signal’s temporal autocorrelation in the time-domain, and (2) the slope of a linear least-square fit in the log-log transformed frequency-domain.

Because temporal autocorrelation drastically increases to unrealistic values without bandpassing, I would like to bandpass the data.

Problem: My problem is that a bandpass filter that I constructed via the code below in Python introduces relatively heavy drops in power near the lowest (0.01 Hz) and highest (0.2 Hz) frequency that destroys the scale-free nature of the data.

# data = time-series data that I load for several subjects

sr = 1/sr # s to Hz
low = 0.01 # min freq
high = 0.2 # max freq

bandpass = sp.signal.butter(N=5, fs=sr, Wn=[low, high], btype="bandpass", analog=False, output="sos")
filtered_data = sp.signal.sosfilt(sos=bandpass, x=data, zi=None)

Within a neuroimaging preprocessing software, I have previously used the same bandpass range, that is, 0.01-0.2 Hz. Interestingly, the results of this neuroimaging software do not produce a loss in power near the lower and upper end of the frequency band, as observed here for sp.signal.butter. I noticed that increasing the order (=N) option of sp.signal.butter, such as from 1 to 5, 10 or even higher slightly reduces this drop in power near the edges, but it does not disappear.

A comparison between non-bandpassed and bandpassed data (using scipy) is shown below. enter image description here enter image description here

Question: Is there a way to modify sp.signal.butter or another way to bandpass data in Python without decreasing the power of the signal near the lower and upper frequency limits of the bandpass filter?

Update: Instead of using sp.signal.butter, I tested a Chebyshev Type II filter via sp.signal.cheby2 and the following settings:

# Bandpass
                sr = 1/sr # s to Hz
                low = 0.01 # min freq
                high = 0.2 # max freq
                bandpass = sp.signal.cheby2(N=1, rs=0.001, Wn=[low, high],
                                            btype="bandpass", analog=False,
                                            output="sos", fs=sr)
                row = sp.signal.sosfilt(sos=bandpass, x=row, zi=None)

This results in the following output. As can be seen, the roll-off now appears to be sharp, that is, the power no longer drops near the lowest and highest frequency. I am not sure if the selected paramters are reasonable from a methodological point of view. At least visually the results look good.

enter image description here

$\endgroup$
2
  • 1
    $\begingroup$ Regarding your update, I don't think you'll achieve the roll-off you want with a 1st order, but more importantly, rs is in decibels, so you're specifying $0.001\tt{dB}$ minimum attenuation in the stopband... basically you have an all-pass filter. $\endgroup$
    – Jdip
    Jul 18, 2023 at 14:13
  • $\begingroup$ You are right and I just realized that. The power spectra without bandpassing and the spectra with rs=0.001 look identical. Also, the temporal correlation with rs=0.001 is way too high, which makes sense given that this filter does not filter out anything. I tried your suggested broader frequency band, that is, Wn=[0.006, 0.27]. But even that results in a loss of the scale-free power spectrum. Probably it is best only to apply bandpassing for the time-domain, but not for the frequency-domain (since I "cut" the power spectra to the desired band by sp.periodogram anyway). $\endgroup$
    – Philipp
    Jul 18, 2023 at 14:19

2 Answers 2

2
$\begingroup$

The cut-off frequencies for Butterworth filters are at $-3\tt{dB}$:

For a Butterworth filter, this is the point at which the gain drops to 1/sqrt(2) that of the passband (the “-3 dB point”).

If you want the gain at these frequencies to be closer to $0\tt{dB}$, extend the pass-band a little:
You might want to experiment a little until you're happy with the gain at the frequencies of interest

enter image description here


Code:

sr = 1

bandpass1 = signal.butter(N=5, fs=sr, Wn=[0.01, 0.2], btype="bandpass", analog=False, output="sos")
[w1, h1] = signal.sosfreqz(bandpass1, 512, fs=sr)
bandpass2 = signal.butter(N=5, fs=sr, Wn=[0.006, 0.27], btype="bandpass", analog=False, output="sos")
[w2, h2] = signal.sosfreqz(bandpass2, 512, fs=sr)

plt.semilogx(w1, 20*np.log10(np.abs(h1)))
plt.semilogx(w2, 20*np.log10(np.abs(h2)))
plt.grid(True)
plt.ylim(-3,0)
plt.xlim(0.005,0.3)
plt.xticks([0.01, 0.2], ['0.01', '0.2'] )
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude response (dB)')
plt.legend(['original', 'extended bandpass'])
plt.show()
$\endgroup$
1
  • 1
    $\begingroup$ Thank you. By increasing the bandpass filter, that is, from 0.01-0.1 to 0.005-0.25, the visible roll-off is indeed not so bad. However, the change in the high-pass filter (0.005 instead of 0.01) drastically increases the signal’s autocorrelation. I would like to use the exact same filter settings for both (1) temporal autocorrelation and (2) frequency-domain measurements (to keep both measurements compareable). Moreover, I now tested a Chebyshev type II filter, the results are much better than with the butterworth filter. I will update my initial post. $\endgroup$
    – Philipp
    Jul 18, 2023 at 13:41
1
$\begingroup$

A Butterworth filter has cutoff frequency gain of -3dB, which is the effect that you are seeing. You can either increase the band limits or use a different type of filter (Chebyshev Type 1, Elliptic), where you can specify the gain at the cutoff directly.

If you want to reduce the out-of-band frequencies more, you can increase the order of the filter and/or the stop-band attenuation (for an Elliptic filter). However, the steeper the filter becomes, the more it will smear out the signal in time and the longer it takes the filtered signal to reach its final state. That's a fundamental tradeoff of filter design.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.