2.2 频域与小波去噪
2.2 频域与小波去噪
本节摘要：把谱图搬进频率域，信号与噪声的"户口"分居两地——缓变的谱峰住在低频，白噪声摊满全频段，设一道截止线就能分家。小波去噪再给这道截止线装上"局部开关"：峰附近少滤、平坦区多滤。本节用可复算的对照实验给出两种手法的边界与代价。
上一节 SG 平滑治好了大部分毛刺，但它的窗口是全程统一的：窄峰嫌窗口大，平坦区嫌窗口小。这一节换一副眼镜看同一份档案——频率域视图，它承接 2.1 的白噪声假设，通往 2.3 基线校正（基线正是频率域里"最低频的住户"），也为第 6 章的参数调优积累第二套工具。
信号与噪声的户口之争
一条光谱放进傅里叶变换里，会分解成不同周期的波动成分。谱峰是缓变的——一个半高宽几十个点的峰，周期长、频率低；白噪声是逐点独立的随机抖动，周期短、频率高，能量摊满整个频轴。去噪在频率域里就是一道户口线：线以下住信号，线以上住噪声。这与 SG 在"窗口域"里做的取舍是同一件事的两种坐标——SG 窗口宽度本质上就是截止频率的时域表达。
先用一份模拟档案做低通滤波实验：三个高斯峰（半高宽 0.42、0.71、0.59）加 0.03 的白噪声，500 个点：
# 傅里叶低通去噪：截止频率扫描 import numpy as np rng = np.random.default_rng(7) x = np.linspace(0, 10, 500) # 波数轴（任意单位） clean = (0.8*np.exp(-0.5*((x-3.2)/0.18)**2) # 主峰，sigma=0.18 + 0.30*np.exp(-0.5*((x-5.6)/0.30)**2) # 次峰 + 0.12*np.exp(-0.5*((x-8.1)/0.25)**2)) # 弱峰 y = clean + rng.normal(0, 0.03, x.size) # 加 0.03 白噪声 def rmse(a, b): return float(np.sqrt(np.mean((a-b)**2))) print(f"raw RMSE={rmse(y, clean):.4f} peak={y.max():.3f}") Y = np.fft.rfft(y) # 实数傅里叶变换 for k in [8, 15, 30]: # 保留的频率分量个数 Yk = Y.copy(); Yk[k:] = 0 # 高频置零 rec = np.fft.irfft(Yk, x.size) # 逆变换回波数域 print(f"k={k:2d} RMSE={rmse(rec, clean):.4f} peak={rec.max():.3f}")
输出：
raw RMSE=0.0284 peak=0.836 k= 8 RMSE=0.0765 peak=0.460 k=15 RMSE=0.0217 peak=0.713 k=30 RMSE=0.0098 peak=0.789
三档截止线画出一条清晰的边界。k=8 时主峰高度从 0.836 被砍到 0.460——截止线画得太低，把窄峰的高频成分连同噪声一起清了户口，谱峰被"削平"；k=30 时 RMSE 降到 0.0098、峰高保住 0.789，但保留的频率分量已经不少，噪声也按比例混了进来。**全谱统一截止的困境就在这里**：档案里既有需要高频的窄峰，又有只配低频的平坦区，一道线照顾不了两边。
图1 频率域里的户口分布与截止线
小波：给截止线装上局部开关
小波变换的思路是把谱图按尺度逐层拆开：粗尺度（低频）给整体轮廓，细尺度（高频）给细节。关键区别在于小波有位置信息——同样一个高频成分，出现在窄峰肩上算"信号"，出现在平坦区算"噪声"。于是阈值不再一刀切，而是逐系数判断：超过阈值的高频系数保留（此处有真结构），低于阈值的置零或收缩（此处只是噪声）。峰附近自动少滤、平坦区自动多滤，正好补上傅里叶全局截止的短板。
# db4 小波 4 层分解 + 通用阈值软收缩 import pywt coeffs = pywt.wavedec(y, 'db4', level=4) # 1 个近似 + 4 层细节 sigma = np.median(np.abs(coeffs[-1])) / 0.6745 # 最细层系数估噪声（MAD 法） uth = sigma * np.sqrt(2*np.log(y.size)) # 通用阈值 coeffs[1:] = [pywt.threshold(c, uth, 'soft') for c in coeffs[1:]] wd = pywt.waverec(coeffs, 'db4')[:y.size] # 重构 i_main, i_weak = np.argmin(np.abs(x-3.2)), np.argmin(np.abs(x-8.1)) flat = slice(0, 100) # x<2 的无峰平坦区 print(f"wavelet RMSE={rmse(wd, clean):.4f}") print(f"main peak: truth={clean[i_main]:.3f} wavelet={wd[i_main]:.3f}") print(f"weak peak: truth={clean[i_weak]:.3f} wavelet={wd[i_weak]:.3f}") print(f"flat zone noise: raw std={y[flat].std():.4f} wavelet std={wd[flat].std():.4f}")
输出：
wavelet RMSE=0.0134 main peak: truth=0.799 wavelet=0.764 weak peak: truth=0.120 wavelet=0.115 flat zone noise: raw std=0.0263 wavelet std=0.0081
四行输出各有含义。平坦区噪声标准差从 0.0263 压到 0.0081，压掉近七成——平坦区的高频系数几乎全是噪声，阈值下无处可逃；主峰高度 0.799 到 0.764、弱峰 0.120 到 0.115，损失都在 5% 以内——峰附近的高频系数幅度大，跨过了阈值被保留。同一个阈值在两种区域产生了两种行为，这就是"局部开关"四个字的数学实现。代码里还有个值得记住的细节：噪声水平不是拍脑袋给的，而是用最细层系数的中位数除以 0.6745 估计——对正态噪声，这个中位数与标准差有确定的换算关系，噪声参数让数据自己报。
工程取舍
- 小波基与层数：db4、sym8 是光谱里的常客；分解层数 3 到 5 层够用，层数过多会把谱峰本身也拆进"细节"里；
- 通用阈值偏保守：它按"最坏噪声点"定线，保峰但压噪不彻底。想要更狠的去噪可换软阈值为硬阈值，代价是重构谱会出现轻微波纹；
- 脉冲尖峰别用它：宇宙射线打出的单点尖峰在傅里叶域是全频段的，先中值滤波再进本节流程（2.1 的老建议依然有效）；
- 频率域是诊断工具，不只是处理工具：怀疑谱里有 50 赫兹电源纹波这类周期性干扰，先做一次傅里叶变换看谱——窄而高的孤立尖峰一眼现形，比盲调参数快得多。
💡 关键直觉：SG 是"统一窗口"，傅里叶是"统一截止"，小波是"随位置变化的截止"。工具越自适应，可调参数越少、越需要理解它替你做的判断。
本节要点回顾
- 频率域户口：谱峰住低频、白噪声摊满全频段，低通滤波即分家；
- 全局截止的两难：截止过低削窄峰（实验里主峰 0.836 跌到 0.460），过高噪声回潮；
- 小波阈值是局部开关：同一阈值下平坦区噪声压掉近七成、峰高损失小于 5%；
- 噪声水平用 MAD 估计：最细层系数中位数除以 0.6745，让数据自己报噪声；
- 周期性干扰先看频谱再动手：傅里叶变换更是诊断利器。
毛刺处理告一段落，档案上还垫着那条缓缓起伏的斜坡。下一节先学手工派：人工挑出无峰区域，用多项式把基线拟合出来再整条减掉。