Bandpass butterworth filter frequencies in scipy

Question

I'm designing a bandpass filter in scipy following the cookbook. However, if I decrecrease the filtering frequencies too much I end up with garbage at high order filters. What am I doing wrong?

from scipy.signal import butter, lfilter

def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a

if __name__ == "__main__":
    import numpy as np
    import matplotlib.pyplot as plt
    from scipy.signal import freqz  
    # Sample rate and desired cutoff frequencies (in Hz).
    fs = 25
    # Plot the frequency response for a few different orders.
    plt.figure(1)
    plt.clf()
    for order in [1, 3, 5, 6, 9]:
        b, a = butter_bandpass(0.5, 4, fs, order=order)
        w, h = freqz(b, a, worN=2000)#np.logspace(-4, 3, 2000))
        plt.semilogx((fs * 0.5 / np.pi) * w, abs(h), label="order = %d" % order)  
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Gain')
    plt.grid(True)
    plt.legend(loc='best')

    plt.figure(2)
    plt.clf()
    for order in [1, 3, 5, 6, 9]:
        b, a = butter_bandpass(0.05, 0.4, fs, order=order)
        w, h = freqz(b, a, worN=2000)#np.logspace(-4, 3, 2000))
        plt.semilogx((fs * 0.5 / np.pi) * w, abs(h), label="order = %d" % order)  
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Gain')
    plt.grid(True)
    plt.legend(loc='best')

    plt.show()

Update: the issue was discussed and apparently solved on Scipy 0.14. However, the plot still looks really bad after Scipy update. What's wrong?

Answer 1

1. Don't use b, a = butter for high-order filters, whether in Matlab or SciPy or Octave. Transfer function format has numerical stability problems, because some of the coefficients are very large while others are very small. This is why we changed the filter design functions to use zpk format internally. To see the benefits of this, you need to use z, p, k = butter(output='zpk') and then work with poles and zeros instead of numerator and denominator.
2. Don't do high-order digital filters in a single stage. This is a bad idea no matter what software or hardware you're implementing them on. Typically it's best to break them up into second-order sections. In Matlab, you can use zp2sos to generate these automatically. In SciPy, you can use sos = butter(output='sos') and then filter using sosfilt() or sosfiltfilt(). This is the recommended way to filter for most applications.