Problem with scipy.signal.butter for designing filter bank I need to design a filter bank of band pass filters (code for implementation is shown below). There is a problem in that if my center frequencies fc are much smaller than the sampling frequency fs then you get the behavior as shown in the first graph:



A temporary solution is to change the sampling frequency to fs = fh*4 where 'fh' is the hi cutoff frequency of each respective band. This gives 



Which is what I want, however because the sampling frequency is no longer the true one (relative to the time series I then want to filter) I cant use it for filtering purposes using lfilter(b,a,x) where x is my time series because all of the bands will be identical since they have the same non-dimensional frequency! My question is what can I do to make this work as expected for a true fs? Presumably there is an issue with finite precision? 

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, lfilter
from scipy.signal import freqz



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


def fbank(fs):
    """
    Calculate third octave filter bank.
    An octave is the interval between two frequencies having
    a ratio of 2:1. Hence for 1/3 octave
    fh = fc*(sqrt(2)^(1/3))
    fl = fc/(sqrt(2)^(1/3))
    """

    # ISO standard center frequencies (bins)
    F0 = [50, 63, 80, 100, 125, 160, 200, 250, 315, 400, 500,
          630, 800, 1000, 1250, 1600, 2000, 2500, 3150, 4000,
          5000, 6300, 8000, 10000, 12500, 16000, 20000]


    plt.figure(figsize=(10,6))
    for fc in F0:
        fh = fc*(2**(1/6.0))
        fl = fc/(2**(1/6.0))
        fs = fh*4 # Temporary solution
        fs = 5e4
        b, a = butter_bandpass(fl,fh,fs)
        w, h = freqz(b, a, worN=2000)
        plt.loglog((fs * 0.5 / np.pi) * w, abs(h), label="fc = %d" % fc)

    plt.xlabel("Frequency (Hz)")
    plt.ylabel("Gain")
    plt.legend(loc='best')
    plt.savefig("test_filter.png", bbox_inches='tight')
    plt.show()
    return 0
Although Im not sure how I can display the frequency response using this method, I have found that it is workable to use second-order section method for filtering. So say we have a signal x with sampling frequency fs_real that needs to filtered using a constant percentage bandwidth (in this case third octaves) then the corresponding power spectrum will be given by the get_cpb function:

def butter_bandpass_sos(lowcut, highcut, fs, order=9):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    sos = butter(order, [low, high], btype='band', output='sos')
    return sos

def get_cpb(x,fs_real):
    """
    Calculate third octave spectrum from pressure time series.
    An octave is the interval between two frequencies having
    a ratio of 2:1. Hence for 1/3 octave
    fh = fc*(sqrt(2)^(1/3))
    fl = fc/(sqrt(2)^(1/3))
    """

    cpb = []

    for fc in F0:
        fh = fc*(2**(1/6.0))
        fl = fc/(2**(1/6.0))
        sos = butter_bandpass_sos(fl,fh,fs_real)
        y = sosfilt(sos,x)
        cpb.append( 20*np.log10(np.std(y)/2e-5) )


    return F0, cpb