Why is implementing a digital LPF with low cutoff frequency but high sampling frequency infeasible?

Question

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:

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:

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!

Answer 1

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.

Answer 2

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:

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.