ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

Python数字滤波器实战:参数设计、scipy实现与边界效应排查

Python数字滤波器实战:参数设计、scipy实现与边界效应排查 滤波器这个词看起来很正经做起来却很考验心态。数字滤波器用来处理一维时间序列把不需要的频率成分压掉、保留有用的部分听起来只是“几句话的事”。可真当自己拿到一段脏数据调完一个低通滤波之后发现输出还是乱、波形开头结尾还会跳盯着图看半天的瞬间脑子里确实会蹦出三个字我的刀呢当然不是真去动刀那是夸张说法。更准确的状态是想深呼吸不太想砸电脑但也确实想找点东西发泄一下。我处理信号数据这些年凡是滤波出问题最后回想起来大多不是滤波器本身不行而是使用方式出了偏差。这篇文章就按我的实际落地流程把数字滤波器的类型选择、参数设计、Python 实现和排查顺序完整拆一遍。不管你是刚开始用 scipy 处理信号的新手还是偶尔要清理传感器数据、音频信号、振动曲线的工程师都可以把这套流程当作一份能照着操作的经验记录。先说个总体结论滤波要稳定跑通不只是填一个butter()函数那么简单。你需要先把采样率、目标频率、滤波器类型、阶数和处理方向这五件事定清楚。下面按顺序展开。1. 滤波器看着不难实际动手后的怪象太多了1.1 最常见的三类“崩溃现场”第一类现场是低频漂移处理翻车。采集到的数据往往自带基线漂移传感器在桌面上放了一夜温度变化、结构松动、零点迁移都会让波形整体慢慢上下走。你会很自然地想到用高通滤波把趋势去掉。结果滤波器跑完之后漂移确实小了一些但信号两端出现一大段上下乱甩的振铃原本正常的起始段被改得面目全非。这个时候人最容易怀疑参数写错。第二类现场是工频干扰。实验室或者工业现场录到的信号经常带 50Hz 或者 60Hz 的周期性噪声。按理说处理这种固定频率干扰陷波器最合适。可实际操作时有人用带通、有人用低通还有人想用移动平均硬拖。一顿操作以后波形上仍然能看到周期性毛刺或者有效信号被削掉一块。问题不在“要不要滤波”而在“用什么形状的滤波响应去切”。第三类现场是相位偏差。你只想对比两个信号里特征的先后顺序直接用因果滤波器跑完波形整体会延迟一段时间。眼睛粗看感觉差别不大但一旦去定位峰值、起点、过零时刻就会发现时间轴对不上。有时候延迟只有几十毫秒但对于高频振动或者瞬态事件分析这点偏差已经足以让结论完全跑偏。这三个现场归根结底反映的是同一个问题滤波不是“选一个函数点运行”。它背后隐藏着采样频率、频率归一化、滤波器类型、阶数、离线或在线处理方式等一系列条件。任何一个条件没有和实际信号对齐输出结果都会变得很怪。1.2 滤波器到底能解决什么问题把话说白数字滤波器描述的是一个线性系统。输入一个离散序列输出另一个离散序列区别在于不同频率成分的幅值被放大或衰减同时相位也会被改变。低通滤波器让低频通过、压掉高频高通滤波器相反带通只保留中间一段带阻则把某一小段频率挖掉。实际场景中它主要解决四类问题。第一去掉已知的周期性干扰。最常见的就是市电工频噪声50Hz 或 60Hz 这种固定频率适合用陷波器定点处理。第二降低随机噪声。传感器本身的底噪、采集链路上的随机扰动如果集中在高频段低通滤波可以把它们压下去。第三把研究关心的频段单独切出来。比如分析旋转机械振动时只想看与转速相关的几百赫兹频段就不需要把全频段都保留。第四做基线漂移修正。数据整体有一个缓慢变化的趋势项不在你关注的频率范围内用高通滤波把这个趋势先拿掉。但边界也很明确。滤波器不是信号增强器它不会把缺失的频段补回来。如果采样率不够信号已经发生混叠那么之后再怎么滤都不可能把真实信息恢复出来。如果传感器前端已经失真采集电路已经饱和滤波器同样无能为力。我见过太多人调了半天滤波参数最后发现原因是硬件采集配置有误高频信号早已混叠到低频段。滤波器是个好工具但它的前置条件是输入信号本身没有坏到不可挽救。2. 动手前先量信号四个关键参数别拍脑袋2.1 采样率决定一切频率参数的坐标任何数字滤波器的频率设置都要先落到采样率这个坐标系里。离散信号能表示的最高频率是采样率的一半也就是奈奎斯特频率。你的截止频率必须小于这个值同时还要用截止频率除以奈奎斯特频率得到归一化频率后再传给滤波器设计函数。举个例子。如果采样率fs 1000Hz奈奎斯特频率是500Hz。你想保留 30Hz 以下的信号归一化截止频率就是30 / 500 0.06。如果采样率是2000Hz同样的 30Hz 目标归一化频率就变成了30 / 1000 0.03。很多人踩坑就踩在这里把截止频率直接除以采样率而不是除以奈奎斯特频率。结果滤波器实际截止频率比预期高了一倍高频噪声自然压不掉。还有一类情况是根本不清楚采样率。文件来自某个采集软件采样率写在配置里没同步过来。这时你用时间轴的单位算频率算出来的物理含义是错的。我的经验是拿到数据第一步先确认采样率、单位、通道顺序没有这些信息先不要写任何滤波代码。2.2 先确定要保留哪个频率区域滤波器的类型选择不能靠“看着不顺眼就滤一下”。你需要先大致知道信号里都有哪些成分以及哪一种成分是干扰。我这里给出一张常用选型表场景推荐滤波器对应 btype数据里高频毛刺多要保留整体趋势低通lowpass基线漂移要移除缓慢趋势高通highpass只关心某个频段比如 10Hz 到 100Hz带通bandpass某个固定频率干扰比如 50Hz 工频陷波或窄带带阻bandstop只想去掉某一段固定频率其余保留带阻bandstop频率区域定下来之后还要看一下目标频段和干扰频段离得有多近。如果目标信号是 45Hz 到 55Hz而干扰恰好是 50Hz这种情况下要非常谨慎。陷波器虽然能挖掉 50Hz但如果带宽设置太宽会把 45Hz 到 55Hz 的有效成分一起削弱。如果目标频段和干扰频段靠得很近但又必须做滤波不要迷信一阶两阶的简单滤波器。你需要更陡的过渡带或者更大的滤波器阶数同时还要接受它带来的相位失真和边缘效应。2.3 阶数不是越大越好默认可以从 4 阶开始很多人调滤波器时第一反应就是把阶数调高觉得阶数越高衰减越快效果越好。这个想法方向是对的但不是没有代价。阶数越高过渡带越窄阻带衰减越快这是好处。但代价同样明显相位失真变大群延迟明显时域波形畸变更大数值上也更容易出现不稳定。特别是在用直接型结构传递系数时高阶滤波器容易出现数值精度问题输出结果看起来像在震荡。所以在代码里我会优先使用二阶节形式也就是outputsos。这种结构把高阶滤波器拆成多个二阶环节数值稳定性比直接传b, a更好批量处理长序列时也更可靠。新手第一次做我一般建议从 4 阶开始。先用 4 阶跑通看效果如果频率响应不够陡再逐步增加到 6 阶、8 阶。每次只加一点点并且对比滤波前后的波形。另外要注意离线做零相位滤波时前向和后向各跑一遍实际等效阶数是原来的两倍。4 阶的 Butterworth 滤波器经过filtfilt处理后实际频率响应接近 8 阶但相位是零。这一点在设置时要提前想清楚。3. 在 Python 里把常用滤波器逐个跑通3.1 先准备一套最小环境做数字滤波最常见的工具组合是 Python 加 NumPy、SciPy配合 Matplotlib 画图。下面这些导入语句足够覆盖大部分场景import numpy as np from scipy import signal import matplotlib.pyplot as plt我的建议是先在虚拟环境里安装 scipy 和 matplotlib再跑下面的样例。如果你本机 scipy 版本较新函数签名可能有一些扩展参数如果代码报错优先看scipy.signal里的函数说明。不要照着网上老版本的代码直接抄版本差异本身就是一个常见坑点。为了验证效果我会先构造一段合成信号。合成信号的好处是你要滤掉的频率和目标频率是已知的结果对不对一眼就能看出来。等合成信号验证通过再换成真实数据。这里用一个 5Hz 目标信号叠加 120Hz 干扰和随机噪声的例子fs 1000.0 # 采样率 t np.arange(0, 5, 1.0 / fs) # 5 秒时间 rng np.random.default_rng(0) x (np.sin(2 * np.pi * 5.0 * t) 0.3 * np.sin(2 * np.pi * 120.0 * t) rng.normal(0, 0.05, len(t)))这样一段信号理想低通滤波后应该能明显看到 5Hz 正弦波而 120Hz 分量被压到很低。3.2 低通滤波器把高频干扰压掉低通是最常用的滤波器。下面是 Butterworth 低通的典型写法cutoff 30.0 order 4 sos signal.butter(order, cutoff / (fs / 2.0), btypelowpass, outputsos) y signal.sosfilt(sos, x)这里sos是二阶节滤波器系数。用outputsos而不用默认的ba是因为高阶 IIR 滤波器直接返回b, a在数值上可能不稳定而二阶节结构更稳。尤其当阶数调高或者数据量很大时这个区别更明显。只跑一次sosfilt是因果滤波输出会有相位延迟。如果是在采集过程中实时处理数据只能用它。但如果已经采集完一整段数据放在本地做离线分析我会优先使用零相位滤波y_zero signal.sosfiltfilt(sos, x)filtfilt会先正向滤波一次再把结果反向滤波一次抵消相位偏移。代价是两端会有边界效应数据很短时尤其明显。使用时要先看两端是否出现异常跳变。3.3 高通滤波器把基线漂移移走处理基线漂移高通滤波器很常见。下面这段构造了一个带有线性漂移的测试信号x_drift (np.sin(2 * np.pi * 5.0 * t) 2.0 * t / t[-1] rng.normal(0, 0.02, len(t))) cutoff_hp 0.5 order 4 sos_hp signal.butter(order, cutoff_hp / (fs / 2.0), btypehighpass, outputsos) y_hp signal.sosfiltfilt(sos_hp, x_drift)高通滤波的截止频率要结合信号的物理含义来定。如果基线漂移变化很慢比如周期超过 10 秒那么截止频率设在 0.1Hz 左右通常够用。如果你把截止频率设成 5Hz那低于 5Hz 的有效低频成分也会一起被去掉滤波后的波形可能看不出原来的低频轮廓。高通滤波做零相位处理时两端跳变往往比低通更明显。因为信号里如果存在很长的趋势项前向滤波开始阶段会因为历史状态未知而产生一段过渡。解决办法有几种延长数据头尾再切掉用padlen参数增加边界填充或者检查是否真的需要零相位。不管用哪种输出后都要单独看前 1 秒和后 1 秒的波形。3.4 带通和带阻精确切出频率窗口如果只关心某个频率区间用带通。比如 10Hz 到 50Hzlow_cut 10.0 high_cut 50.0 sos_bp signal.butter(order, [low_cut / (fs / 2.0), high_cut / (fs / 2.0)], btypebandpass, outputsos) y_bp signal.sosfiltfilt(sos_bp, x)注意两个边界频率都要除以奈奎斯特频率。这个列表里的数值会映射成两个归一化频率点。如果漏掉除法或者低边界大于高边界函数会直接报错或者输出明显不对的系数。带阻与带通相反用于切除一段固定频率。比如要切掉以 50Hz 为中心的窄带干扰除了陷波器也可以用窄带带阻近似。不过 Butterworth 带阻的过渡带相对宽如果只想去掉固定工频陷波器更精准。3.5 陷波器定点处理固定频率干扰陷波器专门对付某一固定频率。最典型的场景是 50Hz 工频干扰。下面用iirnotch实现注意这里传入的是归一化频率f0 50.0 Q 30.0 b_notch, a_notch signal.iirnotch(f0 / (fs / 2.0), Q) y_notch signal.lfilter(b_notch, a_notch, x)Q值决定陷波的带宽。Q 越大陷波越窄只会切除 50Hz 附近的少量频率Q 越小陷波越宽可能会把相邻频率的有用信息一起削弱。先从一个适中的 Q 开始比如 20 到 50再根据频谱观察效果调整。如果 50Hz 附近毛刺非常宽可能不是单频干扰而是别的问题需要先看频谱再决定。另外iirnotch的函数签名在一些旧版本里可能没有fs参数。如果你使用的是新版本也可以直接写成signal.iirnotch(50.0, Q, fsfs)本质上和归一化是同一个意思。遇到接口报错就先用归一化写法通常更通用。3.6 验证频率响应和滤波结果滤波代码能跑不等于滤波结果正确。我习惯从两个角度验证。第一画滤波器的频率响应。以低通为例b_check, a_check signal.butter(order, cutoff / (fs / 2.0)) w, h signal.freqz(b_check, a_check, worN8192) freq_hz w / np.pi * (fs / 2.0) amp_db 20 * np.log10(np.abs(h)) plt.semilogx(freq_hz, amp_db) plt.xlabel(Frequency (Hz)) plt.ylabel(Gain (dB)) plt.show()这条曲线能直观看出截止频率是否准确阻带衰减够不够。一般来说截止频率位置对应 -3dB也就是幅度降到 0.707 左右。如果你发现曲线形状不对先别怀疑绘图的freqz回到归一化频率的计算上检查。第二对比滤波前后的频谱。直接看时域波形容易产生错觉频谱图更客观。你可以用np.fft.rfft计算幅度谱看目标频段的能量是否保留、干扰频段是否明显降低。如果滤波后目标频率也变弱了说明参数选择过于激进要调低阶数或放宽频段范围。4. 输出不对时先按这条链路排查4.1 常见现象和对应原因滤波结果异常现象通常分几类。我先把常见对应关系列出来现象优先排查方向输出和原始信号几乎一样是否忘了给结果变量赋值或截止频率离信号频段太远输出有 NaN 或无穷值输入数据是否包含 NaN 或 inf波形两端上下剧烈跳变滤波边界效应padlen不够或数据太短整段波形明显延迟用了因果滤波sosfilt/lfilter且未做时间对齐滤波后仍有很多高频毛刺截止频率是否按奈奎斯特归一化错误或阶数不足滤完波后有效信号也变了频段窗口太窄过渡带吞掉有效成分程序直接报维度错误输入是二维数组滤波方向或axis参数没指定对这些现象里最容易被误判的是“滤波后出现 NaN”。很多人在这一步反复调滤波阶数甚至还换滤波器类型结果发现原始数据里有几个空洞值。滤波器遇到 NaN 会把空洞一路传播出去处理多少样点都躲不掉。所以排查的第一步永远先看输入数据。4.2 先看数据再看参数最后才怀疑代码我的排查顺序是固定的。第一步检查输入数据。看长度、维度、NaN、inf、常值段。常值段很致命一段电平完全不变的数据在频率上等价于极低频甚至直流。如果用高通滤波常值段会被当成低频趋势去掉结果在常值段边缘产生跳变。数据清理要先处理缺失值和异常跳点不要直接丢给滤波器。第二步画原始数据的频谱。不知道信号里有哪些频率成分就没有办法判断滤波结果对不对。画fft幅度谱以后你会看到几个明显峰低频基线对应最低频的峰工频干扰对应 50Hz 附近目标信号对应自己的频段。先把这些峰标出来写滤波参数时心里就有数了。第三步检查归一化频率。这是 IIR 滤波器最常见的参数错误。在代码里加一行打印print(Nyquist:, fs / 2.0) print(Wn:, cutoff / (fs / 2.0))如果打印出来的Wn不在 0 到 1 之间说明采样率或者截止频率写错了。这个问题发生过很多次尤其当数据来自不同采样率的采集设备时最容易搞混。第四步才看滤波代码本身。用一段已知的合成信号做测试比如叠加一个 10Hz 正弦和一个 500Hz 正弦用 50Hz 截止频率的低通去滤。如果合成信号的测试没过说明代码或参数有问题如果合成信号过了但真实数据仍然很怪那问题大概率在数据质量或者真实信号与目标频段不符合假设。4.3 因果滤波和零相位滤波不能混着用这是很多人容易忽略的一个点。sosfilt和lfilter是因果滤波器适合实时在线处理但会引入相位延迟。sosfiltfilt和filtfilt是零相位滤波适合离线整段处理但不可用于实时流式数据处理。如果你在离线分析时用了sosfilt然后拿滤波结果去和原始信号做特征点对比波形会因为延迟而错位。解决方法是改用sosfiltfilt或者对齐延迟量。反过来如果系统是做在线数据处理的每一帧数据到了就要立即输出结果那就不该用filtfilt。它需要整段数据才能反向滤波在流式场景下要么延迟极大要么根本无法工作。这种情况下应该用因果滤波器并接受一定的相位延迟或者设计一个相位响应已知的滤波器后续再做延迟补偿。很多项目里“滤波结果看起来对但系统就是不稳定”本质上就是这两类滤波方式被用错了场景。5. 批量处理时最怕一个参数打天下5.1 先跑两条别一上来就整个目录跑完如果只是处理单个文件事情其实好办。真正容易出问题的是批处理目录里有几百个数据文件脚本写好以后双击运行过一会儿看输出目录好像都生成了文件。但如果你没人工检查很可能整批结果都带着同样的错误。我见过最典型的场景是一批数据来自多个采样率的采集设备但脚本里写死了fs 1000。结果采样率是 500Hz 的文件实际奈奎斯特频率是 250Hz脚本却按 500Hz 来算截止频率滤波行为完全错误。程序没有报错因为所有数值计算都合法但输出已经不能用。所以批量处理的第一原则是先挑两个有代表性的文件一条一条跑跑完人工看图。确认波形、频谱、两端边界都正常以后再放开整批跑。文件数量越多越要谨慎。5.2 把参数集中管理不要散落在循环里批量滤波时我会把滤波器参数集中写在文件开头或单独配置里。这样调整参数时不需要在几百行循环代码里翻来翻去。下面是一个示例结构from pathlib import Path config { fs: 1000.0, order: 4, low_cut: 0.5, high_cut: None, } if config[high_cut] is None: cutoff config[low_cut] btype highpass else: cutoff [config[low_cut], config[high_cut]] btype bandpass sos signal.butter(config[order], np.asarray(cutoff) / (config[fs] / 2.0), btypebtype, outputsos)然后用循环处理每个文件。循环里至少要做三件事记录输入文件名、记录异常、只在成功时写输出文件。raw_dir Path(./raw) out_dir Path(./filtered) out_dir.mkdir(exist_okTrue) for path in sorted(raw_dir.glob(*.npy)): data np.load(path) filtered signal.sosfiltfilt(sos, data, axis0) out_name out_dir / (path.stem _filtered.npy) np.save(out_name, filtered) print(done:, path.name)5.3 为每个文件的采样率单独设计滤波器如果整批数据的采样率相同可以用同一个sos系数。但采样率不同时必须在循环里分别设计滤波器不能复用外层那个sos。原因是滤波器的频率参数是物理频率转换到归一化频率后的结果。1000Hz 采样率下截止 30Hz 的归一化频率是 0.06换到 500Hz 采样率同样截止 30Hz归一化频率就变成了 0.12。复用同一个sos本质上是在不同的频率坐标系里执行同一个操作结果自然不会一样。正确做法是在每次读取文件时从文件名或元数据里拿到真实采样率再重新调用一次signal.butter。多花一点计算时间但能避免整批数据处理结果失真。批量跑完后不要只看输出文件数量对不对还要随机抽查 3 到 5 个文件。比较原始频谱和滤波后频谱确认没用错采样率也没有把有效频段滤没。如果能自动生成一个汇总图比如每个文件滤波前后频谱的小缩略图排查效率会高很多。6. 把滤波调到“不那么想找刀”的几个收尾习惯6.1 调参时每次只改一个变量滤波参数之间存在联动。阶数影响过渡带截止频率影响保留频段滤波器类型影响整个频响形状。如果同时改三个参数结果变化时你根本不知道是哪一个起了作用。我自己调试时会从一组保守参数开始比如 4 阶 Butterworth 低通截止频率先放在干扰频段和目标频段之间。然后一次只改一个参数先看截止频率是否合适再决定要不要提高阶数。每改一次就看一遍时域波形和频谱留下记录。这个流程看起来慢实际上是最快的。比“随机组合参数”快得多。如果手里文件非常多需要批量试参数就把参数组合写成字典列表逐个组合跑同一份样例数据。只要样例数据选得有代表性自动筛选也能帮你缩小范围。但最终选哪组参数还是要靠人工看输出波形不能只看某个数值指标。数值指标只是辅助信号处理最终要回归到你关心的物理特征是否被正确保留。6.2 落地前做一套固定
返回列表