4.4 实战：信号滤波与频谱分析

4.4 实战：信号滤波与频谱分析

本节摘要：本节把前三节串成一条完整流水线：一段含 50 赫兹工频干扰的振动监测信号，先用 rfft 频谱分析定位干扰，再用带阻与低通滤波器组合去噪，用 filtfilt 保证零相位，最后用频谱对比验证效果；后半部分对一张带噪图像做中值与高斯滤波对比。跑完本节，你将得到一条可复用的"分析—设计—验证"标准流程。

本节地图

阅读完本节，你应当能够：

1. 独立完成"频谱定位 → 滤波器设计 → 滤波 → 频谱验证"的完整流程；
2. 根据干扰特征在带阻、低通之间做出组合决策；
3. 用滤波前后的频谱差异定量评估效果；
4. 对比中值与高斯滤波在图像去噪上的表现；
5. 形成"先分析、再设计、后验证"的工程习惯。

一、场景与数据

假设你在做设备状态监测。传感器采到的振动信号里，主体是转轴的 10 赫兹和 25 赫兹振动，但现场供电线路的 50 赫兹工频干扰串了进来，还叠加了随机噪声。先造一份模拟数据，把问题复现出来：

import numpy as np
from scipy import signal
from scipy.fft import rfft, rfftfreq

fs = 1000
t = np.arange(0, 4, 1/fs)
rng = np.random.default_rng(42)
x = (np.sin(2*np.pi*10*t) + 0.6*np.sin(2*np.pi*25*t)
     + 1.2*np.sin(2*np.pi*50*t) + 0.15*rng.normal(size=len(t)))
​

注意工频干扰的幅度是 1.2，比 10 赫兹和 25 赫兹的幅度都大。画时域波形的话，整条曲线被 50 赫兹的波浪主导，10 赫兹和 25 赫兹的成分几乎看不出来。这正是先做频谱分析的理由：时域看不清的问题，频域一眼定位。

先说两个参数为什么这么定。工频干扰来自电网：中国市电是 50 赫兹，北美是 60 赫兹，仪器一旦接地不良或屏蔽不够，这个频率就会串进测量链路，所以它是振动监测里最常见的干扰源。采样率取 1000 赫兹，是因为奈奎斯特频率 500 赫兹远高于我们关心的 30 赫兹频带，留足了裕量；真实项目里，如果信号里可能存在更高频的成分，采样前还要加抗混叠低通，这是硬件层面的功课。数据时长取 4 秒，频率分辨率就是 0.25 赫兹，足够把 10、25、50 三个频率看得清清楚楚——这三个数字都不是随手写的，背后各有一笔账。

二、第一步：频谱定位

对原始信号做 rfft，得到幅值谱：

X = rfft(x)
freq = rfftfreq(len(x), 1/fs)
​

频谱上 50 赫兹处会出现一个明显高于其他成分的尖峰，10 赫兹和 25 赫兹的峰被它盖住大半，随机噪声则铺满整个频带。这一眼能看出三件事：干扰频率精确已知，需要衰减的频带很窄，适合带阻；有用成分集中在 30 赫兹以下，低通截止频率可以放在 40 赫兹附近；随机噪声是宽带信号，只能靠低通整体压制。方案随之清晰：先带阻干掉 50 赫兹，再低通把 30 赫兹以上的噪声压下去。

读频谱时顺便核对一下幅值：50 赫兹处峰值接近 1.2，10 赫兹处接近 1.0，25 赫兹处接近 0.6，和构造参数一一对应——这说明我们的频谱分析读数是可信的。养成这个习惯：任何频谱分析，先拿已知成分对答案，再信它。对不上就回去查时间轴、查归一化，别急着往下走。这一核对动作的成本不到一分钟，却能让后面所有结论都站得住脚。

三、第二步：滤波器设计

工频干扰频带窄，用带阻最合适。scipy.signal 提供 iirnotch，专门为"掐掉单一频率"设计：

b_notch, a_notch = signal.iirnotch(50, 30, fs)   # 50 赫兹，品质因数 30
​

第二个参数是品质因数 Q，控制凹陷的宽度。工频在实际电网里有轻微漂移，Q 取 20 到 40 比较稳妥：太大，凹陷太窄，频率一漂就掐不中；太小，凹陷太宽，会误伤 40 赫兹附近的信号。之后设计 40 赫兹低通，把残余高频噪声清掉：

sos_lp = signal.butter(4, 40/(fs/2), btype='low', output='sos')
​

低通用 sos 形式，保证数值稳定。两个滤波器，一个解决窄带干扰，一个解决宽带噪声，分工明确。

Q 值具体怎么影响效果？品质因数 30 的带阻，负 3 分贝带宽大约是中心频率除以 Q，也就是 1.7 赫兹左右：50 赫兹处被掐掉，49 赫兹和 51 赫兹基本不受影响。低通阶数取 4 是折中——2 阶过渡带太缓，25 赫兹附近会受损；6 阶以上过渡带更陡，但相位畸变明显，波形变形更重。选参数的本质是在"砍得干净"和"伤得少"之间找平衡点，而平衡点的位置要靠频谱验证来确认。

三种方案的取舍可以先想清楚再动手：

四、第三步：滤波与验证

离线分析用 filtfilt，保证零相位，波形的峰值位置不偏移：

y_notch = signal.filtfilt(b_notch, a_notch, x)
y = signal.sosfiltfilt(sos_lp, y_notch)
​

验证不能靠"看起来干净了"，要回到频域对比。滤波前后分别做 rfft，重点看 50 赫兹处的幅值：理想情况衰减 40 分贝以上，也就是幅值缩到原来的百分之一。同时对比时域波形：10 赫兹和 25 赫兹的周期特征清晰可见，波形没有整体偏移，说明零相位处理起作用了。这一步做下来，去噪效果是数字说话，而不是感觉说话。

⚠️ 常见坑：iirnotch 的第三个参数是采样率，它会在内部完成频率归一化，别再手动除以奈奎斯特频率。传错采样率不会报错，但凹陷会落在错误的位置，频谱验证时一眼就能看出"掐错了频率"。
💡 关键直觉：验证优先看频域。滤波效果在时域上永远"看起来不错"，只有回到频谱，才能确认每个频率成分都按预期被处理——这套分析、设计、验证的闭环，比任何一个单个函数都值钱。

顺带提醒：filtfilt 在数据两端有边缘效应，验证时把首尾各丢 0.2 秒再看，避免把边缘假象当成功效。真实项目里，这两小段本来也不该进入统计。定量指标再加一个：滤波前后的信噪比。用 10 赫兹处的幅值代表信号强度，用 100 到 400 赫兹区间的平均幅值代表噪声底，两者之比就是信噪比。滤波后这个比值通常能提升 20 分贝以上，配合 50 赫兹处的 40 分贝衰减，两份数字合起来才是完整的验收报告。时域波形也要看：对比滤波前后的前一秒，10 赫兹的慢波与 25 赫兹的快波叠加形态应该清晰可辨，且波形起点与原始信号对齐，没有整体平移——对齐这件事，就是零相位处理的直接证据。

五、第四步：图像去噪对比

信号处理完，再来图像。构造一张带高斯噪声和椒盐噪声的合成图，分别用 gaussian_filter 和 median_filter 处理，用均方误差量化谁更接近原图：

from scipy import ndimage

img = np.zeros((256, 256))
img[80:180, 60:200] = 1
noisy = img + 0.3*rng.normal(size=img.shape)
noisy[rng.random(img.shape) < 0.08] = 1
noisy[rng.random(img.shape) < 0.08] = 0

g = ndimage.gaussian_filter(noisy, sigma=1.5)
m = ndimage.median_filter(noisy, size=3)

mse_g = ((g - img)**2).mean()
mse_m = ((m - img)**2).mean()
​

结果通常显示中值滤波在这张图上略胜一筹：椒盐坏点是极端的离群值，高斯滤波的加权平均会被它们拖偏，而中值滤波天然免疫。这再次验证了 4.3 的结论：先判断噪声类型，再选工具，顺序不能反。如果你拿到的是只有高斯噪声的图，结论就会反过来——所以这个对比实验值得亲手跑一遍，感受"工具跟着噪声类型走"的分量。

如果要做正式汇报，均方误差不够直观，可以换算成峰值信噪比（PSNR）：峰值信号的平方除以均方误差，再取 10 倍对数。PSNR 每高 6 分贝，视觉质量就上一个台阶，中值滤波在这张图上的 PSNR 通常会高出几个分贝。另外提醒一句：合成图的原图是已知的，所以能做定量对比；真实照片没有原图，只能靠目测加领域知识判断——这也是为什么实验数据要留一份"标准答案"。图像去噪的结论和信号滤波是相通的：工具的选择由噪声性质决定，参数由数据决定，验证由指标决定。到这里，信号与图像两条线在本章正式汇成一条。

六、流程总结

整条流水线可以画成一张流程图：

这张图的顺序就是执行的顺序，任何一步都可以单独替换：频谱分析可以换 Welch 方法（第 5 章），滤波器可以换 FIR，图像部分可以换成自己的数据。模板的价值在于固定了流程，而不是固定了函数。

选型这一环值得单独画清楚。下面这张图总结了"干扰长什么样，就用什么滤波器"的决策逻辑：窄带干扰找带阻，宽带高频噪声找低通，只留中间频带找带通，低频漂移找高通，复杂场景组合使用、最后统一做零相位和频谱验证。

图 4-4 滤波器选型流程

常见疑问：这个流程能直接复用吗

三个问题几乎每次都会被问到。其一，实时系统怎么办？filtfilt 依赖整段数据，实时场景改用 lfilter，接受相位延迟，或用延迟补偿做修正——实时与零相位不可兼得，这个取舍要在方案阶段就定下来。其二，换一组参数效果差多少？把 Q 从 30 改到 10，凹陷变宽，40 赫兹附近会开始受损；把低通截止从 40 改到 35，25 赫兹成分安然无恙，噪声更干净。参数是连续谱，多跑几次对比就熟了，不必背。其三，音频能这么处理吗？能，但注意两件事：音频的采样率通常是 44.1 千赫，所有归一化频率按它重新算；音频对相位更敏感，滤波器选型时要更谨慎。流程本身不用改，参数跟着数据走。如果干扰源不止一个（比如同时有 50 赫兹工频和 100 赫兹谐波），用带阻数组把两个频率都陷掉，或者评估后直接上低通，以频谱图上残留噪声的形态为准。

要点串联

• 要点一：干扰特征决定方案——窄带干扰用带阻，宽带噪声用低通，两者常组合使用。
• 要点二：iirnotch 是单频抑制的专用工具，品质因数控制凹陷宽度，Q 取 20 到 40。
• 要点三：离线分析优先 filtfilt，零相位保证波形峰值位置不偏移。
• 要点四：验证必须回到频域，看干扰频率处的幅值衰减了多少分贝。
• 要点五：滤波后的数据两端有边缘效应，分析时丢弃首尾一小段。
• 要点六：图像去噪先判断噪声类型，椒盐噪声用中值，高斯噪声用高斯。
• 要点七：定量指标（均方误差、分贝衰减）比肉眼判断可靠。

第 4 章到此收尾。下一章我们换一批数据：统计分布、空间点云与特殊函数——把信号里学到的"先分析、再动手"带到统计检验与空间算法中。