Scipy数字滤波器通带内信号衰减问题求助

我使用高通滤波器处理信号时，发现滤波后的信号在通带内出现明显衰减，无法理解原因，特寻求解释。

输出图显示滤波后通带内的信号功率有明显下降，放大图更清晰展示了这一现象。

完整Python可运行代码如下：

# %% modules
import numpy as np
import scipy as sc
import matplotlib.pyplot as plt

# %% load data
data = np.load("data.npz")
t = data['t'].squeeze()     # (s)       time vector
dt = np.diff(t)[0]          # (s)       sample period
fs = 1/dt                   # (Hz)      sampling rate
y = data['y'].squeeze()     # (null)    signal

# %% create filter and filter data
# define filter function
def highpass(data: np.ndarray, cutoff: float, sample_rate: float, poles: int = 5):
    sos = sc.signal.butter(poles, cutoff, 'highpass', fs=sample_rate, output='sos')
    filtered_data = sc.signal.sosfiltfilt(sos, data)
    return filtered_data, sos

# filter the data
cutoff = 0.1                                    # (Hz)      cutoff frequency
poles  = 100                                    # (null)    filter order
yF, sosfilt = highpass(y, cutoff, fs, poles)

# %% compute power spectrum
fout1,psd  = sc.signal.welch(y.squeeze(),fs)
fout2,psdF = sc.signal.welch(yF.squeeze(),fs)

# %% plot
fig,ax = plt.subplots(3,1)
fig.tight_layout()

# time domain
axh = ax[0]
axh.plot(t, y)
axh.set(xlabel='Time (s)',
          ylabel='Amplitude')
axh.grid()

# filter characteristics
axh = ax[1]
w, h = sc.signal.sosfreqz(sosfilt, worN=8000,fs=fs)
db = 20 * np.log10(abs(h))
yUB = np.max(db[np.isfinite(db)])
yLB = np.min(db[np.isfinite(db)])
axh.plot(w, db)
axh.fill_betweenx([yLB,yUB],cutoff,fout1[-1],
                    color='g',
                    alpha=0.1)
axh.set(xlabel='Frequency (Hz)',
        ylabel='Filter Freq. Response (db)')
axh.grid()

# power spectrum
axh = ax[2]
yUB = np.max(np.fmax(psd,psdF))
yLB = np.min(np.fmin(psd,psdF))
axh.plot(fout1, psd)
axh.plot(fout2, psdF)
axh.fill_betweenx([yLB,yUB],cutoff,fout1[-1],
                    color='g',
                    alpha=0.1)
axh.set_yscale('log')
axh.set(xlabel='Frequency (Hz)',
          ylabel='Power spectrum')
axh.legend(["Original","Original Filt.","Passband"],loc="lower right")
axh.grid()

plt.savefig("output.png",format="png")
plt.show()

问题原因分析 #

1. 过高的滤波器阶数
   你使用了100阶的Butterworth高通滤波器，虽然Butterworth滤波器以通带最大平坦为特点，但阶数过高时，通带内靠近截止频率（0.1Hz）的区域，幅度响应会出现明显衰减。从Butterworth幅频响应公式来看：
   |H(jω)| = 1 / sqrt(1 + (fc/ω)^(2n))
   其中n为滤波器阶数，fc为截止频率，ω为信号角频率。当ω接近fc时，(fc/ω)接近1，随着n增大，这个项对幅度的影响会被放大，导致靠近截止频率的通带信号出现可见衰减。

2. 零相位滤波的二次放大效应
   你使用了sosfiltfilt进行滤波，该函数会对信号进行正向和反向两次滤波，以实现零相位失真，但代价是幅频响应变为原滤波器的平方。这意味着通带内的衰减会加倍，进一步放大了高阶滤波器带来的幅度损耗。

解决方案 #

• 降低滤波器阶数：无需使用100阶，建议尝试10-20阶，既能保证足够的阻带滚降速率，又能避免通带内过度衰减。
• 调整截止频率：如果信号中需要保留的分量集中在远高于0.1Hz的频段，可以适当提高截止频率，让需要保留的信号远离衰减区域。
• 更换滤波器类型：如果需要更严格的通带平坦度，可以考虑使用Chebyshev I型滤波器（允许通带微小波纹，换取更陡的滚降），或根据需求选择合适的滤波器类型。

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