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

Question

基于带通滤波器here，我尝试使用下面的代码制作一个多带滤波器。然而，滤波后的信号接近于零，这会影响绘制频谱时的结果。每个频段的滤波器系数是否需要归一化？请有人建议我如何修复过滤器？

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


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


def multiband_filter(data, bands, fs, order=10):
    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_list = [butter_bandpass(lowcut, highcut, fs, order=10) for lowcut, highcut in cut_off]

    # 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_list[i], worN=2000)
        ax[1].semilogy(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_filtered = sosfilt(sos, toy_signal)

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

    ax[4].plot(t_toy_signal, toy_signal_filtered)
    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)
    Y = fft(toy_signal_filtered, n=N)

    # fig.set_size_inches((10, 4))
    ax[5].plot(np.arange(N) / N * fs, 20 * np.log10(abs(X)), 'r-', label='FFT original signal')
    ax[5].plot(np.arange(N) / N * fs, 20 * np.log10(abs(Y)), '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()

Answer 1

对于您的情况的频率响应，您可能不应该使用对数标度图。然而 ax.semilogy 在 y 轴上绘制了对数缩放图。

也许尝试下面的代码可以使图表更直观地阅读：

for i in range(len(cut_off)):
    w, h = sosfreqz(sos_list[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_ylim(-0.1,1.5)