如何用Scipy为时间序列样本添加高低通滤波器并解决频率参数报错

错误原因 #

当你在scipy.signal.butter中指定fs参数时，Wn需要传入实际物理频率值，且必须小于采样率的一半（即Nyquist频率）。你的采样率是0.1 Hz，Nyquist频率为0.05 Hz，而你设置的Wn=0.1已经等于采样率，超过了Nyquist频率的上限，因此触发错误。

修正步骤 #

1. 调整Wn值：将Wn设置为小于0.05 Hz的实际频率（比如0.02 Hz，根据你的需求选择合适的截止频率）。
2. 使用正确的滤波函数：不要直接用np.convolve应用滤波器系数，应该用scipy.signal.lfilter——Butterworth返回的ba系数是为线性滤波器设计的，lfilter能正确处理IIR滤波器的递归结构，避免边界和滤波逻辑错误。
3. 分别处理前后样本：拆分数据集，对前1000个样本做低通滤波，后1000个做高通滤波，再合并回原数据。

修正代码示例 #

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

# 假设你的原始数据已经加载：time_data_ISLL_11_21, yf_ISLL_11_21_irfft
fs = 0.1  # 采样率1/10 Hz
nyquist = fs / 2  # Nyquist频率0.05 Hz

# 1. 对前1000个样本做低通滤波
low_cutoff = 0.02  # 实际截止频率，小于0.05 Hz
b_low, a_low = signal.butter(2, low_cutoff, btype='low', analog=False, output='ba', fs=fs)
y_lowpass = signal.lfilter(b_low, a_low, yf_ISLL_11_21_irfft[:1000])

# 2. 对后1000个样本做高通滤波
high_cutoff = 0.01  # 实际截止频率，小于0.05 Hz
b_high, a_high = signal.butter(2, high_cutoff, btype='high', analog=False, output='ba', fs=fs)
y_highpass = signal.lfilter(b_high, a_high, yf_ISLL_11_21_irfft[-1000:])

# 3. 合并滤波后的样本与原始中间数据
y_filtered = np.copy(yf_ISLL_11_21_irfft)
y_filtered[:1000] = y_lowpass
y_filtered[-1000:] = y_highpass

# 绘图查看结果
plt.plot(time_data_ISLL_11_21, y_filtered)
plt.show()

额外说明 #

• 如果不想用实际频率，也可以不指定fs参数，此时Wn需要是归一化到Nyquist频率的值（范围0到1）。比如对于0.02 Hz的截止频率，归一化后是0.02 / 0.05 = 0.4，此时代码可以写成：signal.butter(2, 0.4, btype='low', analog=False, output='ba')。
• 若需要更平滑的边界处理，可以考虑使用signal.filtfilt（零相位滤波），但注意它会处理数据两次，可能改变信号相位特性（如果你的傅里叶变换对相位敏感，需谨慎使用）。

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