简介：一个基于MATLAB的谐波小波滤波算法脚本，专注于从复杂信号中提取特定频率成分并完成滤波，适用于电力系统谐波检测、结构振动分析及声学信号处理等场景。小波分析具备时频局部化特性，相比传统傅里叶变换能更精确定位信号的频率与时间信息；脚本围绕谐波小波这一分支展开，涵盖数据预处理、连续/离散小波变换、频率选择、阈值滤波、逆变换重构与结果评估等完整流程，并支持通过调整尺度或频率等参数优化滤波效果。压缩包仅含1个m文件，大小524B，轻量简洁，适合有一定MATLAB基础、需要快速上手小波滤波或特定频率提取方法的工程师与研究人员参考，也可作为高校信号处理课程的小型教学示例。目前已有470人学习下载，通过运行脚本可直接观察特定谐波频率的提取效果，便于对照理论学习或基于自身信号数据进行二次开发。
1. 谐波小波为什么比“滤波+傅里叶”更适合特定频率提取
电网谐波分析、轴承故障诊断、声学共振检测这类场景里，最常见的需求是从一串含噪信号里把某个频率分量完整抠出来，同时保留它的幅值包络随时间的变化。传统做法是带通滤波器加傅里叶变换，但带通滤波器有两个硬伤：中心频率和带宽互相牵制，频率漂移时系数要重新设计；而傅里叶变换只能告诉你信号里有这个频率，给不出它出现在哪段时间、强度怎么波动。谐波小波滤波直接绕开了这两点——它在频域里构造一个“盒式”窗口，把目标频段内的分量原样保留，频段外全部置零，再用逆变换回到时域。这意味着“特定频率提取”被转换成“选一段尺度范围做系数掩码”，窗口位置和宽度是解耦的，频率漂移时不用重新设计滤波器，只需要平移尺度区间。harmonic_3.rar 里的 harmonic_3.m 就是这套思路的 MATLAB 最小实现，适合正在做谐波分离、振动特征提取或者小波课程设计的读者。
2. 谐波小波的盒式频谱与频率选择原理：从尺度到频段的映射关系
2.1 谐波小波的频域构造：为何频带内全通、频带外全阻
谐波小波（Harmonic Wavelet）是 Newland 在 1993 年前后提出的一类复值小波，它的特点是频域响应函数是理想盒状。以最基础的实值谐波小波为例，其傅里叶变换在频域满足：
当 2π ≤ |ω| ≤ 4π 时，Ψ̂(ω) = 1，其余为 0也就是说，它在频域里是一段“全通、陡截止”的窗口，没有旁瓣、没有过渡带。这个性质是传统 FIR/IIR 带通滤波器做不到的——后者无论如何设计，过渡带和旁瓣衰减之间都存在折中。谐波小波把这种折中转移到小波系数域：时域滤波器的阶数对应到小波分析中的尺度层数，层数越多，频域分得越细，但时域支撑越长，边缘效应越明显。
用一个具体例子来说明。假设信号的采样率 fs = 1000 Hz，目标频率是 50 Hz 的基波和 150 Hz 的 3 次谐波。传统做法是两个带通滤波器分别提取，再分别做希尔伯特变换求包络。谐波小波的做法不同：把信号变换到小波系数矩阵，每一行系数对应一个频段，找到 50 Hz 和 150 Hz 各自所在的行，把其他行系数置零，最后逆变换。整个过程只操作系数矩阵，不碰滤波器系数，频率选择变成“查表定位”。
2.2 尺度、中心频率与带宽：把“特定频率”翻译成小波参数
要操作小波系数，核心是把目标频率换算成尺度（scale）或伪频率（pseudo-frequency）。MATLAB 的 cwtft 使用对数尺度序列，每个尺度 s 对应的伪频率由小波中心频率 fc 决定：
f = fc / (s * Δt)其中 Δt = 1/fs 是采样间隔，fc 是所选小波的中心频率（Morlet 小波默认 fc = ω0 / (2π) ≈ 0.9549，ω0 = 6）。反过来，给定目标频率 f_target，需要的尺度是 s = fc / (f_target * Δt)。
harmonic_3.m 这类脚本里，最关键的参数就是尺度序列的选择。在 MATLAB 中我一般用 1/fs 到 1 的对数间隔序列，覆盖从 Nyquist 频率到最低分析频率的全范围。确定目标频段 [f1, f2] 后，对应的尺度区间是 [s1, s2] = [fc/(f2Δt), fc/(f1Δt)]。注意尺度与频率是反比关系，频率越高尺度越小。
下表给出 fs = 1000 Hz、Morlet 小波（fc ≈ 0.9549）时几个典型频段对应的尺度范围，方便核对计算是否正确：
| 目标频段 (Hz) | 中心频率 (Hz) | 尺度下限 s_min | 尺度上限 s_max | 尺度区间跨度 |
|---|---|---|---|---|
| 45 - 55 | 50 | 17.36 | 21.22 | 3.86 |
| 145 - 155 | 150 | 6.16 | 6.58 | 0.42 |
| 295 - 305 | 300 | 3.13 | 3.24 | 0.11 |
| 495 - 505 | 500 | 1.89 | 1.93 | 0.04 |
可以看到，频率越高，对应尺度区间越窄。这就是谐波小波提取特定频率的第一个坑：低频段尺度间隔大，掩码容易做；高频段尺度间隔小，对数尺度序列的分辨率可能不够，导致掩码把旁边的频率也圈进来。
2.3 与 Morlet、db 小波、小波包的边界与取舍
谐波小波和 Morlet 小波、Daubechies 小波、小波包经常被放在一起比较，但它们解决的不是同一个问题。
Morlet 小波是复数小波，频域响应是高斯形状而非盒状，适合时频分析和瞬时频率提取，但因为频带边缘缓慢衰减，做“特定频率提取”时会产生泄漏——目标频段附近的噪声会被部分保留；Daubechies 小波是正交小波，重构完美，但频域是振荡衰减的，频带边界不锐利，而且离散小波变换的频带是二进划分的，只能提取类似 [fs/4, fs/2] 这样的固定频段，无法任意指定 45-55 Hz；小波包解决了二进频带过粗的问题，可以细分到任意频段，但节点选择逻辑复杂，需要先做完整的小波包分解再按频率排列节点。
谐波小波的优势在于它的盒状频域响应正好服务于“提取某一段频率”这个需求。harmonic_3.m 里可以同时使用 cwtft（连续小波变换）和 wavedec（离散小波变换），前者用于精确定位频段，后者用于验证滤波效果。两者的取舍原则是：需要保留信号的时变包络，用 cwtft；只需要频段分离和重构，用小波包或 wavedec。
3. harmonic_3.m 核心流程重构：cwtft 频带掩码与系数重建
3.1 信号结构与预处理
harmonic_3.m 面对的信号通常是采样率固定的离散序列，包含基波、若干次谐波和噪声。脚本的第一步是确定信号长度和采样率，然后构造时间轴。常见做法是直接在脚本开头用load或sin合成测试信号。无论数据来源是什么，预处理都要做两件事：去直流分量（避免小波变换后低频系数被直流淹没）和信号长度对齐到 2 的幂次（便于 wavedec 的层数计算，不过 cwtft 不要求）。
% harmonic_3.m 重构示例 - 预处理部分 fs = 1000; % 采样率 1000 Hz t = 0:1/fs:1-1/fs; % 1 秒时长 N = length(t); % 构造信号: 50Hz 基波 + 150Hz 三次谐波 + 随机噪声 sig = 1.0 * sin(2*pi*50*t) + 0.4 * sin(2*pi*150*t + pi/4); sig = sig + 0.1 * randn(size(t)); % 添加高斯白噪声 sig = sig - mean(sig); % 去除直流分量这段代码里fs直接决定所有频率映射关系，改采样率必须同步调整后续所有尺度参数；去直流用mean而非高通滤波器，是为了避免滤波器的相位延迟影响后续小波变换的时间定位精度。
3.2 使用 cwtft 做连续小波变换：参数与解释
harmonic_3.m 里如果要精确提取 50 Hz 和 150 Hz，我会用cwtft而不是新版cwt函数。原因是cwtft直接返回小波系数矩阵和对应的尺度序列、伪频率序列，可以方便地对系数做任意掩码操作；而新版cwt返回的是 CWT 对象，重建函数icwt的掩码操作接口更受限。
% harmonic_3.m 重构示例 - 小波变换与频率定位 scales = 1:100; % 尺度序列, 覆盖频率范围约 9.5 Hz ~ 954.9 Hz cwtS = cwtft({sig, 1/fs}, 'scales', scales, 'wavelet', 'morl'); freqs = centfrq('morl') ./ (scales * 1/fs); % 计算伪频率 % 定位目标频段对应的尺度索引 f_target = 50; % 目标频率 50 Hz [idx_min, ~] = min(abs(freqs - f_target)); scale_idx = find(abs(freqs - f_target) == idx_min);cwtft的第一个输入是元胞数组{sig, 1/fs}，这是它的固定调用格式，表示“信号 + 采样间隔”；'wavelet', 'morl'指定 Morlet 小波，centfrq('morl')返回 Morlet 的中心频率约为 0.9549。尺度序列选 1:100 是保守做法，实际中可以用logspace(log10(1), log10(N/2), 200)生成 200 个对数间隔的尺度，频段定位更精细。
这段代码的边界条件需要注意：freqs的最大值是centfrq('morl') * fs（尺度为 1 时），最小值受最大尺度限制。如果目标频率超出这个范围，min(abs(freqs - f_target))会定位到端点尺度，提取结果是错误的。
3.3 频带掩码与逆变换重建
获得系数矩阵后，特定频率提取就是一次矩阵行操作。找到目标频率对应的尺度索引，保留该索引附近的若干行系数，其余行置零，然后调用icwtft做逆变换。
% harmonic_3.m 重构示例 - 频带掩码与重建 cfs = cwtS.cfs; % 小波系数矩阵 (尺度 x 时间) bandwidth = 5; % 频带半宽, 单位: 尺度索引 mask = zeros(size(cfs)); mask(scale_idx-bandwidth:scale_idx+bandwidth, :) = 1; cfs_masked = cfs .* mask; % 应用掩码 % 重建时域信号 cwtS_recon = cwtS; cwtS_recon.cfs = cfs_masked; sig_filtered = icwtft(cwtS_recon); % 去除虚部的数值残留 sig_filtered = real(sig_filtered);掩码的bandwidth参数直接控制提取频段的宽度。等价地，在频率域中，这相当于保留以 50 Hz 为中心、宽度约为2 * bandwidth * Δf的频段（Δf 为相邻尺度的频率间隔）。实际信号频率存在波动时（比如电网频率在 49.8-50.2 Hz 之间变化），把bandwidth设大一点是有效的；但设得太大，旁边的 150 Hz 或其他干扰会被卷入。我通常先用 3 个小尺度索引做初步提取，观察重建信号频谱中相邻频率的衰减情况，再决定是否加宽。
icwtft的重建质量取决于两个因素：尺度序列是否覆盖了信号的主要能量频段，以及掩码是否破坏了小波系数矩阵的完整性。如果尺度序列最小时频率低于目标频率，重建信号的幅值会偏小，因为信号能量有一部分落在未覆盖的尺度上。
3.4 wavedec 二进频带与小波包节点的补充说明
cwtft适合精细提取，但计算量大。harmonic_3.m 若追求速度，可以用wavedec配合小波包。离散小波变换把信号分解成近似系数和细节系数，每一层细节对应一个二进频带。以 fs = 1000 Hz 为例，第 1 层细节对应 250-500 Hz，第 2 层对应 125-250 Hz，第 3 层对应 62.5-125 Hz。要提取 50 Hz，落在第 4 层细节（31.25-62.5 Hz）。
% 使用 wavedec 做频段逼近 wname = 'db8'; level = 4; [C, L] = wavedec(sig, level, wname); % 提取第4层细节系数 (31.25-62.5 Hz) detail4 = wrcoef('d', C, L, wname, level);db8是 Daubechies 小波系中支撑较长的基函数，频域局部性好，适合谐波分离；wrcoef从分解系数中重构指定层细节，不改变其他层系数。注意这里的频段是 31.25-62.5 Hz，比 50 Hz 单频宽得多，如果信号里还有 40 Hz 的干扰，wavedec会一并保留，所以它只能作为粗糙分离手段，精细提取必须回到cwtft。
4. 参数标定与排错：尺度间隔、边缘效应、阈值与采样率之间的权衡
4.1 尺度序列的对数间隔 vs 线性间隔
连续小波变换中尺度序列的选择直接决定频率分辨率。线性尺度 1:100 对低频段（大尺度）分辨率差，对高频段（小尺度）分辨率好；对数尺度则让每个频段的相对分辨率一致，适合宽频分析。harmonic_3.m 里如果目标频率是 50 Hz 这种低频，用对数尺度更能保证 45-55 Hz 的频段内至少有 5 个以上尺度索引可用，否则掩码圆滑度不够，重建信号会出现波纹。
% 对数尺度序列生成 num_scales = 200; scales_log = logspace(log10(2), log10(N/2), num_scales); freqs_log = centfrq('morl') ./ (scales_log * 1/fs);经验值是让目标频段内至少包含 8-12 个尺度点。检查方法是计算目标频段上下限对应尺度索引之差，差小于 5 就增加num_scales。高频段（比如 500 Hz）在 200 个尺度下可能只有 2-3 个点，这时候要么增加尺度数到 400，要么改用小波包做该频段的提取。
4.2 边缘效应与采样率对频段上限的约束
小波变换在信号两端会产生边缘畸变，因为小波在边界处没有足够的信号支撑。cwtft默认使用周期性延拓，谐波信号恰好首尾不连续时，边缘效应会非常明显，表现为重建信号两头出现异常波动。排错方法如下：
% 观察边缘效应: 对比原信号与重建信号的端点差异 figure; plot(t(1:50), sig(1:50), 'b'); hold on; plot(t(1:50), sig_filtered(1:50), 'r'); legend('原始信号', '提取后信号');如果端点误差超过信号幅值的 5%，就要考虑剔除边缘。常见做法是把有用信号段裁剪掉前 10% 和后 10%，或者用wkeep保留中间部分。另外，采样率决定了可提取频率的上限。根据 Nyquist 定理，cwtft能分析的频率上限是 fs/2，但执行掩码操作时尺度为 1 的系数对应频率接近 fs/2，此时小波只有几个采样点支撑，时域分辨率极差，提取出的波形畸变严重。所以实际可用的频率上限建议不要超过 fs/5，谐波提取时我会把注意力放在 fs/10 以下。
4.3 阈值选择与小波系数收缩
harmonic_3.m 中滤波操作如果包含降噪，会涉及系数阈值处理。与小波硬阈值相比，软阈值在去除低幅噪声时更平滑，但会压缩信号幅值。谐波提取场景下目标频率幅值通常远大于噪声，先用硬阈值剔除小系数，再对保留的系数做掩码，效果更直接。
% 对每个尺度计算阈值 (基于噪声标准差) sigma = median(abs(cfs(:))) / 0.6745; thr = sigma * sqrt(2 * log(N)); cfs_wth = wthresh(cfs, 's', thr); % 软阈值 cfs_wth = cfs_wth .* mask; % 再应用频带掩码参数说明：median(abs(cfs(:))) / 0.6745是稳健的噪声标准差估计，来源于 Donoho 的小波降噪理论；0.6745 是标准正态分布的中位数与标准差之比；thr是通用阈值，随信号长度 N 增大而增大。软阈值对重建信号的幅值有压缩作用，如果后续要做幅值定量分析（比如谐波幅值监测），建议改用'h'硬阈值。
4.4 与常见误用方式的对比：滑动窗口、低通滤波和阈值降噪
谐波小波滤波经常被拿来与滑动窗口滤波、卡尔曼滤波和普通低通滤波对比，容易踩的坑是选错工具。滑动平均窗在抑制白噪声时表现尚可，但它的频域响应是缓慢振荡衰减的旁瓣，对邻近频率的抑制能力不足；卡尔曼滤波需要建立信号的状态空间模型，谐波成分如果多于两个，状态维度和调参工作量迅速膨胀，不适合脚本化的快速滤波。
下表总结了不同滤波手段在“提取 50 Hz 谐波、抑制 60 Hz 邻频干扰”场景下的表现：
| 方法 | 频带形状 | 邻频抑制 | 相位失真 | 适用场景 |
|---|---|---|---|---|
| 二阶低通滤波 | 缓慢衰减 | 差 | 有 | 去高频噪声 |
| 滑动平均窗 | 旁瓣振荡 | 中 | 有 | 时域平滑 |
| 谐波小波掩码 | 盒状陡截止 | 好 | 小 | 特定频率提取 |
| 小波软阈值 | 依赖小波基 | 中 | 无 | 信号降噪 |
谐波小波的优势不是“更先进”，而是它的频率响应形状与“提取某一频率”的需求恰好匹配。对应到 harmonic_3.m 的调试思路上，如果提取结果混入邻频干扰，首先检查掩码宽度和尺度序列密度；如果提取结果幅值偏小，检查尺度序列是否完整覆盖目标频段；如果波形两端畸变，检查边缘效应。
5. 谐波小波提取效果的验证：合成信号测试与频谱残差分析
5.1 构造带谐波和噪声的基准信号
验证谐波小波滤波是否成功，不能用真实数据直接下结论，因为真实信号的“真值”未知。正确做法是构造一个频率、幅值、相位完全已知的合成信号，让谐波小波提取后再对比。
% 构造已知真值的验证信号 fs = 1000; t = 0:1/fs:2-1/fs; N = length(t); f1 = 50; f3 = 150; f5 = 250; % 基波、3次、5次谐波 A1 = 1.0; A3 = 0.4; A5 = 0.2; sig_truth = A1*sin(2*pi*f1*t) + A3*sin(2*pi*f3*t + pi/4) ... + A5*sin(2*pi*f5*t + pi/3); sig_noise = sig_truth + 0.15 * randn(size(t)); % 用谐波小波提取 50 Hz 分量(沿用第3章的 cwtft 流程) sig_50hz = extract_harmonic(sig_noise, fs, 50, 3); % 封装函数extract_harmonic内部的带宽参数半宽取 3 个尺度索引，对于对数尺度序列在 50 Hz 附近大约对应 ±2 Hz 的频段。提取后计算三个指标：信噪比提升（提取信号与真值的相关系数）、幅值误差（提取信号包络的均值与 A1 的差）、频谱残差（170 Hz 处是否有伪峰值）。
5.2 频谱泄漏检查与残余谐波量化
提取效果的量化核心是看目标频率之外还有多少残余能量。将提取前后的信号分别做功率谱分析，对比目标频段外的能量占比，这叫频谱泄漏检查。如果掩码宽度过大，150 Hz 和 250 Hz 的谐波能量会泄漏进 50 Hz 提取结果中，表现为重建信号频谱在偏离 50 Hz 的位置出现异常凸起。
% 频谱对比与泄漏量化 [pxx_orig, f_axis] = pwelch(sig_noise, hann(500), 250, 512, fs); [pxx_filt, ~] = pwelch(sig_50hz, hann(500), 250, 512, fs); % 计算45-55Hz内的能量占比 band_idx = f_axis >= 45 & f_axis <= 55; orig_band_power = sum(pxx_orig(band_idx)) / sum(pxx_orig); filt_band_power = sum(pxx_filt(band_idx)) / sum(pxx_filt); % 计算残余谐波比 (150Hz处能量与50Hz处能量的比值) ratio_orig = pxx_orig(f_axis >= 145 & f_axis <= 155) ... / pxx_orig(band_idx); ratio_filt = pxx_filt(f_axis >= 145 & f_axis <= 155) ... / pxx_filt(band_idx);pwelch的窗口长度取 500 个样本，频率分辨率约为 1 Hz，足以区分 50 Hz 和 45 Hz 的泄漏成分。band_power从原始信号里的约 0.25 提升到提取后的 0.9 以上，说明 50 Hz 分量保留完整；ratio_filt相对于ratio_orig下降至少两个数量级，说明邻频谐波抑制到位。如果ratio_filt只下降一个数量级，多半是掩码宽度太大，把 150 Hz 的部分能量圈进了 50 Hz 频段。
5.3 频率漂移时的动态频带跟踪技巧
实际工程信号里的谐波频率很少是理想稳定的。电网频率会在 49.8-50.2 Hz 之间摆动，机械振动信号的转频也会缓慢变化。固定掩码在这种情况下会丢失部分信号能量导致幅值波动。常见做法是先对信号做一次功率谱估计确定实际频率中心，再把掩码中心对准实测值，而不是直接使用理论频率。
一个轻量级的实现是先用findpeaks在功率谱中找到目标频率附近的最大峰，以其位置重新计算尺度索引：
% 动态中心频率检测 [pk_mag, pk_loc] = findpeaks(pxx_orig, 'MinPeakHeight', ... 0.1*max(pxx_orig), 'SortStr', 'descend'); f_detect = f_axis(pk_loc(1)); % 主峰频率 % 以检测到的频率为中心重新计算尺度掩码 scale_center = round(centfrq('morl') / (f_detect * (1/fs))); mask = zeros(size(cfs)); mask(max(1, scale_center-4):scale_center+4, :) = 1;MinPeakHeight设为最大峰值的 10%，避免把噪声峰当成谐波。这种方法把掩码中心动态锁定到实际频率上，带宽只负责覆盖漂移范围，提取信号的幅值波动可以控制在 2% 以内。harmonic_3.m 在落地到实际项目时，把第 3 章的固定掩码替换为这个动态版本，滤波的适应能力会有明显提升。
本文还有配套的精品资源，点击获取
