5
$\begingroup$

I am new to signal processing and I want to implement an LPF using SciPy. In order to do so, I used the following python code from: here

The code itself looks like:

import numpy as np from scipy.signal import butter, lfilter, freqz import matplotlib.pyplot as plt def butter_lowpass(cutoff, fs, order=5): nyq = 0.5 * fs normal_cutoff = cutoff / nyq b, a = butter(order, normal_cutoff, btype='low', analog=False) return b, a def butter_lowpass_filter(data, cutoff, fs, order=5): b, a = butter_lowpass(cutoff, fs, order=order) y = lfilter(b, a, data) return y # Setting standard filter requirements. order = 6 fs = 30.0 cutoff = 3.667 b, a = butter_lowpass(cutoff, fs, order) # Plotting the frequency response. w, h = freqz(b, a, worN=8000) plt.subplot(2, 1, 1) plt.plot(0.5*fs*w/np.pi, np.abs(h), 'b') plt.plot(cutoff, 0.5*np.sqrt(2), 'ko') plt.axvline(cutoff, color='k') plt.xlim(0, 0.5*fs) plt.title("Lowpass Filter Frequency Response") plt.xlabel('Frequency [Hz]') plt.grid() # Creating the data for filteration T = 5.0 # value taken in seconds n = int(T * fs) # indicates total samples t = np.linspace(0, T, n, endpoint=False) data = np.sin(1.2*2*np.pi*t) + 1.5*np.cos(9*2*np.pi*t) + 0.5*np.sin(12.0*2*np.pi*t) # Filtering and plotting y = butter_lowpass_filter(data, cutoff, fs, order) plt.subplot(2, 1, 2) plt.plot(t, data, 'b-', label='data') plt.plot(t, y, 'g-', linewidth=2, label='filtered data') plt.xlabel('Time [sec]') plt.grid() plt.legend() plt.subplots_adjust(hspace=0.35) plt.show()

As it can be seen in the figure, a 6th order LPF with fs = 30 Hz and a cutoff frequency of 3.667 Hz filters the signal nicely:

enter image description here

In my application, the signal has a sampling frequency of 1.5 MSPS (1.5 MHz, the ADC has to sample with such a high rate, I cannot change that) and I want to use an LPF with 100 Hz cutoff frequency to filter my incoming signal.

When I change the corresponding two parameters in this test code just to see what would happen with the test data:

fs = 1.5e6 cutoff = 100

then the result becomes completely uninterpretable: enter image description here

Testing it by an empirical way, I could find out that the problem is somewhat related to the cutoff/fs ratio. When it gets too small, strange things happen. However, I do not understand why this is happening. Also, any advice related to how to implement a 100 Hz cutoff LPF on a signal with a high sampling frequency would be appreciated.

Thank you very much in advance!

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ fc/fs ratio may be part of the problem but also order of the filter may play part in it. Try by splitting the filter to second order sections (SOS). $\endgroup$ Commented Apr 4, 2023 at 12:28
  • 2
    $\begingroup$ There is something called the "cosine problem" which makes it hard to compute coefficients with single-precision floating point. $$ \cos(\omega_0) \ = \ 1 - 2 \sin^2 \left( \frac{\omega_0}{2} \right) $$ where $$ \omega_0 \triangleq 2 \pi \frac{f_c}{f_s} $$ All of the data that defines the cutoff frequency are contained in the bits that fall offa the edge when $\omega_0$ gets really small. $\endgroup$ Commented Apr 5, 2023 at 2:08
  • 1
    $\begingroup$ I have the "cosine problem" when I wake up in the morning before having my coffee. It may help with Robert's good comment that the poles for the 2nd order biquad (SOS) filter are given by the polynomial $1 - 2rcos(\omega_o)z^{-1}+r^2z^{-2}$ where $r$ is the radius to the pole (hiqh Q means near the unit circle). Thus the 2nd coefficient defines the resonant frequency (pole location) in the filter. $\endgroup$ Commented Apr 5, 2023 at 3:09

2 Answers 2

Reset to default
8
$\begingroup$

Also, any advice related to how to implement a 100 Hz cutoff LPF on a signal with a high sampling frequency would be appreciated.

Easy: use second order representation (output = 'sos') and sosfilt() instead.

However, I do not understand why this is happening.

The problem is using transfer function representation ([b,a]) for IIR filters which is generally a bad idea. The filter is defined by its poles and zeros and if the cutoff frequency is very low (as compared to the sample rate), the poles become very close to the unit circle, which is the limit for stability.

Converting from poles & zeros to transfer function and back is a numerically tricky and at some point even double precision floating point isn't good enough. Specifically finding the roots of a polynomial where the roots are very close together is difficult.

You can see some of the problem simply by inspecting your filter coefficients. The numerator coefficients b are in the order of $10^{-16}$ which is a hint that something unhealthy is happening here. A more formal test is to convert [b,a] to poles and zeros. You'll find poles outside the unit circle (due to numerical noise) which causes the instability.

Improve this answer
$\endgroup$
0
7
$\begingroup$

In addition to Hilmar's good comments in his own answer, see this question as well as Rick Lyon's interesting blog post which shows similar pole pattern graphs, demonstrating why very low cut-off filters are a challenge even in the second-order section structure, when implemented in fixed point. I copied the relevant plot here showing the pattern of pole locations for the classic SOS using any of the Direct Form standard structures:

poles

Note that a filter with a low cutoff relative to the sampling rate will require a pole close to $z=1$ on the z plane inside the unit circle. This is the area in the plot above where the possible pole locations once quantized becomes very sparse. We see that the possible error only increases as the ratio of sampling frequency to cut-off frequency increases- so increasing the sampling rate only makes things worst. That said this is resolved by using even more precision (number of bits) within the filter structure to meet filter performance requirements than we might otherwise expect with typical FIR structures.

The other way around this (assuming the resources used for the number of bits required becomes more of an issue than using more dedicated multipliers) other structures to be considered are the Normal (or Coupled) Form by Rader and Gold, or alternate structures proposed by Agarwal and Burrus.

For low pass filters with cutoffs significantly smaller than the sampling rate, it is often just as efficient to implement the filter as a decimating FIR filter, which also has the benefit of avoiding all the possible issues with IIR filters (instability, quantization noise growth, limit cycles, overflow oscillations, etc). I have demonstrated such an efficiency comparison as well as additional considerations here at this post where a 40 Hz low pass filter running at a 48KHz clock rate was implemented with 8.45 real multiplications per clock cycle (compared to the 15 real multiplications per clock cycle require for the 6th order IIR needed to get the same magnitude response). This comes at the cost of memory storage. As mentioned, the time delay in the demonstrated implementation would be easily resolved for the case of a low pass filter by decimating a minimum phase FIR filter instead of linear phase for the implementation.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks for the detailed explanation! I would upvote your answer, but don't have enough reputation for it... $\endgroup$ Commented Apr 4, 2023 at 13:27
  • $\begingroup$ @gvg No problem! I am glad it was helpful $\endgroup$ Commented Apr 4, 2023 at 15:02

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.