Why does a higher sampling frequency mess up my bandpass filter?

Question

I was designing a bandpass filter in python using some of the scipy.signal modules.

I am plotting the frequency response of my filter to verify that my desired frequency is in the passband. However, when I make the sampling frequency larger, the frequency response of my filter gets completely messed up.

I currently have a signal that I sample >300KSPS and am trying to create a bandpass filter for some fairly low frequencies (1-100hz). Could somebody explain why this happens?

For example, the following code yields this frequency response:

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

N = 3
fs = 10000.0
low = 100.0
high = 150.0

nyq = fs * 0.5

Wn = [low/nyq, high/nyq]

b, a = scipy.signal.butter(N, Wn, btype='bandpass')
w, h = scipy.signal.freqz(b, a, worN=round(fs/2))

ax = plt.subplot(121)
ax.set(title='filter frequency response',
       xlabel='frequency [hz]',
       ylabel='gain',
       xlim=(low/2, high*2))

ax.plot((nyq / np.pi) * w, abs(h), label='filter freq response')
ax.axvline(x=125.0, linestyle='--', alpha=0.5, c='black', label='f=125.0')
ax.axhline(y=np.sqrt(0.5), linestyle='--', alpha=0.5, c='black', label='sqrt(0.5)')
ax.grid()
ax.legend()

Then when I change my sampling frequency to 100000.0, the response turns out to be this:

fs = 100000.0

EDIT: Outputting the filter as second-order sections and using scipy.signal.sosfreqz yielded the proper filter. See code below for the modified lines:

sos = scipy.signal.butter(N, Wn, btype='bandpass', output='sos')
w, h = scipy.signal.sosfreqz(sos, worN=round(fs/2))

Answer 1

As others said in the comments, this looks like numerical error. 3rd-order filters are not typically prone to this, but the higher your sampling frequency, the closer the poles move to +1:

You might benefit from splitting this into second-order sections. The easiest way is to use output='sos' and sosfreqz() and sosfilt(), which handle the splitting automatically.

These second-order section functions should be the first thing you try, actually. They didn't exist in scipy until recently, so they aren't the default, for backwards-compatibility, but should always work better than output='ba', lfilter, etc.

Answer 2

With an IIR filter, the higher the ratio is of sampling frequency to bandpass filter frequency, the closer the poles are to the complex unit circle in the Z-transform. Poles close to the unit circle can produce numerical instabilities (ill conditioned computations) when computing the filters response. This is because any numerical quantization or "rounding errors" can get amplified in the arithmetic recurrences if there is only a tiny difference in location between the pole and the unit circle, even if the pole is "supposed" to be inside.

A better solution is to first downsample the signal before filtering to use an IIR filter that is possibly of a lower order and more stable. You can upsample back up, if needed, after the filter has been applied.