MATLAB小波高通与低通滤波实战:从原理到代码实现
1. 从一个实际信号处理需求说起如果你手头有一段采样率不低的振动信号、一段音频或者某类传感器采集的时序数据想把它拆成快变部分和慢变部分分别研究那小波分析里的高通滤波器和低通滤波器基本是绕不开的一对工具。我最早接触这对概念是在处理一组机械振动数据的时候——低频段藏着设备的整体运行趋势高频段则对应冲击、摩擦这类瞬态事件两者混在一起看什么都看不清必须先把频带拆开。这篇内容就是围绕用 MATLAB 做小波高通/低通滤波这件事展开的。它适合三类人一是刚学完小波变换理论、但不知道怎么落到代码上的学生二是做信号处理、故障诊断、音频分析相关工作的工程师三是想用 MATLAB 快速验证滤波效果、不想从零推导公式的实践派。我会把滤波器组的设计思路、MATLAB 里几套常用实现路径、参数怎么定、结果怎么验证以及我自己踩过的坑尽量讲透。需要先明确一点小波滤波和传统的 FIR/IIR 滤波不是一回事。传统滤波器是在频域上切一刀小波滤波是在**尺度scale**上做分解把信号拆成近似分量approximation和细节分量detail近似分量对应低频细节分量对应高频。所谓小波低通滤波器和小波高通滤波器本质就是这对分解滤波器组。理解这一点后面的代码才不会写歪。2. 小波高通/低通滤波器的核心原理拆解2.1 为什么小波能同时做低通和高通传统傅里叶变换把信号拆成正弦波的叠加频域分辨率高但完全没有时间信息。小波变换换了个思路用一个可以伸缩、平移的小波基函数去和信号做内积。尺度大的时候小波被拉宽对应低频尺度小的时候小波被压缩对应高频。这个伸缩的过程天然就把信号分到了不同频带。在离散小波变换DWT里这个分解是通过一对滤波器实现的低通滤波器Low-pass Filter保留信号的低频成分输出近似系数Approximation Coefficients记作 cA。高通滤波器High-pass Filter保留信号的高频成分输出细节系数Detail Coefficients记作 cD。这两个滤波器不是随便设计的它们必须满足**正交镜像滤波器QMF**关系也就是高通滤波器是低通滤波器的镜像翻转加调制。MATLAB 里常用的dbNDaubechies、symNSymlets、coifNCoiflets这些族本质上就是提供了不同长度、不同消失矩的 QMF 对。2.2 分解与重构滤波不是单向的很多人第一次用dwt会懵我拿到了 cA 和 cD然后呢这里要区分两个动作分解Decomposition信号 → 低通 → 下采样 → cA信号 → 高通 → 下采样 → cD。下采样是因为滤波后带宽减半按奈奎斯特准则可以丢掉一半样本这样总数据量不变。重构ReconstructioncA 和 cD 分别上采样、经过重构滤波器、相加还原出原信号。如果重构滤波器设计得当这个过程是**完美重构Perfect Reconstruction**的也就是不丢信息。所以滤波在小波里通常有两种含义一种是只保留 cA 或只保留 cD 再重构得到纯低频或纯高频信号另一种是修改某些系数比如阈值去噪再重构。前者就是标题里说的低通滤波和高通滤波。2.3 单层 vs 多层滤波器的级联单层 DWT 只把信号分成两个频带。如果继续对 cA 做分解就得到多层结构这叫小波包分解或多分辨率分析MRA。层数越多低频被切得越细。层数怎么定一个经验公式是最大分解层数 floor(log2(N / (L - 1)))其中 N 是信号长度L 是小波滤波器的长度。比如db4的 L8信号长度 1024那最大层数大约是 floor(log2(1024/7)) ≈ 7 层。实际用的时候一般取 3~5 层就够了层数太多边界效应会很明显。3. MATLAB 里实现小波滤波的几条路径3.1 路径一dwt / idwt 手动控制这是最裸的方式适合你想完全掌控分解和重构过程。核心函数就两个dwt做单层分解idwt做单层重构。% 单层分解 [cA, cD] dwt(signal, db4); % 只保留低频低通滤波 signal_low idwt(cA, zeros(size(cD)), db4); % 只保留高频高通滤波 signal_high idwt(zeros(size(cA)), cD, db4);注意idwt的长度问题重构出来的信号长度可能比原信号多几个点因为滤波器有延迟。通常用idwt(cA, cD, db4, length(signal))指定输出长度或者事后用wkeep裁剪。3.2 路径二wavedec / waverec 多层处理多层分解用wavedec它返回一个系数向量 C 和一个记录长度的向量 L。% 3层分解 [C, L] wavedec(signal, 3, db4); % 提取第3层近似系数 cA3 appcoef(C, L, db4, 3); % 提取第1层细节系数 cD1 detcoef(C, L, 1); % 只保留近似系数重构低通 A3 wrcoef(a, C, L, db4, 3); % 只保留第1层细节重构高通 D1 wrcoef(d, C, L, db4, 1);wrcoef这个函数特别好用它可以直接把某一层的近似或细节单独重构回原信号长度不用你手动拼系数。做低通滤波就取a做高通滤波就取d。3.3 路径三wden 自动去噪式滤波如果你做滤波的目的是去噪wden是更省事的选择。它内部自动完成分解、阈值处理、重构。% 启发式阈值去噪 signal_denoised wden(signal, heursure, s, one, 3, db4);这里的参数含义heursure是阈值规则s是软阈值one表示阈值不随层数缩放3是分解层数db4是小波基。软阈值会把系数往零收缩重构出来的信号更平滑硬阈值保留大于阈值的系数原值保边能力更强但可能有振铃。3.4 三条路径怎么选路径适用场景控制粒度代码量dwt/idwt单层、教学、快速验证高少wavedec/waverec多层、需要分层分析中中wden去噪为主、不想调参低极少我个人的习惯是探索阶段用dwt快速看效果正式处理用wavedecwrcoef因为分层重构对分析各频带贡献特别直观。4. 完整实操从信号生成到滤波验证4.1 构造一个可验证的测试信号光看代码没感觉得有个已知成分的信号来验证。我构造一个包含低频趋势和高频瞬态的合成信号fs 1000; % 采样率 1000 Hz t 0:1/fs:1-1/fs; % 1秒数据 f_low 5; % 低频 5 Hz f_high 120; % 高频 120 Hz % 低频正弦 高频正弦 一个冲击 sig_low sin(2*pi*f_low*t); sig_high 0.5*sin(2*pi*f_high*t); impulse zeros(size(t)); impulse(500) 2; % 第500点加冲击 signal sig_low sig_high impulse;这个信号的好处是低频、高频、瞬态三者分离得很清楚滤波后一眼就能看出哪部分被保留、哪部分被滤掉。4.2 用 wavedec 做 4 层分解wname db4; level 4; [C, L] wavedec(signal, level, wname); % 各层细节系数能量占比 for i 1:level cD detcoef(C, L, i); energy_ratio(i) sum(cD.^2) / sum(signal.^2); fprintf(第%d层细节能量占比: %.4f\n, i, energy_ratio(i)); end跑出来你会发现120 Hz 的成分主要落在第 2~3 层细节里冲击能量则分散在多个高频层。这个能量分布就是判断该保留哪层的依据。4.3 低通滤波只保留近似系数% 只保留第4层近似重构 sig_lowpass wrcoef(a, C, L, wname, level); % 对比 figure; subplot(3,1,1); plot(t, signal); title(原始信号); subplot(3,1,2); plot(t, sig_lowpass); title(低通滤波结果); subplot(3,1,3); plot(t, signal - sig_lowpass); title(被滤掉的高频部分);低通结果应该只剩 5 Hz 正弦120 Hz 和冲击都被滤掉了。第三张图就是高通残差它应该和 120 Hz 正弦加冲击高度吻合。4.4 高通滤波只保留细节系数% 把所有层的细节系数重构并相加 sig_highpass zeros(size(signal)); for i 1:level sig_highpass sig_highpass wrcoef(d, C, L, wname, i); end figure; plot(t, sig_highpass); title(高通滤波结果);这里有个细节wrcoef(d, ..., i)只重构第 i 层细节把所有层加起来才是完整的高通结果。如果你只想要某个特定频带比如只要 120 Hz 那段就单独取对应层不要全加。4.5 参数选择的计算过程小波基怎么选我一般看两个指标消失矩和支撑长度。消失矩越高对多项式类信号的逼近越好但支撑越长边界效应越明显。db4有 4 阶消失矩支撑长度 8是个折中选择。如果信号很平滑用db2或sym3就够如果信号有尖锐突变db6以上保边更好但计算量上升。分解层数怎么定前面给了公式但实际还要看你想分离的频带。假设采样率 1000 Hz想分离 100 Hz 以上的成分那第几层对应 100 Hz每层细节覆盖的频带是 [fs/2^(i1), fs/2^i]。第 2 层是 [125, 250] Hz第 3 层是 [62.5, 125] Hz。所以 120 Hz 落在第 3 层。想精确分离就得让目标频率落在某一层的中心附近。5. 常见问题与排查技巧实录5.1 重构信号长度对不上这是最高频的问题。idwt和waverec的输出长度经常比输入多几个点原因是滤波器有群延迟。解决办法有两个一是重构时指定长度waverec(C, L, wname)通常能自动对齐但idwt需要手动指定二是用wkeep裁剪首尾。注意不要用resample或简单截断来凑长度那会引入相位偏移滤波结果就不可信了。5.2 边界效应导致首尾失真DWT 在信号两端会做延拓延拓方式不对就会在首尾产生明显振荡。MATLAB 默认用对称延拓sym大多数情况够用。如果首尾失真严重可以改用per周期延拓信号本身近似周期时在信号两端各补一段数据再滤波事后裁掉降低分解层数我处理振动信号时习惯在两端各镜像延拓 100 个点滤完再裁掉边界失真能压下去一大半。5.3 高通结果里混进了低频如果你用wrcoef(d, ...)重构后发现还有低频残留八成是分解层数不够低频成分还没被完全分到近似系数里。增加层数或者检查小波基的消失矩是否足够。另一个可能是你只取了部分层的细节漏掉了某些层。5.4 软阈值和硬阈值怎么选阈值类型优点缺点适用软阈值结果平滑、无振铃幅值被压缩去噪为主硬阈值保幅值、保边可能振铃特征提取半软阈值折中需调参精细处理我一般先用软阈值看整体效果如果发现有用信号的幅值被压得太狠再换硬阈值或半软阈值。5.5 常见问题速查表现象可能原因排查方向重构长度不对滤波器群延迟指定输出长度或用 wkeep首尾振荡边界延拓不当换延拓模式或补数据高通残留低频层数不足增加分解层数低通残留高频小波基消失矩低换 db6/sym8去噪后信号变钝软阈值压缩改硬阈值或降阈值计算太慢层数或基太长降层数、换短基6. 几个容易被忽略的实操心得第一个心得先画频谱再动手。很多人上来就wavedec结果层数选错、基选错白忙活。正确顺序是先用fft看信号频谱确认目标频带在哪再反推层数和小波基。这一步花两分钟能省半小时调试。第二个心得能量占比是最好的层数指示器。前面代码里算的energy_ratio不是摆设它能告诉你信号能量主要分布在哪几层。如果某一层能量占比超过 60%那层大概率是主成分滤波时重点保护。第三个心得高通滤波不要盲目全加。很多人把 1 到 N 层细节全加起来当高通结果但如果你的目标是提取某个特定频带全加会把无关频带也带进来。正确做法是先看各层能量只加目标层。第四个心得验证滤波效果要用残差。低通结果和原信号相减得到的残差应该等于高通结果理想情况下。如果两者对不上说明重构过程有问题可能是系数被误改或长度没对齐。这个交叉验证我每次都会做。第五个心得小波基不是越高级越好。db10听起来比db4厉害但支撑长度 20边界效应和计算量都上去了。对于大多数工程信号db4到db6是甜点区。音频信号可以试试sym系列因为它的相位特性更接近线性。7. 从单层到多层的进阶玩法如果你已经把单层滤波玩熟了可以试试这几个进阶方向。小波包分解Wavelet Packet普通 DWT 只对低频继续分解高频不再动。小波包对低频和高频都继续分解频带划分更细。MATLAB 里用wpdec和wprcoef。适合需要精细频带分离的场景比如把 0~500 Hz 切成 8 个等宽子带。平稳小波变换SWT普通 DWT 有下采样导致它不是时不变的——信号平移几个点分解结果就变了。SWT 去掉下采样每个尺度输出和原信号等长时不变性好适合做特征提取和突变检测。MATLAB 里用swt和iswt。代价是数据量变大计算变慢。多分辨率分析MRA可视化用mra函数可以一次性得到各层的重构分量画出来就是经典的金字塔图低频在上、高频在下一眼看清各频带贡献。做报告和论文时特别好看。% MRA 分解 mra_sig mra(signal, 4, db4); % mra_sig 是一个矩阵每行对应一层重构分量结合阈值做自适应滤波固定阈值对不同信噪比的信号效果差异很大。可以先用一段纯噪声估计噪声水平再按阈值 sigma * sqrt(2*log(N))计算这就是经典的通用阈值公式。wden的sqtwolog规则用的就是这个。8. 一个完整的可复现脚本框架把前面的东西串起来给你一个可以直接改参数用的脚本骨架%% 参数区 fs 1000; wname db4; level 4; thr_rule heursure; thr_type s; %% 信号准备 t 0:1/fs:1-1/fs; signal sin(2*pi*5*t) 0.5*sin(2*pi*120*t); signal(500) signal(500) 2; %% 分解 [C, L] wavedec(signal, level, wname); %% 低通 sig_lp wrcoef(a, C, L, wname, level); %% 高通 sig_hp zeros(size(signal)); for i 1:level sig_hp sig_hp wrcoef(d, C, L, wname, i); end %% 去噪 sig_dn wden(signal, thr_rule, thr_type, one, level, wname); %% 验证 residual signal - sig_lp; err norm(residual - sig_hp) / norm(residual); fprintf(低通残差与高通结果相对误差: %.6f\n, err); %% 绘图 figure; subplot(4,1,1); plot(t, signal); title(原始); subplot(4,1,2); plot(t, sig_lp); title(低通); subplot(4,1,3); plot(t, sig_hp); title(高通); subplot(4,1,4); plot(t, sig_dn); title(去噪);这个脚本里那个err计算就是前面说的交叉验证。正常情况下它应该在 1e-10 量级如果大于 1e-3说明重构环节有问题回去检查系数有没有被误改、长度有没有对齐。9. 关于参数调优的一点个人体会调小波滤波参数这件事没有万能公式但有优先级。我的顺序是先定层数再定小波基最后调阈值。层数决定频带划分是最粗粒度的控制小波基决定时频局部化特性影响保边和平滑的平衡阈值只影响去噪强度是最后微调的。层数定错后面怎么调都别扭。我见过有人用 8 层分解处理 200 点的信号结果每层就剩几个系数重构出来全是边界效应。记住那个公式层数上限受信号长度和滤波器长度共同约束别硬上。小波基的选择上如果你不确定db4先跑一遍看结果再决定要不要换。sym系列在相位上更友好coif系列在消失矩和支撑长度的平衡上做得不错但计算稍慢。实际工程里db和sym占了绝大多数场景。阈值这块软阈值是默认选择除非你明确知道需要保幅值。heursure规则在信噪比未知时比较稳sqtwolog在噪声水平已知时更准。如果去噪后信号明显变钝先把软阈值换硬阈值试试还不行就手动把阈值乘个 0.5~0.8 的系数。最后说个我踩过的坑不要对已经滤波过的信号再做小波滤波。有人为了滤得更干净把低通结果再低通一次结果信号被过度平滑有用特征全没了。小波滤波一次到位效果不够就回去调参数别叠加。