如何改进Scipy.signal.butter带通巴特沃斯滤波器曲线开头结果？

一、开头失真的根源 #

你遇到的开头失真是线性滤波器的暂态响应：滤波器初始状态为0，处理信号时需要一定时间收敛到稳态，这段时间内的输出会偏离真实信号，这是巴特沃斯滤波器的固有特性。

二、改进开头失真的实用方法 #

1. 零相位滤波（离线处理首选） #

用scipy.signal.filtfilt替代常规的lfilter，它会对信号正向滤波后再反向滤波，既抵消相位偏移，又让开头的暂态被反向处理的信号覆盖，几乎消除开头失真。适合电流测量这类离线后处理场景。

from scipy.signal import butter, filtfilt
import numpy as np

def butter_bandpass_filtfilt(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a

# 调用示例：假设raw_current是你的电流测量数据
b, a = butter_bandpass_filtfilt(lowcut=10, highcut=100, fs=1000, order=4)
filtered_current = filtfilt(b, a, raw_current)

2. 信号前置延拓 #

在原始信号开头添加一段延拓数据，让滤波器在处理真实数据前进入稳态，之后切除延拓部分。常用延拓方式包括镜像延拓、平稳段复制：

from scipy.signal import butter, lfilter

# 先定义常规带通滤波器
def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a

b, a = butter_bandpass(lowcut=10, highcut=100, fs=1000, order=4)
raw_current = np.random.randn(1000)  # 替换为你的真实数据

# 镜像延拓开头：取前N个点的镜像，N取滤波器阶数的5-10倍
N = 40  # 4阶滤波器取40个点
extended_data = np.concatenate([raw_current[:N][::-1], raw_current])
filtered_extended = lfilter(b, a, extended_data)
# 切除延拓部分，得到修正后的信号
filtered_current = filtered_extended[N:]

3. 降低滤波器阶数 #

高阶滤波器的暂态响应更长（极点更多，收敛慢），如果对滤波陡度要求不极端，可适当降低阶数（比如从5阶降到3阶），缩短失真区域长度。

三、失真区域位置的计算 #

经验估算 #

巴特沃斯滤波器的暂态长度（样本数）大约为 (3~5)*滤波器阶数（当截止频率远低于采样频率时）；若要换算成时间，再除以采样频率即可。比如4阶滤波器、采样频率1000Hz，暂态长度约为12-20个样本，对应0.012-0.02秒。

精确计算（通过脉冲响应） #

计算滤波器的脉冲响应，找到响应衰减到稳态最大值1%以内的位置，这个位置就是失真区域的结束点：

from scipy.signal import impulse
import numpy as np

# 用之前得到的滤波器系数b,a
impulse_resp, _ = impulse((b, a), n=2000)  # n取足够大的样本数
max_resp = np.max(np.abs(impulse_resp))
# 找到第一个响应绝对值小于最大响应1%的索引
distortion_end_idx = np.argmax(np.abs(impulse_resp) < 0.01 * max_resp)
# 失真区域为信号的0到distortion_end_idx索引范围
print(f"失真区域结束于第{distortion_end_idx}个样本")

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