我需要设计一个带通滤波器的滤波器组（实现代码如下所示）。存在一个问题，如果我的中心频率fc远小于采样频率fs，那么您会得到第一张图中所示的行为：

一个临时解决方案是将采样频率更改为fs = fh*4“fh”是每个相应频段的高截止频率。这给

这就是我想要的，但是因为采样频率不再是真实的（相对于我想要过滤的时间序列）我不能使用它来过滤使用lfilter(b,a,x)where xis my time series 因为所有波段都是相同的因为它们具有相同的无量纲频率！我的问题是我该怎么做才能使这项工作按预期进行fs？大概存在有限精度的问题？

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