★ FEATURED ARTICLE 说实话刚接触“小波分析”这四个字的时候我和很多人一样第一反应是这又是哪个数学分支里的高深术语但真正用起来之后才发现小波分析本质上就是一套“既能看清森林、又能看清树木”的信号处理工具。无论你是做故障诊断、脑电信号处理、地震资料分析还是金融时序预测只要手里攥着一堆非平稳信号想同时知道“什么时间发生了什么频率成分”小波分析就是绕不开的那条路。这篇入门介绍我尽量不堆公式、不甩定理先帮你把小波分析到底在解决什么问题、它的核心思路是什么、上手时最该关注哪几个参数讲清楚。适合刚接触信号处理的研究生、刚转行做算法的工程师以及所有被傅里叶变换“全局频率”坑过的人。读完这篇你至少能看懂小波变换的时频图能用Python跑通一次完整的离散小波分解并且知道该怎么跟别人解释小波分析为什么比傅里叶“更聪明”。1. 为什么傅里叶变换不够用1.1 傅里叶变换只能告诉你“有什么”却说不清“在哪里”经典的傅里叶变换本质上是在做一件事把一段信号拆解成一系列不同频率的正弦波的叠加。对一个平稳信号来说这招非常管用比如一个50Hz的工频正弦波傅里叶变换后频谱上就是一根干净的谱线。但问题来了真实世界里的信号几乎都不平稳。拿机械振动信号举例轴承出现局部缺陷时每一次滚珠路过缺陷位置都会产生一个短促的冲击脉冲这个脉冲在频域上会激起一个很宽的频带。用傅里叶变换看这个信号的频谱你能看到高频成分确实存在但你完全不知道这些高频成分究竟出现在哪个时刻。换句话说傅里叶变换把时间信息彻底“平均”掉了。打个比方傅里叶变换就像你走进一家餐厅只告诉你这桌菜的总热量是2000大卡却不告诉你哪道菜是油炸的、哪道菜是清蒸的。很多工程问题恰恰需要知道“哪道菜有问题”比如故障诊断必须定位到某个时间点出现了冲击语音识别必须知道某个音素在哪个时间段发声。这时候单纯的傅里叶变换就无能为力了。1.2 短时傅里叶变换的妥协与局限为了把时间信息找回来工程师们很自然地想到了一个办法把信号切成一段一段对每一段分别做傅里叶变换这就是短时傅里叶变换STFT。它的思路很直观选定一个固定长度的窗函数在时间轴上滑动每滑动一步就对窗内的信号做一次傅里叶变换最后得到一张“时间-频率”二维谱图。听起来不错但STFT有一个天生的硬伤窗长固定意味着频率分辨率和时间分辨率无法兼得。窗选得越短时间定位越准但频率分辨率越差窗选得越长频率分辨得越清楚但时间上就越来越“模糊”。海森堡不确定性原理在这里给出了一个硬约束时间分辨率和频率分辨率的乘积存在下限你不可能同时把两者都做到无限好。STFT只是在“时间分辨率”和“频率分辨率”之间做了一次固定比例的妥协一旦窗长定了整个时频平面上所有位置的“分辨率”就都一样了。但在真实信号中低频成分往往持续时间长、需要高的频率分辨率高频成分往往是瞬态冲击、需要高的时间分辨率。用同一把尺子去量所有频率成分显然不够聪明。于是小波分析就顺理成章地登场了。2. 小波分析的核心思想多分辨率分析2.1 “变焦镜头”式的观察方式小波分析最核心的思想可以用“多分辨率分析”来概括。它的思路是用一组可伸缩、可平移的基函数去“适配”信号的不同局部特征。低频部分用又宽又平的小波去匹配获得更好的频率分辨率高频部分用又窄又尖的小波去匹配获得更好的时间分辨率。整个时频平面上不同位置自动选取不同形状的观察窗口实现了“自适应变焦”。你可以把小波理解成一把“数学显微镜”想看信号的轮廓时把镜头拉远看到的是整体趋势想看某个瞬间的细节时把镜头拉近看到的是局部突变。而且这个过程是连续可调的你随时可以在“看全貌”和“看细节”之间切换。这就是为什么小波分析特别适合处理那些带有瞬态突变、边缘跳变、局部奇异的信号。2.2 小波基函数怎么来的平移与伸缩小波基函数并不是只有一种而是由一个“母小波”经过平移和伸缩得到的一族函数。母小波是一个均值为0、在有限区间内迅速衰减的波形。把这个母小波记作ψ(t)那么通过对它做平移和伸缩就能生成一族基函数伸缩用尺度因子a控制a越大小波被拉伸得越宽频率越低a越小小波被压缩得越窄频率越高。平移用平移因子b控制决定小波在时间轴上的位置。连续小波变换CWT做的就是把信号和这一族小波基函数逐个做内积得到一组系数C(a,b)。系数越大说明在该尺度、该位置附近信号与这个小波波形的相似程度越高。把所有系数按照“尺度-时间”展开画出来就是一张小波时频图。这里有一个新手比较容易绕晕的点小波的“尺度a”并不直接等于频率它和频率之间是一个近似反比的关系。对一个中心频率为f0的母小波在采样率为fs的信号中尺度a对应的“伪频率”大约是f f0 * fs / a。也就是说尺度数值越小对应的实际频率越高。这个换算关系在做时频图时非常关键很多人在Python里画CWT图发现横轴和纵轴对不上多半就是忘了做这个换算。2.3 连续小波变换到底在算什么连续小波变换的公式写出来其实就是一个内积运算但它的含义得展开说。小波基函数在每一步都会被“挪到”信号的不同位置然后跟信号对应位置的那一段做相似度比较。如果信号局部波形长得像这个小波内积结果就大不像结果就小。实际计算时尺度a从最小到最大遍历一遍平移因子b从信号开头走到结尾最终得到一张二维系数图。这张图的横轴是时间纵轴是尺度或换算后的频率颜色代表系数幅值。举个例子一段含有20Hz和50Hz成分的信号20Hz成分持续时间较长在时频图上表现为一条横贯中低频的亮带50Hz如果只在某个时间段出现则表现为那个时间段内的一条局部亮带。一眼就能看出频率成分随时间的变化。在Python里做CWT非常方便PyWavelets库一行代码就能搞定import numpy as np import pywt import matplotlib.pyplot as plt fs 1000 t np.linspace(0, 1, fs) # 构造一个 20Hz 持续全程 50Hz 仅中间 0.4~0.6s 出现的信号 sig np.sin(2 * np.pi * 20 * t) sig[400:600] np.sin(2 * np.pi * 50 * t[400:600]) scales np.arange(1, 128) coefs, freqs pywt.cwt(sig, scales, morl, sampling_period1/fs) plt.pcolormesh(t, freqs, np.abs(coefs), shadinggouraud) plt.ylabel(频率 (Hz)) plt.xlabel(时间 (s)) plt.show()运行这段代码你能直观看到20Hz那条亮带从头亮到尾而50Hz那条亮带只在0.4秒到0.6秒之间出现。这就是小波分析相对傅里叶变换最直观的优势时频定位能力。3. 离散小波变换与多分辨率分解实战3.1 从连续到离散为什么要离散化连续小波变换虽然直观但有个现实问题它把尺度a和平移因子b都当作连续变量处理计算量巨大而且会产生大量冗余信息。实际工程里我们几乎不会直接计算CWT而是使用离散小波变换DWT对尺度和平移都进行离散化采样最常见的是取a 2^jb k * 2^jj和k都是整数。这种“二进离散”方式既保留了多分辨率分析的核心能力又大幅压缩了计算规模。更妙的是离散小波变换可以借助滤波器组来实现这就是Mallat算法。它把分解过程拆成两步先用一个低通滤波器提取信号的“近似”成分低频轮廓再用一个高通滤波器提取“细节”成分高频变化。每一步之后都对信号做一次二抽取也就是隔一个点取一个让数据量减半。这样不断重复每一层得到上一层的近似系数和细节系数。3.2 一层分解到底做了什么滤波器组视角很多人第一次接触DWT时容易被“低通、高通、下采样”这些词绕晕我用人话拆一遍。假设你有一段长度N的离散信号它经过一对互补的滤波器低通高通输出长度仍然是N的两路信号。低通输出保留了原始信号里的低频趋势高通输出保留了高频细节。紧接着做二抽取两路信号长度都变成N/2。这时低通输出的那一路叫“近似系数”本质上是原信号的下采样模糊版本高通输出的那一路叫“细节系数”记录的是原信号在高频部分的变化信息。两层分解就是对上一层的近似系数再做一次同样的低通/高通滤波和下采样。这个过程非常像把一张图片不断缩小每缩小一次既保留缩略图近似又记录下缩小时丢掉的细节细节。一级一级往下分解就得到一棵“小波分解树”。在Python里用PyWavelets完成三层分解的代码极其简单import pywt # 生成一段测试信号 t np.linspace(0, 1, 1000, endpointFalse) sig np.sin(2 * np.pi * 10 * t) 0.5 * np.sin(2 * np.pi * 50 * t) np.random.randn(1000) * 0.2 # 三层离散小波分解 coeffs pywt.wavedec(sig, db4, level3) cA3, cD3, cD2, cD1 coeffs print(cA3.shape, cD3.shape, cD2.shape, cD1.shape)这里的cA3是第三层近似系数cD3、cD2、cD1分别对应第三、第二、第一层细节系数。重构时用pywt.waverec(coeffs, db4)就能还原出原始信号在无损耗处理的前提下。这也是小波去噪、小波压缩的基础操作。3.3 层数怎么定分解层数与最大可分解层数分解层数并不是越多越好。层数越多低频越“粗”但每多分解一层数据长度就减半一次能分解的最大层数和信号长度直接相关。PyWavelets提供了现成的估算函数pywt.dwt_max_level它会根据信号长度和小波滤波器长度自动算出最大允许层数。实际选层时通常要结合信号的采样率和目标频段来定。举个例子采样率1000Hz奈奎斯特频率是500Hz。第一层细节系数对应250-500Hz第二层对应125-250Hz第三层对应62.5-125Hz。如果你想分析的是10-20Hz这个频段的成分那就要分解到第四层或第五层让这个频段落在某一层的近似或细节系数里。层数选得太少目标频段的特征混在较浅层的系数里看不清楚层数选得太多不仅计算量增加边界效应也会越来越明显。4. 小波基怎么选没有最好的只有最合适的4.1 常见小波族和它们的性格差异小波基的选择是新手最容易纠结的问题。答案是没有绝对正确的小波只有最适合当前信号特征的小波。不同小波族的时频特性差异很大我总结了几个最常用的供你按需选择。小波族简称特点典型场景Haar小波haar结构最简单类似方波紧支撑、正交但不连续教学示例、快速原型验证Daubechies小波dbN正交、紧支撑N越大消失矩越高、越平滑计算量越大通用信号分析、去噪、压缩Symlets小波symNdb小波的改进版近似对称相位畸变更小对相位敏感的生物医学信号Coiflets小波coifN比db更对称消失矩更高图像处理、需要对标更对称的场景Biorthogonal小波biorNr.Nd可对称、可精确重构但正交性被放弃图像压缩、信号重构如果你拿不准用哪个我的惯用做法是先用db4或sym4做一个快速分解观察细节系数里是否清晰分离出目标特征如果效果不理想再逐个换小波族对比。不需要一上来就迷信“最强小波”大多数工程场景下db族和sym族已经够用。4.2 消失矩与紧支撑的取舍选择小波时有两个关键指标需要理解消失矩和紧支撑性。消失矩决定了小波对多项式信号的“钝感力”。一个具有N阶消失矩的小波对N-1阶以下的多项式信号做内积时结果恒为0。这个性质在实际应用中非常重要因为它意味着小波系数只反映信号的“非多项式”成分也就是真正的高频突变或细节变化。消失矩越高对平滑信号的抑制能力越强越能突出奇异点但代价是小波波形变得越长时间定位能力变差边界效应也更明显。紧支撑性则决定了小波在时域上“只占多大范围”。支撑越短小波越“局部”对瞬时突变的时间定位越准适合检测脉冲类信号支撑越长小波越“发散”频域选择性越好但时间定位变模糊。Haar小波支撑最短但消失矩只有1阶频域表现一般db8支撑较长消失矩达到8阶频域分辨更细但对瞬态事件的定位就不如Haar灵敏。所以选小波本质上是做一次“时间分辨率vs频率分辨率”的平衡跟STFT选窗长的逻辑类似只是小波把它变成了可以随尺度变化的可调机制。4.3 边界效应处理信号两端时最容易踩的坑所有基于滤波器组的变换都躲不开边界效应离散小波变换也不例外。信号两端在进行滤波时滤波器窗口会“伸到”信号范围之外这时候必须决定“信号外面到底是什么”。PyWavelets提供了多种扩展模式以下是几种常见的zero补零。最简单但两端会引入突变通常只用于没有边界要求的场景。symmetric镜像对称对大多数自然信号比较友好视觉效果平稳是默认常用选项。reflect反射类似symmetric但镜像方式不同。periodization周期延拓保证DWT分解后各层系数长度相等是做多层分解和重构时最省心的选择。我自己的经验是如果信号本身是周期性的比如旋转机械的振动信号首选periodization如果信号两端走势平稳可以用symmetric。但一定不要默认用补零否则分解结果两端会多出明显的高频伪影并且随着层数增加伪影范围会逐渐向中间扩散。5. 小波去噪实操从原理到代码一次走通5.1 去噪原理小波系数为什么能区分噪声和有效信号小波去噪是目前工业界用得最多的方向之一它背后的逻辑其实很朴素有效信号在小波域里的系数通常幅值较大、分布稀疏而高斯白噪声在小波域里依然是高斯白噪声系数幅值普遍较小。所以只要找到合适的阈值把小于阈值的系数视为噪声置零或收缩再把剩余系数重构回去就能得到去噪后的信号。这里最经典的估计方法是Donoho提出的VisuShrink阈值阈值T sigma * sqrt(2 * ln(N))其中N是信号长度sigma是噪声标准差。sigma通常用第一层细节系数的中位绝对偏差来估计公式是sigma median(|cD1|) / 0.6745。这个0.6745来自正态分布的性质它能保证估计出来的标准差不受信号本身幅值影响只反映噪声水平。阈值选定后有两种常见的处理方式硬阈值系数绝对值小于等于阈值直接置零大于阈值保持不变。优点是保留细节缺点是重构信号容易出现局部抖动和不连续。软阈值系数绝对值小于阈值置零大于阈值的向零点收缩一个阈值量。优点是重构更平滑缺点是会把真实信号里较大的系数也“削”掉一点。实际工程里如果信号是振动冲击信号想保留冲击的尖锐性我倾向于用硬阈值如果是平滑的生物电信号用软阈值效果更自然。5.2 完整去噪流程与代码下面我用一个仿真信号演示完整去噪流程一段由低频正弦波叠加瞬态冲击的信号再人为加入高斯白噪声然后通过“分解-阈值-重构”三步完成去噪。import numpy as np import pywt # 1. 生成仿真信号低频正弦 冲击成分 噪声 fs 1000 t np.linspace(0, 2, 2 * fs, endpointFalse) clean np.sin(2 * np.pi * 5 * t) impulse np.zeros_like(t) impulse[300:330] 0.8 impulse[900:930] -0.6 impulse[1500:1530] 0.5 clean clean impulse np.random.seed(42) noisy clean 0.15 * np.random.randn(len(t)) # 2. 小波分解 wavelet db4 level 4 coeffs pywt.wavedec(noisy, wavelet, modesymmetric, levellevel) cA4, cD4, cD3, cD2, cD1 coeffs # 3. 噪声标准差估计与阈值计算 sigma np.median(np.abs(cD1)) / 0.6745 thresh sigma * np.sqrt(2 * np.log(len(noisy))) print(估计噪声标准差:, sigma, 阈值:, thresh) # 4. 对每一层细节系数做软阈值处理 coeffs_thresh [cA4] for i in range(1, len(coeffs)): coeffs_thresh.append(pywt.threshold(coeffs[i], thresh, modesoft)) # 5. 重构去噪信号 denoised pywt.waverec(coeffs_thresh, wavelet, modesymmetric) # 6. 查看去噪前后的信噪比 def snr(orig, est): noise orig - est return 10 * np.log10(np.sum(orig**2) / np.sum(noise**2)) print(原始含噪信号信噪比: %.2f dB % snr(clean, noisy)) print(去噪后信号信噪比: %.2f dB % snr(clean, denoised))从输出结果可以看到去噪后的信噪比通常会提高10dB以上。如果你还想更精细可以对不同层的细节系数分别设置不同的阈值甚至用每层单独估计的sigma来计算各自的阈值。这样做的原因是不同层的噪声能量分布并不完全均匀特别是经过滤波器后较低层系数里的噪声占比往往更高。5.3 阈值去噪的局限性阈值去噪也不是万能的。当信号本身的细节成分的幅值和噪声水平接近时无论阈值怎么选都会在“误杀真实信号”和“放任噪声残留”之间摇摆。另外阈值法对非高斯噪声的抑制效果有限比如脉冲噪声、有色噪声这种情况下需要考虑更复杂的贝叶斯估计或基于模型的方法。另一个容易被忽略的问题是阈值去噪前必须先把信号的直流分量或趋势项去掉。否则大的趋势项会把分解后的近似系数撑得很大阈值估算时容易被带偏。实际工业数据往往有缓慢的基线漂移我一般会先对数据做一次均值移除或高通预滤波再进入小波去噪流程。6. 小波分析常见问题与排查技巧6.1 时频图上频率轴对不上怎么办用pywt.cwt做连续小波变换时如果直接用scales数组作为纵轴看到的图像会以尺度为单位。要把它换算成频率就必须提供两个参数母小波的中心频率f0和采样周期sampling_period。PyWavelets内部已经有预设的中心频率表直接指定sampling_period后pywt.cwt的返回结果里会自动带上换算后的频率数组freqs。很多人在网上找了一堆代码直接抄了coeffs没接freqs最后画图时用尺度轴当频率轴整张图的物理含义全错了。如果你用的是其他库或自己写的CWT经验公式是对应频率约等于母小波中心频率乘以采样率再除以尺度。每次做图前检查一下最大频率是否接近奈奎斯特频率是验证换算是否正确的最快方法。6.2 信号长度不是2的整数次幂很多教材里的DWT示例都喜欢用长度为2的整数次幂的信号因为这样每一层分解都能整除以2系数长度一目了然。但真实数据哪有这么听话PyWavelets的做法是在分解时根据选择的边界模式自动处理长度不整除的情况。例如在symmetric模式下滤波器输出后下采样时遇到奇数长度会通过延拓来补全下一层系数长度可能不完全等于上一层的一半。如果你发现wavedec返回的系数长度跟你预期的不一致不用慌只要分解和重构使用的mode一致waverec就能正确还原信号。多层分解遇到奇数长度时我的建议是用modeperiodization它能保证各层系数长度的整齐关系后续处理逻辑会更清晰。6.3 实时处理中如何避免“边界污染”小波分解是基于整段信号的传统做法在实时流式数据处理里并不友好。每来一个新样本就重新对整段信号做一次DWT计算量浪费严重不说边界处的系数还会不断变化导致结果不稳定。一种常见的工程做法是“分块处理”把信号按固定长度比如4096点切片每片独立做小波分解和重构但块与块之间需要留出重叠区域避免边界效应扩散到有效区域。另一种做法是用滑窗增量更新的方式但对初学者来说实现门槛较高。我建议先用分块重叠的方式跑通实时链路后续再逐步优化延迟和计算量。保留上一块的最后一层近似系数作为下一块的初始状态是一种可行的折中方案具体实现细节要看你的实时性要求。6.4 常见问题速查表问题现象可能原因排查方法时频图全是亮点、看不出结构纵轴直接用了尺度而非频率确认传入sampling_period检查freqs数组范围重构信号两端严重失真边界模式选择不当更换为symmetric或periodization去噪后信号反而变“软”了软阈值削掉了真实冲击幅值改硬阈值或保留较大系数分解层数报错数据长度不满足该小波在该层数的要求缩短层数或对信号做补零延拓不同小波分解结果差异巨大小波基和信号特征不匹配对比db4、sym5、bior3.5等族的分解结果CWT计算特别慢尺度数量取得过多缩小scales范围或用对数间距覆盖目标频段7. 从入门到应用小波分析的下一步如果这篇你已经看下来并且亲手跑过一遍CWT时频图、DWT分解和阈值去噪那么小波分析的地基就算是打好了。接下来可以根据你的应用场景往两个方向深入。一个是往“特征提取”方向走。比如在故障诊断里用DWT把振动信号分解成多个频带系数再从每个频带里计算能量占比、均方根值、峭度等统计特征喂给分类器做故障识别。这套流程在工业领域非常成熟关键点在于怎么根据故障特征频率选择分解层数和小波基让故障信息尽可能集中到某几个频带里。另一个是往“小波包分析”方向走。DWT只对低频部分做逐层分解高频细节不再细分对某些高频信息丰富的信号会丢失不少信息。小波包分析把高频部分也继续分解得到更完整的频带划分对信号压缩、图像处理、语音增强等场景更有优势。它的思路和DWT几乎一致只是多了一棵“完整二叉树”。我个人在实际操作中体会最深的一点是小波分析最重要的不是会调库、会画图而是理解尺度、频率、层数、边界模式这几个参数之间的联动关系。参数调得好不好直接决定你是从信号里“看到”了特征还是“脑补”出了特征。入门阶段多拿几组典型信号跑一跑用自己熟悉的场景去验证输出结果是否合理比背一百条理论公式更管用。希望这篇基础介绍能帮你把这条路顺利走通。
阅读完成 · 觉得有帮助？
