如何用scipy.signal.butter实现多带通滤波器？解决滤波信号趋近零问题

问题描述 #

基于巴特沃斯带通滤波器实现方案编写多频段滤波器代码后，滤波输出信号幅值接近零，无法正常绘制频谱。咨询是否需要对各频段滤波器系数做归一化，以及修复该滤波器的方法。

原实现代码 #

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 = []
    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]

    # 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
    f7, f8 = 1000, 1000.1
    cut_off = [(f1, f2), (f3, f4), (f5, f6), (f7, f8)]

    t_toy_signal, toy_signal = get_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)

    # Plot the frequency response
    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()

    # 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
    # toy_signal_filtered = sosfilt(sos, toy_signal)
    toy_signal_filtered = np.sum([sosfilt(sos, toy_signal) for sos in sos_list], axis=0)

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

    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='upper right')

    plt.tight_layout()
    plt.show()

    plt.figure()
    plt.plot(np.arange(N) / N * fs, 20 * np.log10(abs(X)), 'r-', label='FFT original signal')
    plt.plot(np.arange(N) / N * fs, 20 * np.log10(abs(Y)), 'g-', label='FFT filtered signal')
    plt.xlim(xmax=fs / 2)
    plt.ylim(ymin=-20)
    plt.ylabel(r'Power Spectrum (dB)', fontsize=fontsize)
    plt.xlabel("frequency (Hz)", fontsize=fontsize)
    plt.grid()
    plt.legend(loc='upper right')
    plt.tight_layout()
    plt.show()

问题原因 #

原代码中把所有频段的sos滤波器系数堆叠后用sosfilt处理，相当于让信号依次通过每个带通滤波器，最终输出只会保留所有滤波器共同允许通过的频段（即各频段的交集）。由于目标频段互不重叠，最终输出自然接近零。

另外，不需要对sos系数做额外归一化——scipy.signal.butter生成的巴特沃斯滤波器已经是归一化处理的，通带内增益为1。

修复方案 #

正确的多频段滤波逻辑是：分别对每个频段做带通滤波，再将各频段的滤波结果相加，以此保留所有目标频段的信号成分。

修改代码中滤波部分的实现：

# 替换原串联滤波的代码
toy_signal_filtered = np.sum([sosfilt(sos, toy_signal) for sos in sos_list], axis=0)

结果说明 #

• 采用求和方式后，滤波信号的频谱会完整保留所有目标频段的成分，幅值恢复正常，频谱图显示效果符合预期
• 若误用np.mean替代np.sum，信号幅值会被平均稀释，虽然能看到目标频段，但整体幅值偏低
• 单窄带测试时，直接用单个带通滤波器处理即可得到正常结果

内容的提问来源于stack exchange，提问作者Thoth