如何使用 scipy.signal.butter 实现多带通滤波器

您的滤波信号会相互“乘法”，因为您是按系列计算的。在我下面的代码中，我更改了您的代码，以便它逐个计算您的filtered_signals，然后将它们一一绘制在您的图表上。我希望这就是你试图实现的目标：

from scipy.signal import butter, sosfreqz, sosfilt
from scipy.signal import spectrogram
import matplotlib
import matplotlib.pyplot as plt
from scipy.fft import fft
import numpy as np
import scipy


def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    print(high)
    sos = butter(order, [low, high], analog=False, btype='band', output='sos')
    return sos


def multiband_filter(data, bands, fs, order=10):
    sos_list = []
    for lowcut, highcut in bands:
        sos = butter_bandpass(lowcut, highcut, fs, order=order)
        scalar = max(abs(fft(sos, 2000)))
        sos = sos / scalar
        sos_list += [sos]

    # sos_list = [butter_bandpass(lowcut, highcut, fs, order=order) for lowcut, highcut in bands]

    # Combine filters into a single filter
    sos = np.vstack(sos_list)

    # Apply the multiband filter to the data
    y = sosfilt(sos, data)

    return y, sos_list


def get_toy_signal():
    t = np.arange(0, 0.3, 1 / fs)

    fq = [-np.inf] + [x / 12 for x in range(-9, 3, 1)]

    mel = [5, 3, 1, 3, 5, 5, 5, 0, 3, 3, 3, 0, 5, 8, 8, 0, 5, 3, 1, 3, 5, 5, 5, 5, 3, 3, 5, 3, 1]
    acc = [5, 0, 8, 0, 5, 0, 5, 5, 3, 0, 3, 3, 5, 0, 8, 8, 5, 0, 8, 0, 5, 5, 5, 0, 3, 3, 5, 0, 1]

    toy_signal = np.array([])

    for kj in range(len(mel)):
        note_signal = np.sum([np.sin(2 * np.pi * 440 * 2 ** ff * t)
                              for ff in [fq[acc[kj]] - 1, fq[acc[kj]], fq[mel[kj]] + 1]], axis=0)

        zeros = np.zeros(int(0.01 * fs))
        toy_signal = np.concatenate((toy_signal, note_signal, zeros))

    toy_signal += np.random.normal(0, 1, len(toy_signal))

    toy_signal = toy_signal / (np.max(np.abs(toy_signal)) + 0.1)
    t_toy_signal = np.arange(len(toy_signal)) / fs

    return t_toy_signal, toy_signal


if __name__ == "__main__":

    fontsize = 12
    # Sample rate and desired cut_off frequencies (in Hz).
    fs = 3000

    f1, f2 = 100, 200
    f3, f4 = 470, 750
    f5, f6 = 800, 850
    cut_off = [(f1, f2), (f3, f4), (f5, f6)]
    #cut_off = [(f1, f2), (f3, f4)]
    #cut_off = [(f1, f2)]
    # cut_off = [f1]

    t_toy_signal, toy_signal = get_toy_signal()
    # toy_signal -= np.mean(toy_signal)
    # t_toy_signal = wiener(t_toy_signal)

    fig, ax = plt.subplots(6, 1, figsize=(8, 12))
    fig.tight_layout()

    ax[0].plot(t_toy_signal, toy_signal)
    ax[0].set_title('Original toy_signal', fontsize=fontsize)
    ax[0].set_xlabel('Time (s)', fontsize=fontsize)
    ax[0].set_ylabel('Magnitude', fontsize=fontsize)
    ax[0].set_xlim(left=0, right=max(t_toy_signal))
    sos_lists=[]
    for cut in cut_off:
        print(cut)
        sos_list = [butter_bandpass(lowcut, highcut, fs, order=10) for lowcut, highcut in [cut]]
        sos_lists.append(np.vstack(sos_list))
    # Combine filters into a single filter
    sos = np.vstack(sos_list)



    # w *= 0.5 * fs / np.pi  # Convert w to Hz.
    #####################################################################
    # First plot the desired ideal response as a green(ish) rectangle.
    #####################################################################

    # Plot the frequency response
    for i in range(len(cut_off)):
        w, h = sosfreqz(sos_lists[i], worN=2000)
        ax[1].plot(0.5 * fs * w / np.pi, np.abs(h), label=f'Band {i + 1}: {cut_off[i]} Hz')

    ax[1].set_title('Multiband Filter Frequency Response')
    ax[1].set_xlabel('Frequency [Hz]')
    ax[1].set_ylabel('Gain')
    ax[1].legend()
    # ax[1].set_xlim(0, max(*cut_off) + 100)

    #####################################################################
    # Spectrogram of original signal
    #####################################################################

    f, t, Sxx = spectrogram(toy_signal, fs,
                            nperseg=930, noverlap=0)
    ax[2].pcolormesh(t, f, np.abs(Sxx),
                     norm=matplotlib.colors.LogNorm(vmin=np.min(Sxx), vmax=np.max(Sxx)),
                     )

    ax[2].set_title('Spectrogram of original toy_signal', fontsize=fontsize)
    ax[2].set_xlabel('Time (s)', fontsize=fontsize)
    ax[2].set_ylabel('Frequency (Hz)', fontsize=fontsize)

    #####################################################################
    # Compute filtered signal
    #####################################################################

    # Apply the multiband filter to the data
    toy_signal_filtereds=[]
    for sos in sos_lists:
        toy_signal_filtered = sosfilt(sos, toy_signal, )
        toy_signal_filtereds.append([toy_signal_filtered])


    #####################################################################
    # Spectrogram of filtered signal
    #####################################################################

    f, t, Sxx = spectrogram(toy_signal_filtered, fs,
                            nperseg=930, noverlap=0)

    ax[3].pcolormesh(t, f, np.abs(Sxx),
                     norm=matplotlib.colors.LogNorm(vmin=np.min(Sxx),
                                                    vmax=np.max(Sxx))
                     )

    ax[3].set_title('Spectrogram of filtered toy_signal', fontsize=fontsize)
    ax[3].set_xlabel('Time (s)', fontsize=fontsize)
    ax[3].set_ylabel('Frequency (Hz)', fontsize=fontsize)
    for toy_signal_filtered in toy_signal_filtereds:
        ax[4].plot(t_toy_signal, toy_signal_filtered[0])
    ax[4].set_title('Filtered toy_signal', fontsize=fontsize)
    ax[4].set_xlim(left=0, right=max(t_toy_signal))
    ax[4].set_xlabel('Time (s)', fontsize=fontsize)
    ax[4].set_ylabel('Magnitude', fontsize=fontsize)

    N = 1512
    X = fft(toy_signal, n=N)
    Ys=[]
    for toy_signal_filtered in toy_signal_filtereds:
        Y = fft(toy_signal_filtered[0], n=N)
        Ys.append(Y)

    # fig.set_size_inches((10, 4))
    ax[5].plot(np.arange(N) / N * fs, 20 * np.log10(abs(X)), 'r-', label='FFT original signal')
    for bigY in Ys:
        ax[5].plot(np.arange(N) / N * fs, 20 * np.log10(abs(bigY)), 'g-', label='FFT filtered signal')
    ax[5].set_xlim(xmax=fs / 2)
    ax[5].set_ylim(ymin=-20)
    ax[5].set_ylabel(r'Power Spectrum (dB)', fontsize=fontsize)
    ax[5].set_xlabel("frequency (Hz)", fontsize=fontsize)
    ax[5].grid()
    ax[5].legend(loc=0)

    plt.tight_layout()
    plt.show()

对于频率响应，在您的情况下，您可能不应该使用对数标度图。 但是，在 y 轴上绘制了带有对数缩放的图。ax[1].semilogy

更简单和推荐的方法是沃伦在评论中写道的。只需计算单独带通滤波信号的总和即可。

话虽如此，对于想要创建和应用单个多波段滤波器的人来说，他可以尝试通过组合滤波器来实现这一点：

• 低通（切断最后一个通滤波器以上的所有内容），
• 高通（将所有低于第一个通滤波器的东西都切断），
• N-1 带阻滤波器，其中 N 是通带数（用于切割滤波器之间的部件）。

不过，要让它稳定（注意过滤顺序）可能很困难，而且要让它变得陡峭更难。

觉得很有趣，自己试了一下：

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


def multi_band_filter(bands, subfilter_order=5):
    # high-pass filter
    nyq = 0.5 * fs
    normal_cutoff = bands[0][0] / nyq
    b, a = butter(subfilter_order, normal_cutoff, btype='highpass', analog=False)
    all_b = [b]
    all_a = [a]

    # band-stop filters
    for idx in range(len(bands) - 1):
        normal_cutoff1 = bands[idx][1] / nyq
        normal_cutoff2 = bands[idx+1][0] / nyq
        b, a = butter(subfilter_order, [normal_cutoff1, normal_cutoff2], btype='bandstop', analog=False)
        all_b.append(b)
        all_a.append(a)

    # low-pass filter
    normal_cutoff = bands[-1][1] / nyq
    b, a = butter(subfilter_order, normal_cutoff, btype='lowpass', analog=False)
    all_b.append(b)
    all_a.append(a)

    # combine filters:
    combined_a = all_a[0]
    for a in all_a[1:]:
        combined_a = np.convolve(a, combined_a)
    combined_b = all_b[0]
    for b in all_b[1:]:
        combined_b = np.convolve(b, combined_b)

    return combined_b, combined_a


bands = [[400, 700], [1000, 1500]]

fs = 8000
time = np.arange(0, 1 - 0.5/fs, 1/fs)
signal_to_filter = np.sum([np.sin(2 * np.pi * (freq + 0.01 * np.random.random()) * time + np.pi*np.random.random()) for freq in range(10, 3800)], axis=0)

b, a = multi_band_filter(bands)
filtered_signal = lfilter(b, a, signal_to_filter)
original_spectrum = fft(signal_to_filter)
filtered_signal_spectrum = fft(filtered_signal)

plt.figure(figsize=(16, 10))
plt.plot(np.linspace(0, fs, len(original_spectrum)), np.abs(original_spectrum), color='b')
plt.plot(np.linspace(0, fs, len(filtered_signal_spectrum)), np.abs(filtered_signal_spectrum), color='orange')
plt.xlim([0, 4000])
plt.show()

SOS 版本

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


def multi_band_filter(bands, subfilter_order=5):
    # high-pass filter
    nyq = 0.5 * fs
    normal_cutoff = bands[0][0] / nyq
    sos = butter(subfilter_order, normal_cutoff, btype='highpass', analog=False, output='sos')
    all_sos = [sos]

    # band-stop filters
    for idx in range(len(bands) - 1):
        normal_cutoff1 = bands[idx][1] / nyq
        normal_cutoff2 = bands[idx+1][0] / nyq
        sos = butter(subfilter_order, [normal_cutoff1, normal_cutoff2], btype='bandstop', analog=False, output='sos')
        all_sos.append(sos)

    # low-pass filter
    normal_cutoff = bands[-1][1] / nyq
    sos = butter(subfilter_order, normal_cutoff, btype='lowpass', analog=False, output='sos')
    all_sos.append(sos)

    # combine filters:
    combined_sos = np.vstack(all_sos)
    return combined_sos


bands = [[400, 700], [1000, 1500]]

fs = 8000
time = np.arange(0, 1 - 0.5/fs, 1/fs)
signal_to_filter = np.sum([np.sin(2 * np.pi * (freq + 0.01 * np.random.random()) * time + np.pi*np.random.random()) for freq in range(10, 3800)], axis=0)

sos = multi_band_filter(bands)
filtered_signal = sosfilt(sos, signal_to_filter)
original_spectrum = fft(signal_to_filter)
filtered_signal_spectrum = fft(filtered_signal)

plt.figure(figsize=(16, 10))
plt.plot(np.linspace(0, fs, len(original_spectrum)), np.abs(original_spectrum), color='b')
plt.plot(np.linspace(0, fs, len(filtered_signal_spectrum)), np.abs(filtered_signal_spectrum), color='orange')
plt.xlim([0, 4000])
plt.show()

w, h = freqz(b, a)
freq_domain = np.linspace(0, fs/2, len(w))
plt.figure(figsize=(16, 10))
plt.plot(freq_domain, 20 * np.log10(abs(h)), 'b')
plt.show()

如您所见，过滤器的斜率不是很陡。