2
$\begingroup$

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()

bandpassfilter-1

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

fs = 100000.0

bandpassfilter-2

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))
|improve this question
$\endgroup$
  • 3
    $\begingroup$ this is happening using double precision? doesn't scipy.signal.freqz() use double precision? because, when increasing the sample rate to a high value turns a nice frequency response into a crappy frequency response, that usually is because of the cosine problem. this is because $$ \cos\big(\tfrac{\omega}{f_\mathrm{s}}\big) \approx 1$$ is so close to 1 that much of the precision regarding $\tfrac{\omega}{f_\mathrm{s}}$ is lost. $\endgroup$ – robert bristow-johnson Jan 22 at 19:39
  • $\begingroup$ i have an answer at another question about plotting frequency response "manually" (instead of using freqz()) that deals with this cosine problem. $\endgroup$ – robert bristow-johnson Jan 22 at 19:41
  • 1
    $\begingroup$ You should try to split your filter into an order-2 filter and an order-1 filter and then cascade them. You should get a better conditioned filter. You simply need to add the frequency response of both filters. $\endgroup$ – Ben Jan 22 at 21:33
  • 2
    $\begingroup$ @Ben An easier way is to use output='sos' and sosfreqz() and sosfilt(), which does that automatically. It should be the first thing you reach for when filtering. $\endgroup$ – endolith Jan 22 at 22:30
  • 1
    $\begingroup$ thank you all, i was able to get it using sosfreqz() and sosfilt(). I will update the OP $\endgroup$ – khuynh Jan 22 at 22:42
3
$\begingroup$

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:

pole-zero plot near +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.

|improve this answer
$\endgroup$
2
$\begingroup$

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.

|improve this answer
$\endgroup$
  • $\begingroup$ Is there any any benefit to downsampling, applying the filter, and then resampling compared to using cascaded second-order sections and applying that to the original signal? I assume computation is lower with a downsampled signal. $\endgroup$ – khuynh Jan 22 at 22:44
  • 1
    $\begingroup$ The question about possible benefits appears to be suitable for its own separate question here on the dsp stackexchange. $\endgroup$ – hotpaw2 Jan 23 at 0:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.