Scipy butter filter - remove DC offset w/ high sample rate

Question

I am having a hard time figuring out how to employ a high pass filter to remove the DC offset of my data signal with the "scipy butter" function because my sample rate is quite high. The crux of my issue seems to be that my sample rate Fs is very high and my cut off frequency is very low (ideally just the DC offset component, but could be 1 to 10Hz with a sample rate of 125ksps).

Here is my example code:

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt


# inital parameters
Fs = 125000         # sample rate (if changed to 125sps looks better)
n = 9               # filter order
Fc = 10             # cut off freq (Hz)
Nyq = Fs/2.         # Nyquist freq (1/2 Fs)
Fcc = float(Fc)/float(Nyq)  # normalized cut off freq


# find butterworth xfer functions
sos = signal.butter(n, Fcc, btype='hp', output='sos', analog=False)
# find freq/mag
w, h = signal.sosfreqz(sos, worN=1500, fs=Fs)
#convert to log scale
h_db = 20*np.log10(np.abs(h))

# plot filter response
fig1, ax1 = plt.subplots(1, 1)
ax1.plot(0.5*Fs*(w/np.pi), h_db)
ax1.set_title('db log')
plt.show()

I think this is working like it is supposed to because dropping the sample rate down to 125 instead of 125,000 yields a response that looks reasonable. I think I'm missing something fundamental here with the Fc vs Fs.

Here is the response plot:

Here is a screenshot of the data I'm trying to apply this filter to:

As you can see there is a DC offset I'd like to filter out of my data file with a high pass filter. Right now I'm just trying to get the filter response to behave so I'm not worried about processing the data yet - I just added the data file screenshot for reference. I don't need anything fancy. So my questions are the following:

1. Why do I have ringing and such an aggressive response?
2. How does the low cut off frequency and high sample rate affect the filter response?
3. Why did Scipy go from the b,a approach in the butter() filter to the "sos"? I"m curious why they switch to this new sos appraoch. My understanding is that it is more accurate and/or prevents errors in the b,a approach. I'm unsure if this is just bolting second order transfer functions together though cascading simpler filters or what.
4. Also, you can notice (if you print(w) to see my freq axis) that my frequencies seem to run from 0 to 2*pi. Is that normal or should it be from 0 to pi?
5. If question #4 is is "yes this should be 0 to 2*pi" than should line #23 in script:

actually be

ax1.plot(0.5*Fs*(w/(2*np.pi)), h_db)

so that my frequency range in the final result/plot runs from 0 to 1/2*Fs?

Thanks for any help!

Answer 1

Why do I have ringing and such an aggressive response?

Because you designed a really aggressive filter, which is borderline unstable. All the poles are less than 1/1000 away from the unit circle.

How does the low cut off frequency and high sample rate affect the filter response?

The lower the cutoff, the higher the order, and the higher the sample rate, the more it will ring (in terms of samples).

Why did Scipy go from the b,a approach in the butter() filter to the "sos"?

Because in transfer function representation the filter would actually BE unstable, i.e. the numerical precision of 64-bit floating point would NOT be enough to run this filter as a high order polynomial.

I"m curious why they switch to this new sos appraoch. My understanding is that it is more accurate and/or prevents errors in the b,a approach.

NO, NO, NO. Transfer function representation is the worst and should be avoided in favor of poles, zeros or second order section. Once you are in transfer function, you often can't get back to poles and zeros. It requires finding the roots of a high order polynomial which is a numerically ill defined problem.

Also, you can notice (if you print(w) to see my freq axis) that my frequencies seem to run from 0 to 2*pi. Is that normal or should it be from 0 to pi?

If you use a function, it's always a good idea to read the documentation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.sosfreqz.html#scipy.signal.sosfreqz

Since you pass in the sample rate as an argument, the function already returns the frequency axis. By your attempt to normalize it again, you get complete non-sense. The dead give away here is the scale factor $10^9$ on your frequency axis. Your current axis tops out at 1.25 GHz.

If question #4 is is "yes this should be 0 to 2*pi" than should line #23 in script:

It should simply be ax1.plot(w, h_db)

The question you didn't ask:

What I should be doing ?

What @Ben said: use a first order filter AND either subtract the mean from the signal or seed the filter state with the signal mean.

Answer 2

1. The order of your filter is 9, which is high. For DC-removal purposes, an order-1 IIR filter is often used. Order-1 IIR filters don't ring. You could tune the $\alpha$ value to suit your need. $H(z) = \frac{1-z^{-1}}{1-\alpha z^{-1}}$ You could also perform a DC-block removal. In this case, simply compute the average of your signal and subtract it from your signal. https://www.embedded.com/dsp-tricks-dc-removal/

2. A low cut-off frequency usually means a longer transient. Especially since the order of your filter is 9. A higher sampling frequency means a longer transient, i.e. it will take more samples for your filter to settle. However, the settling time in seconds should be approximately the same even if you vary the sampling frequency.

3. SOS means second-order sections. The filter is converted into a cascade of order-2 IIR filter. These filters are more numerically stable. https://www.dsprelated.com/freebooks/filters/Series_Second_Order_Sections.html

4. and 5 I'm not an expert in Scipy.