Potential problem with causal filters

[britta-wstnr]

Description of the problem

While filtering a rather long recording, I came across (to me) unexpected behavior in causal filters: they do not remove the DC offset for me. I compared the behavior to the zero-phase filter and there the results are as expected. (Baseline correction is not an option for my use case.)

I discussed this briefly with @larsoner already, who suggested to test on a synthetic data example:

Steps to reproduce

import mne
import numpy as np

data = np.zeros((100, 10000))
data[:] = np.linspace(-1, 1, 100)[:, np.newaxis]

info = mne.create_info(100, 1000, ch_types='mag')

# causal filter
evoked = mne.EvokedArray(data, info, tmin=0.0)
evoked.filter(0.5, None, phase='minimum')
evoked.plot();

# non-causal filter
evoked = mne.EvokedArray(data, info, tmin=0.0)
evoked.filter(0.5, None)
evoked.plot();

Expected results

The zero-phase non-causal filter behaves as expected:

Actual results

The causal minimum phase filter does not behave as epxected:

This is in line with what I observe on actual MEG data.

Additional information

I am in a fresh conda env and on main with scipy version 1.11.4.
@larsoner predicts this to be a problem with scipy.signal.minimum_phase and that is indeed the difference in our filter.py code between causal and non-causal filters.
I can look into this more - but could use some guidance on how to decide how/where a fix would be necessary (if this is indeed in need of a fix).

[cbrnr]

Can you post the filter design specs for both causal and non-causal filters (including amplitude and phase responses and/or impulse responses)? Maybe there's something wrong with our design (e.g. filter order too small to achieve the desired properties etc.)? It might also be that our input (e.g. the linear-phase filter) does not fulfill the requirements (e.g. the docs mention an odd number of taps)?

If the only difference is scipy.signal.minimum_phase, it might also be worth trying to design the filter with only SciPy (it could be a bug in this function). It might also be worth trying to replicate the example in the SciPy docs, and see how it behaves when you add more and more MNE stuff (starting with your artifical toy data maybe).

[britta-wstnr]

Thanks @cbrnr - I finally had time to get back to this and this is what I found:

1. Filter characteristics are different between minimum and zero phase filters in MNE

2. This can be reproduced using scipy only

(still using MNE for plotting, also please note that I did not use double the filter length for the minimum phase filter as MNE does -- to be very pedantic, I actually used double for the zero phase, as these were the numbers I had on hand.)

3. The problem seems to be in scipy - filters can be assimilated

(This is using MNE-Python for filter design again, just like 1.)

However,

this does weirdly enough not fully fix the output for me - I don't quite understand why yet.

The change

I commented out the following line of scipy.signal.minimum_phase: https://github.com/scipy/scipy/blob/main/scipy/signal/_fir_filter_design.py#L1281

From the docstring of the function I gather that @larsoner followed the recipe in this link for implementation: http://dspguru.com/dsp/howtos/how-to-design-minimum-phase-fir-filters

However, the addition of this scaling seems to be questionable, as discussed here: https://www.dsprelated.com/showthread/comp.dsp/20912-1.php

I further compared with the filter implemented in FieldTrip - which does not have any scaling added in this step either.

[cbrnr]

Thanks @britta-wstnr for investigating. I agree that the implementation might be incorrect, the best way to make sure is to directly consult e.g. Oppenheim/Schafer.

Regarding the plots, I'm not really sure if the amplitude plot explains the totally weird filter result. It does attenuate below 4 Hz, but much less than the zero-phase filter. This, however, could be due to filtering twice in the latter.

In addition, the impulse response is not visible - what is happening here with the minimum phase filter? Also, why is the impulse shortly after second 0.8 for the zero-phase filter? Should the impulse not be exactly at zero (hence zero-phase)?

[cbrnr]

PS: It looks like Fieldtrip uses the "hilbert" method of converting to minimum phase (and you are referring to "homomorphic"), did you try this with SciPy?

https://github.com/fieldtrip/fieldtrip/blob/master/preproc/private/minphaserceps.m

[britta-wstnr]

Hi @cbrnr, thanks for checking this out.
I did try the Hilbert implementation early on with not much success, but have not compared since. I think the FieldTrip implementation is the homomorphic filter approach, too?

[cbrnr]

I think the FieldTrip implementation is the homomorphic filter approach, too?

Yes, you are right, but it is missing the log (which as you and others have noted might be the correct thing to do)!

Still, it would be interesting to know if the Hilbert method works.

[britta-wstnr]

I don't think the FieldTrip implementation is missing a log - it is only missing the 0.5 (see line 57, https://github.com/fieldtrip/fieldtrip/blob/master/preproc/private/minphaserceps.m#L57). In fact, if I delete the 0.5 from the scipy version, the output is the same as from the FieldTrip function.

I re-ran the Hilbert one (same parameters as above, just method='hilbert' instead of the default homomorphic:

EDIT: I wanted to add, I unfortunately do not have access to the referenced book by Oppenheim. I looked around a bit for publications (papers) on this, but couldn't find much that had the basic math. If anyone could point me to accessible literature that would be great!

[cbrnr]

I don't think the FieldTrip implementation is missing a log - it is only missing the 0.5

Of course, apparently my brain can't handle MATLAB code anymore 😄.

I have the book somewhere, I can take a look (next week).

Can you post the filtered signal for the three variants (i.e. original SciPy homomorphic, modified/corrected SciPy homomorphic, SciPy hilbert)?

[britta-wstnr]

apparently my brain can't handle MATLAB code anymore 😄.

Hehehe, wouldn't blame you for that @cbrnr

Here you go:

Maybe it'd be better to actually work with an example that has slow frequencies instead of DC offset - not so much to see here ...

[cbrnr]

I can't interpret these outputs. Why do the minimum phase filters just show horizontal lines? Initially, I thought that this was the problem, but they all look the same!

[britta-wstnr]

Yes, I think because there is a problem with all of them: they don't remove the DC offset. (The example was made to show exactly this behavior - which I think is unexpected/wrong.)

[britta-wstnr]

I had hoped the 0.5 constant would fix this - but while it makes the filter characteristics look better / closer to the zero phase filter, it does apparently not fix the offset issue.

[cbrnr]

I can't explain why there is such a large difference, but apparently the zero-phase filter seems to have much better stop band attenuation, even though this is not visible in the amplitude responses. Did you try doubling the filter order? Alternatively, if you have a high offset, would removing the offset before filtering be an option?

[larsoner]

Okay I think it has to do with the fact that scipy.signal.minimum_phase halves the filter order, which is in the SciPy docs, but should be optional. It makes the attenuation at DC two orders of magnitude worse. I'm working on a PR to SciPy for this, which we can fairly easily backport/vendor until our min version reaches the next SciPy release.

Testing script

import numpy as np
from scipy.io import savemat, loadmat
from scipy.signal import firwin, minimum_phase, freqz
import matplotlib.pyplot as plt

fs = 1000.
n_taps = 1001
cutoff = 10 / fs  # 10 Hz
h = firwin(n_taps, cutoff, width=cutoff, window="hann", pass_zero=False)  # 1001
savemat("test.mat", dict(h=h))
h_min = minimum_phase(h)  # 501
h_2 = np.convolve(h_min, h_min[::-1])  # restore to linear phase and length
assert len(h_2) == len(h)
h_min_matlab = loadmat("h_min.mat")["h_min"][0]  # using minphaserceps.m

fig, ax = plt.subplots(layout="constrained")
kwargs = dict(worN=10000, fs=fs)
w, H = freqz(h, **kwargs)
ax.plot(w, 20 * np.log10(np.abs(H)), label=f"H (μ={h.mean():0.1e}", zorder=5, lw=1)
w, H_min = freqz(h_min, **kwargs)
ax.plot(w, 20 * np.log10(np.abs(H_min)), label=f"H_min (μ={h_min.mean():0.1e})", lw=1)
w, H_2 = freqz(h_2, **kwargs)
ax.plot(w, 20 * np.log10(np.abs(H_2)), label=f"H_min^2 (μ={h_2.mean():0.1e})", lw=1)
w, H_min_matlab = freqz(h_min_matlab, **kwargs)
ax.plot(w, 20 * np.log10(np.abs(H_min_matlab)), label=f"H_min_matlab (μ={h_min_matlab.mean():0.1e})", zorder=3, lw=2)

ax.legend()
ax.set(xlim=(0, 20), ylim=(-100, 10), xlabel="Frequency", ylabel="Power (dB)")

[larsoner]

Fixed by scipy/scipy#19706, will backport to mne.fixes once it lands