Matlab连续小波变换时频分析：从原理到轴承故障诊断

简介分享一个MATLAB连续小波变换源码包面向信号处理学习者与科研人员帮助在MATLAB环境下直接实现CWT并对非平稳信号进行频谱分析。资源包含cwt1d、cwt2d、cwt3d_layer三个m文件分别覆盖一维信号、二维图像以及分层三维数据的小波变换处理代码思路清晰便于二次修改与嵌入到自己的项目中适合作为课程设计、毕业设计或论文实验的参考实现。整个压缩包共3个文件大小仅4KB轻量便携目前已有1023人浏览学习。文件内涉及的尺度参数设置、Morlet小波选择、小波系数图绘制与频谱信息提取等要点可帮助读者快速理解连续小波变换从原理到编码落地的过程节省摸索时间。1. 用连续小波变换做频谱分析先搞清楚它和 FFT 的差别在哪拿到一段振动信号、音频或者生理电信号第一反应通常是打开 matlab 直接fft()。但 FFT 给出的是一个全局平均的频谱它把所有时间点的能量折在了一起信号里某个时刻突然出现的冲击、调频成分或者短暂谐振会被平均掉甚至在频谱图上根本看不见。连续小波变换Continuous Wavelet Transform, CWT)解决的是这个问题——它在时间轴和频率轴上同时展开信号输出一张时频图既能看出哪些频率成分存在也能看出它们什么时候出现、持续时间多长。做振动频谱分析、故障诊断、脑电节律分析的人真正需要的其实是这种频率随时间变化的可视化结果。matlab 里实现 CWT 并不复杂核心函数就是一个cwt()。但很多人第一次上手就卡在三个地方小波基怎么选、尺度向量怎么设、输出的coefs矩阵和频率轴的对应关系是什么。这三个问题不解决跑出来的图要么看不出特征要么频率轴标注得完全不对。这篇文直接从理论铺到代码把从原始信号到输出频谱图的最小流程拆开讲参数逐个对齐最后把最容易翻车的边界效应和时间窗问题单独拿出来处理。适合已经会用fft()、想进一步做时频分析的人也适合做振动频谱图分析、想看懂公开代码里cwt参数的人。2. CWT 的时间-尺度原理与 matlab 里的核心函数2.1 连续小波变换的数学表达和物理含义连续小波变换的定义式是$C(a,b) \frac{1}{\sqrt{a}} \int_{-\infty}^{\infty} x(t) \psi^*(\frac{t-b}{a}) dt$其中a是尺度因子对应频率的倒数b是平移因子对应时间位置$\psi(t)$ 是小波母函数。和短时傅里叶变换STFT用固定窗长不同CWT 的窗长随尺度自动变化高频时小波被压缩时间分辨率高、频率分辨率低低频时小波被拉伸频率分辨率高、时间分辨率低。这个特性使得 CWT 特别适合处理频率跨度很大的信号比如机械振动里既有几千赫兹的冲击成分、又有几十赫兹的转频成分。注意尺度a和频率f之间存在反比关系但具体换算公式依赖小波母函数的中心频率不是简单取倒数。用 matlab 的cwt()时直接用freqs cwtfreqparams()或者输出参数里的f向量即可千万别自己手写尺度转频率。2.2 matlab 里 CWT 的三种调用方式现在 matlab 的 CWT 主要在 Wavelet Toolbox 里核心函数是cwt()。它有三种常见用法% 方式一只传信号自动选参数快速看结果 [cfs, frq] cwt(x); % 方式二指定采样频率frq 输出真实频率(Hz) [cfs, frq] cwt(x, fs); % 方式三指定小波基和频率范围做精细控制 [cfs, frq] cwt(x, fs, voicesperoctave, 16, ... wavelet, amor, FrequencyLimits, [0.5, 100]);cfs是复矩阵行对应尺度/频率列对应时间点每个元素是复数模值代表该时频点上的能量强度。frq是和行对应的频率向量画时频图时直接作为y-axis使用。参数说明fs采样频率单位 Hz决定了频率轴的上限奈奎斯特频率 fs/2。不传默认是 1此时所有频率都是归一化的看频谱分析结果时要格外小心。voicesperoctave每倍频程内的频率细分数量默认 10 或 12。调大到 16 或 32 会让时频图更平滑但计算量会增加。wavelet小波基类型。最常用的三个是amor复 Morlet适合看相位和振荡信号、morse默认通用性强、bump频率局部性好适合尖峰成分。FrequencyLimits限制分析的最低和最高频率超过范围的cfs不计算能显著加快运行速度。fs 1000; % 采样率 1000 Hz t 0:1/fs:1-1/fs; % 1 秒时间轴 % 构造一个 10Hz 正弦 250Hz 冲击的混合信号 x sin(2*pi*10*t) 0.8*sin(2*pi*250*t).*exp(-((t-0.5)/0.01).^2); [cfs, frq] cwt(x, fs, voicesperoctave, 16); % 画时频图 t_axis 0:1/fs:1-1/fs; imagesc(t_axis, frq, abs(cfs)); axis xy; % 让频率轴从低到高显示 ylim([0 500]); xlabel(时间 (s)); ylabel(频率 (Hz)); title(CWT 时频图); colorbar;2.3 小波基的选择逻辑不是越复杂越好很多刚开始做时频分析的人会纠结该选amor还是morse。我的经验是先用默认的morse根据效果再决定。Morse 小波是 matlab 从 R2016b 起的默认选择它在时间和频率分辨率之间取了一个平衡折中适合大多数振动信号和生物信号。如果你需要提取瞬时相位——比如分析脑电的相位同步那优先用复值小波amor它的实部和虚部正交相位信息更干净。如果信号里的成分是很窄的频带脉冲比如齿轮箱的啮合频率bump小波的频率局部性最好能把两个靠得很近的频谱峰分开。3. 从时频图到频谱分析的完整 matlab 实现流程3.1 构造测试信号验证 CWT 是否正确的最短路径做频谱分析之前强烈建议先跑一个已知答案的合成信号。因为信号是你自己造的理想结果长什么样你心里完全有数如果 CWT 跑出来的图对不上那一定是你参数设错了趁早排查。常见的做法是构造一个频率随时间线性变化的 chirp 信号外加一个短时冲击clear; clc; close all; % 基础参数 fs 2000; % 采样率 2000 Hz T 2; % 信号时长 2 秒 N T * fs; % 总采样点数 t (0:N-1) / fs; % 时间轴 % 信号1chirp 信号频率从 50 Hz 线性扫到 300 Hz f0 50; f1 300; x1 chirp(t, f0, T, f1, linear); % 信号2在 1.2 秒处的短时冲击 x2 2 * sin(2*pi*200*t) .* exp(-((t-1.2)/0.005).^2); % 叠加成总信号并加入一点高斯白噪声模拟真实采集 rng(42); % 固定随机种子保证可复现 x x1 x2 0.05 * randn(size(t));逻辑说明chirp 信号的频率是连续变化的适合检验 CWT 的时间-频率对应关系。冲击信号用来验证 CWT 对瞬态事件的捕捉能力——FFT 里这种短时特征几乎不可能看得到。加入噪声是模拟真实传感器采集的环境顺便测一下 CWT 的抗噪能力。注意rng(42)这一行固定随机数生成器种子这样每次运行噪声序列完全一致方便对比参数改动前后的效果。不固定种子的话同样的代码每次跑出来的图都会有细微差异不利于调试。3.2 CWT 参数设置采样率、倍频程数和频率边界的配合拿到信号后CWT 参数怎么设这三个参数是关键。参数名推荐值说明fs实际采集的采样频率必须和信号真实情况一致否则频率轴全部漂移voicesperoctave16默认 10 够用16 更平滑32 计算量明显增大FrequencyLimits从关注的最低频率到 fs/2设得越窄计算越快图上频率分辨率越好% CWT 主程序 [cfs, frq] cwt(x, fs, ... voicesperoctave, 16, ... wavelet, morse, ... FrequencyLimits, [20, fs/2]); % 提取幅度谱复数 cfs 取模转成分贝值更直观 amp abs(cfs); amp_db 20 * log10(amp / max(amp(:)) eps); % 绘制时频图 figure(Color, w, Position, [100 100 1200 400]); imagesc(t, frq, amp_db); axis xy; colormap(jet); % 工程上常用 jet 色带冷色低、暖色高 clim([-40, 0]); % 只显示最高以下 40 dB 的成分 xlabel(时间 (s), FontSize, 12); ylabel(频率 (Hz), FontSize, 12); title(连续小波变换时频图幅度归一化, FontSize, 14); h colorbar; ylabel(h, 幅度 (dB), FontSize, 12);这里clim([-40, 0])是一个值得强调的细节。如果直接用线性幅度abs(cfs)画图能量强的高频冲击成分会完全压制低频的 chirp 信号图上只能看到一条亮线细节全被淹没。转成 dB 单位并把范围限制在最高值以下 40 dB相当于给画面做了一个动态范围压缩这是工程上处理频谱图的标准做法。3.3 从 CWT 结果中提取单一时间点的频谱切片时频图整体看趋势但如果要定量分析某个时刻的频谱构成怎么办比如想知道 0.8 秒处有哪些频率成分这时可以从cfs矩阵里抽一列出来画频谱图% 取出 0.8 秒最近的时间索引 [~, idx_t] min(abs(t - 0.8)); % 该时刻的频谱切片取这一时间列的所有频率点 slice_amp abs(cfs(:, idx_t)); slice_amp_db 20 * log10(slice_amp / max(slice_amp) eps); figure(Color, w); semilogy(frq, slice_amp); % 对数幅度更清晰 grid on; xlabel(频率 (Hz)); ylabel(幅度); title(sprintf(t %.2f s 时刻的频谱切片, t(idx_t))); xlim([20 fs/2]);逻辑说明cfs矩阵的列索引和时间轴一一对应通过min(abs(t - 0.8))找到最近的列号这一列就是该时刻的频域切片。它相当于短时傅里叶变换里某一帧的频谱但窗口类型和长度是由高频截止、低频截止动态决定的比 STFT 固定窗更适合宽频分析。提示如果信号里有 50 Hz 工频干扰这里会看到 50 Hz 处有一个稳定的窄峰如果观察到频率在漂移的谱峰那就是机器转动频率随负载变化。3.4 用 CWT 做频谱分析时的频率轴校准matlab 的cwt()输出的frq是归一化频率还是绝对频率取决于你有没有传fs。一个很容易踩的坑采样率是 1000 Hz信号里实际包含 250 Hz 成分但画的图上峰值出现在 0.25归一化频率数值上不明显容易误导判断。传了fs后frq直接就是 Hz 单位250 Hz 就标在 250 Hz 上。多花半秒钟传参省掉一小时的换算排查时间。还要注意不管是cwt还是cwtft返回的frq都是频率中心点的向量不是频率范围区间画图时直接作为 y 轴即可不需要位移半格。4. 边界效应、低分辨率区和参数微调CWT 最常踩的坑4.1 边界效应信号两端的数据要谨慎解释CWT 在计算每个尺度的小波系数时小波在时间轴上滑动靠近端点的时候小波有一部分伸到了信号外面。matlab 内部会根据Boundary参数决定如何处理边缘默认是对信号做对称扩展但这只是一种数值补丁端点附近的系数仍然不可靠。% 对比不同边界处理方式 [cfs_a, frq_a] cwt(x, fs, Boundary, symmetric); [cfs_b, frq_b] cwt(x, fs, Boundary, periodic); [cfs_c, frq_c] cwt(x, fs, Boundary, reflective); % 在时频图中把边缘区域标出来 figure(Color, w); subplot(1,3,1); imagesc(t, frq_a, abs(cfs_a)); axis xy; title(symmetric); subplot(1,3,2); imagesc(t, frq_b, abs(cfs_b)); axis xy; title(periodic); subplot(1,3,3); imagesc(t, frq_c, abs(cfs_c)); axis xy; title(reflective);边界处理参数的影响在低频段尤其明显低频对应大尺度小波被拉长覆盖的时间范围大边缘区域占比自然更高。如果你做的是比较小波系数能量、量化瞬时频率这类定量分析建议干脆把前面和后面各floor(10*scale_max)个点当作无效区只分析中间段。尺度最大值可以从frq(1)反推大概时间跨度duration 1/frq(1)秒。4.2 低频段的高频分辨率幻觉CWT 里低频段的频率分辨率好、时间分辨率差这是它和 STFT 最大的区别也是新手最容易误读的地方。举个例子10 Hz 和 12 Hz 两个成分在 CWT 图上可能会看到两条靠得很近的水平亮线但在某个时间点上它们的幅度值可能完全混在一起分不开。正确做法是要区分两个接近的频点看时频图上的俯视结构要精确定位某个事件的时刻看高频段的尖锐响应。信号特征用 STFT 合适用 CWT 更合适平稳正弦多频叠加是是但没明显优势频率调制、扫频信号一般是chirp 轨迹清晰短时冲击、瞬态故障较差是时间定位准确超低频长时间成分较难选窗长是低频分辨率高4.3 盲目调voicesperoctave带来的性能问题voicesperoctave控制频域采样密度默认 10 时每个倍频程里有 10 个频率点。调大到 32cfs的行数会增加约 3 倍计算时间显著上升。对一段 10 秒、采样率 48 kHz 的信号默认参数跑下来可能只要两三秒调到 32 可能要半分钟以上。实际项目经验先用默认参数快速预览整体时频结构确认频率范围设置没问题后再决定是否增大voicesperoctave。对于最终出图的场景——比如频谱分析报告或者论文插图——16 已经是视觉上足够平滑的值。另一个技巧是用cwtfreqparams()获取当前参数下的频率轴分布帮助判断频率分辨率是否满足需要wp cwtfreqparams(sv, 16); % wp 返回一个结构体包含频率轴、尺度、小波参数等这个方法适合调试时快速确认frq的范围和分布密度是否合理不用反复跑完整cwt。4.4 信号长度不足时的处理方案CWT 对信号长度没有硬性下限但太短的话低频段几乎没有有效数据。假设采样率 1000 Hz你想分析低到 1 Hz 的成分对应周期 1 秒至少需要几个周期的信号长度才可信比如 5 秒以上。如果只有 0.5 秒的数据那 1 Hz 的分析结果基本全是边界伪影。遇到这种短数据通常的做法有两个一是提高最低分析频率到FrequencyLimits的下限比如 10 Hz 以上二是对信号做零填充但这会显著增加计算量很多情况下并不值得。稳妥方案是重新审视数据采集配置让记录时长至少覆盖最低目标频率的 510 个周期。5. 用 CWT 找回 FFT 看不到的瞬态特征一个轴承振动案例把前面所有内容落到一个具体场景里轴承故障振动信号。这种信号的典型特征是低频转频成分每分钟转速对应的频率持续存在同时伴随高频的周期性冲击故障特征频率冲击本身有衰减振荡。用 FFT 看频谱只能看到宽泛的隆起很难定位冲击的发生时刻和重复间隔——而这两个信息恰好是判断轴承故障阶段的关键。% 模拟轴承内圈故障振动信号 fs 12000; T 1; t (0:fs*T-1) / fs; % 转频 30 Hz故障特征频率 120 Hz fr_rot 30; fr_fault 120; x_bearing sin(2*pi*fr_rot*t); % 转频成分 % 每 1/fr_fault 秒出现一次冲击冲击频率 800 Hz衰减系数 150 for k 0:fr_fault-1 tk k / fr_fault; idx find(t tk t tk 0.02); x_bearing(idx) x_bearing(idx) ... 0.8 * sin(2*pi*800*(t(idx)-tk)) .* exp(-150*(t(idx)-tk)); end % 加噪声降低信噪比到比较真实的地步 x_bearing x_bearing 0.1 * randn(size(t)); % CWT 分析 [cfs_b, frq_b] cwt(x_bearing, fs, voicesperoctave, 32, ... FrequencyLimits, [20, 2000]); figure(Color, w, Position, [100 100 1100 700]); % 第一行时域波形 subplot(3,1,1); plot(t, x_bearing); xlabel(时间 (s)); ylabel(幅值); title(轴承振动时域波形); xlim([0 1]); % 第二行FFT 频谱 subplot(3,1,2); Xf fft(x_bearing); f_fft (0:length(Xf)/2-1) * fs / length(Xf); plot(f_fft, 2*abs(Xf(1:length(Xf)/2))/length(x_bearing)); xlabel(频率 (Hz)); ylabel(幅值); title(FFT 频谱); xlim([0 2000]); % 第三行CWT 时频图 subplot(3,1,3); amp_b abs(cfs_b); amp_b_db 20*log10(amp_b/max(amp_b(:)) eps); imagesc(t, frq_b, amp_b_db); axis xy; colormap(jet); clim([-35, 0]); xlabel(时间 (s)); ylabel(频率 (Hz)); title(CWT 时频图可见周期性冲击); colorbar;运行这三行对比FFT 频谱图上你会在 30 Hz 处看到转频峰、在 800 Hz 附近看到一个鼓包但冲击的周期性、每秒多少次、有无间歇完全看不出来。CWT 时频图上则是另一番景象水平方向有一条贯穿始终的 30 Hz 亮线同时每隔约 0.0083 秒1/120 Hz会出现一条竖直的高频亮带衰减过程清晰可见。在这个案例里CWT 真正帮你做的事情有两个。一是直接数图中的竖条亮带数量快速验证故障特征频率二是看亮带的幅度变化如果故障冲击的幅度时大时小并呈规律性波动往往对应轴承某一损伤点在负载区内外交替。这一层提取能力是 FFT 做频谱分析完全给不出来的也正是matlab实现连续小波变换对信号做频谱分析、而不是简单用 FFT 处理的核心理由。进一步量化冲击间隔可以把 CWT 高频带的幅值包络提取出来再做一次 FFT得到包络谱的故障特征峰这属于包络分析的标准流程。一句话总结这个方法的价值CWT 不是要替代 FFT而是负责把 FFT 抹平的瞬态时间信息找回来。先用 CWT 粗看全局确认感兴趣的时频区域再用 FFT 做该区域内的精细频谱计算两种工具配合使用才算把频谱分析这件事做完整。本文还有配套的精品资源点击获取