简介面向MATLAB用户与信号处理学习者的连续小波变换与频谱分析资源包聚焦非平稳信号的时频联合分析适合科研、故障诊断与教学演示场景。压缩包仅4KB共3个m文件分别为一维连续小波变换、二维连续小波变换及分层小波计算脚本代码结构简洁便于直接调用或在此基础上修改尺度参数、小波类型并输出小波系数图。已有1023人学习下载适合需要快速上手CWT、理解小波系数含义及进行频谱特征提取的读者。资源围绕Morlet等小波函数的使用、cwt函数调用、时频图像绘制与功率谱分析展开可辅助完成信号去噪、模式识别等实验是入门连续小波变换并延伸至工程应用的实用工具。1. 连续小波变换不是要替代 FFT而是补上时频这层信息做频谱分析时FFT 是最先想到的工具但它的前提是信号在分析窗内平稳。真实信号比如电机启动、齿轮磨损早期、电力系统电压暂降频率成分随时间在变FFT 只告诉你有哪些频率不告诉你这些频率什么时候出现。短时傅里叶变换能给出时频图但窗口一固定低频和高频分辨率就被锁死。连续小波变换CWT通过压缩和平移母小波不再用固定窗口而是把时间-频率分解拆成时间-尺度分解再换算成频率低频用宽窗获得高频率分辨率高频用窄窗获得高时间分辨率。这是它在非平稳信号频谱分析里不可替代的原因。下面用 MATLAB 的 cwt 函数把 CWT 从概念推到可复现的频谱分析流程覆盖参数设置、瞬时频率提取和常见坑。2. 连续小波变换的数学骨架与 MATLAB 最小可跑代码2.1 从尺度和平移理解连续小波变换别死在积分公式上连续小波变换的定义式是W(a,b) (1/sqrt(a)) ∫ x(t) ψ*((t-b)/a) dta 是尺度b 是平移。初学者容易卡在这个积分上做频谱分析时只需要把它想成母小波 ψ(t) 是一段能量集中在局部时间的小波形尺度 a 决定它被拉伸还是压缩平移 b 决定它滑到信号的哪个位置。把 ψ 在每个尺度、每个位置跟信号做内积得到一张时频图横轴是时间纵轴是频率或尺度颜色是内积的模值。那个1/sqrt(a)是能量归一化系数用来保证不同尺度的小波在扫描时携带相同能量。没有它大尺度波形的系数会被整体抬高频谱分析会误判成低频段能量异常增强。理解到这层直接跳到 MATLAB 里跑通比反复推导积分快得多也够支撑后面的参数调试。2.2 用 matlab 的 cwt 函数3 行代码画出信号的时频图新版本调用 cwt 不需要手写滤波器组参数传采样频率就行。下面生成一个 30 Hz 扫到 100 Hz 的线性调频信号叠加一个 60 Hz 正弦时长 2 秒采样率 1000 Hzfs 1000; % 采样率 1000 Hz t 0:1/fs:2; % 0 到 2 秒共 2001 个采样点 sig sin(2*pi*(30*t 20*t.^2)) 0.4*sin(2*pi*60*t); cwt(sig, fs); % 不接收返回值时直接弹出时频图这段代码里sin(2*pi*(30*t 20*t.^2))的相位是30t 20t²对时间求导得到瞬时频率30 40tHz正好从 30 Hz 线性爬到 110 Hz。cwt(sig, fs)不接收返回值时MATLAB 直接弹图纵轴已经是 Hz颜色是小波系数模值沿时频平面的分布。图上能看到一条斜向上的亮带和一条水平 60 Hz 的亮带这就是读谱的基本形态。要拿数据继续分析改成接收返回值[wt, f] cwt(sig, fs); % wtlength(f) 行length(sig) 列wt每一行对应f里的一个频率点每一列对应一个采样时刻f的单位是 Hz。提示老版本里 cwt 需要手写尺度向量新语法在cwt(sig, fs)一步完成写新代码时不要混用老语法。2.3 小波基选型Morse、复 Morlet、Mexican Hat 怎么配信号cwt 默认用解析 Morse 小波它对大多数工程信号不需要额外调参。观测目标是连续变化的瞬时频率时换成amor也就是解析 Morlet 小波时间定位更锐利代价是频率方向稍微模糊。如果目标是抓突变比如局部放电脉冲、敲击回波用mexh这类实小波更直接但它是实数小波输出没有相位信息频谱分析时只适合看幅值。参数值小波类型适合场景注意点morse默认解析 Morse 小波通用非平稳信号频谱分析调参少边界效应中等amor复 Morlet 小波瞬时频率估计、振动信号时间分辨率高频率分辨率略降mexh实 Mexican Hat 小波冲击检测、突变定位无相位信息只读幅值选型时不用纠结太久先跑默认morse如果时频图上的频率线被抹得看不清或者你想精确读某一条脊线的瞬时频率再换amor对比。振动频谱分析场景里amor出现的频率比我预期高很多因为工程关心的往往是转频及其谐波的时间演化而不是某个孤立频点的绝对功率。2.4 频谱分析时读的是系数模值不是相位小波系数是复数幅值代表该时刻该频率成分与母小波的相似程度相位则和信号波形相位耦合在一起。做频谱分析时绝大多数情况只看模值或模值平方相位留给专门的相位分析。读谱时我喜欢先把功率分布算出来再取每一时刻的峰值频率[wt, f] cwt(sig, fs); powerMap abs(wt).^2; % 时频功率分布 [~, maxIdx] max(powerMap, [], 1); % 沿频率方向取最大 instFreq f(maxIdx); % 每个时刻的主频估计 plot(t, instFreq);powerMap每一列是一个时刻的功率随频率分布max沿第一个维度取得到每个采样时刻能量最强对应的频点索引转成频率就是瞬时主频的粗略估计。这个方法和第四章的tfridge思路一致先用手写版本理解逻辑再用工具箱函数做严谨版。3. 连续小波变换的 4 个必调参数频率分辨率、范围与边界处理3.1 VoicesPerOctave频率网格加密还是稀疏cwt 的频率轴不是线性等间隔而是按倍频程对数分布。MATLAB 用VoicesPerOctave控制每个倍频程里放几个频率点默认 10。需要精细读频时我会调到 24 或 32代价是系数矩阵行数变多内存和时间线性上涨。用下面的循环直观感受频率网格密度for v [4 10 24 48] [~, fv] cwt(sig, fs, VoicesPerOctave, v); fprintf(Voices%2d 频率点数%4d 中位频点间隔%.3f Hz\n, ... v, length(fv), median(diff(fv))); end对数网格下高频段间隔大低频段间隔小median(diff(fv))反映整体密度。对 1000 Hz 采样率、2 秒信号VoicesPerOctave从 10 提到 48频率点数从几十涨到两百多运行时间明显变长。我做深度学习特征时经常回到 10 甚至更低因为高密度时频图像对分类器并不总是更友好反而容易过拟合。3.2 FrequencyLimits 与采样率频率范围不是你随便填的cwt 的频率上限由采样率和小波带宽共同决定最高有效频率到不了fs/2因为小波在接近奈奎斯特频率时已经不能完整成形实际有效范围通常在采样率的四成以下。较新的 MATLAB 版本支持直接指定频率范围[wt, f] cwt(sig, fs, FrequencyLimits, [1 120]);这样滤波器组只生成 1 到 120 Hz 之间的频率点系数矩阵行数减少内存占用下降。如果你的版本不支持这个参数用全频带计算画图后加一行ylim([1 120])限定显示范围两种做法的读谱效果基本一致只是后者慢一点。3.3 边界效应信号两端那些竖条纹怎么来的小波扫描到信号首尾时有一部分波形伸出数据范围。MATLAB 默认做反射延拓避免直接截断造成的突变但延拓出来的内容和真实信号不一致结果就是时频图两端出现竖直的亮条或暗条。越低频小波越长受影响的范围越大。处理办法是直接裁掉两端margin round(0.1 * numel(t)); % 舍弃两端约 10% 的时间点 coreIdx (margin 1):(numel(t) - margin); sigCore sig(coreIdx); tCore t(coreIdx); [wtCore, fCore] cwt(sigCore, fs, VoicesPerOctave, 24);裁边之后时频图两端干净很多代价是时间长度缩短了 20%。对长记录信号无所谓对短脉冲信号就要权衡边界伪影和有效时长只能选一头。真实数据里我一般先不裁边做一次 cwt 看整体确认关注的事件不在两端再决定要不要裁。3.4 参数对比表什么时候改哪个参数参数默认改大改小副作用VoicesPerOctave10频率读数更精细计算更快数据量线性增长FrequencyLimits自动关注窄带覆盖全频带高频细节丢失裁边比例无边界干净保留时长时间信息减少采样率 fs无提高上限降低存储时间分辨率变化这几个参数不是独立的。提高VoicesPerOctave让频率网格变密但数据量涨裁掉边界的本质是牺牲时间长度换掉伪峰去噪也类似只取关注频带能压低无关分量。做频谱分析前先把信号长度、采样率、关注频带写在注释里再决定参数比逐个盲试要快得多。4. 工程频谱分析实战用连续小波变换提取振动信号瞬时频率4.1 场景启动过程振动数据导入 CSV 后怎么做电动机或泵启动阶段转速从零爬升振动主频跟着变化。这个阶段做 FFT频谱是一片带宽模糊的包络读不出转频随时间怎么变。常见做法是先读 CSV再看时间单位最后跑 cwtdata readmatrix(pump_startup.csv); % 第一列时间第二列振动加速度 t data(:, 1); acc data(:, 2) - mean(data(:, 2)); % 去直流分量 fs round(1 / median(diff(t))); % 用采样间隔中位数反推采样率readmatrix会自动跳过 CSVT 表头里的文本行第一列是时间、第二列是振动加速度时这种读法最省事。median(diff(t))比直接用t(2)-t(1)稳因为记录仪偶尔丢点会导致个别间隔异常中位数能避开这些尖刺。提一句常见错误如果时间列单位是毫秒fs会被算成 1000 倍的实际采样率时频图的频率轴整体偏高先画一段原始信号确认fs量级再做 cwt。4.2 用 tfridge 从 cwt 系数中提取瞬时频率转频的时变曲线可以通过脊线提取得到工具箱里的tfridge干的就是这个事[wt, f] cwt(acc, fs, VoicesPerOctave, 32); [fridge, ~] tfridge(wt, f); % 每个时刻取能量最大的频率 plot(t, fridge); xlabel(时间 (s)); ylabel(瞬时频率 (Hz));tfridge默认在每一列里找模值最大的频率点输出fridge的时间长度和wt列数一致可以直接和t对齐画图。如果脊线跳变太多加惩罚项平滑[fridgeS, ~] tfridge(wt, f, 1, 0.05); % 第4个参数越大曲线越平滑第三个参数是脊线数量第四个是频率方向上的惩罚系数。0.05 是常用起点太大把真实的频率变化也压平了太小噪声毛刺全留下。启动过程转频爬升是缓变信号0.05 到 0.1 之间通常合适。4.3 频带能量随时间变化故障诊断的时频特征瞬时频率能看出转频轨迹但很多故障特征不在主频上而在某个固定频带的能量抬升。比如齿轮啮合频率附近的边带能量会随着磨损程度缓慢增强。用频带能量时只需对系数矩阵做索引求和band f 90 f 110; % 关注频带 eband sum(abs(wt(band, :)).^2, 1); % 每个时刻该频带总能量 plot(t, eband);band是逻辑索引选出频率向量里落在 90 到 110 Hz 的行sum(..., 1)沿频率方向求和得到一个时间序列。这个序列描述的是「该频带能量随时间怎么变」连续几个周期都比基线高就值得警惕。实际工程里我会同时提取多条频带的能量曲线拼成特征矩阵送给分类器或时序模型。特征提取方式适用场景瞬时频率曲线tfridge(wt, f)转频监测、转速估计频带能量时序sum(abs(wt(band,:)).^2, 1)窄带故障特征时频图像imagesc后重采样深度学习 CNN 输入时频图作为图像输入是深度学习做故障诊断的常见路线先把 cwt 结果画成图再缩放到统一尺寸喂 CNN。如果只用标量特征做时序分类上面的频带能量和瞬时频率就够了bilstm 这类模型直接吃这些序列。4.4 在与 FFT 的对照中你就知道 cwt 补了哪块同一段启动数据直接做 FFT频谱是整个时段的平均转频从 50 爬到 150 Hz结果是一个宽包络峰值大约对应平均转速既看不出爬升方向也看不出某个频率什么时候被激发。代码是这样的N numel(acc); win hann(N); % 汉宁窗抑制泄漏 spec abs(fft(acc .* win)); fAxis (0:N-1) * fs / N; plot(fAxis(1:floor(N/2)), spec(1:floor(N/2))); xlabel(频率 (Hz));FFT 适合回答「整体有哪些频率成分」cwt 适合回答「这些成分什么时候出现、怎么变化」。振动频谱图怎么分析常规流程就是先 FFT 找异常频段再用 cwt 定位异常出现的时刻和演化路径。两步配合比单用任何一个都完整。5. 连续小波变换用不稳的 3 个坑和一个先跑谱线再谈模型的验证技巧5.1 版本与工具箱cwt 报错先查这两个新代码请用新语法cwt(x, fs)。如果报Undefined function cwt多半是工具箱缺失命令行里敲ver(signal) ver(wavelet)返回空说明对应工具箱没安装cwt属于 Wavelet Toolboxfft属于 Signal Processing Toolbox。另一个常见问题是老代码直接用cwt(x, scales, wname)新版本里这种语法已经不被推荐报错提示会让你改写成滤波器组方式。看到参数数量不对先检查是不是新旧语法混用。5.2 低信噪比时频图上看不到线怎么办低信噪比下系数被噪声托底默认颜色映射会把弱信号压成一片暗色。把系数转成 dB 显示对比度会明显改善imagesc(t, f, 20*log10(abs(wt) / max(abs(wt(:))))); axis xy; % y 轴从小到大 colorbar;分母用全局最大值显示的是相对强度弱信号从 -30 dB 变成可辨认的亮斑。如果还是看不清换amor小波试一次复 Morlet 的时间定位更锐利对短时弱信号比默认 Morse 更敏感。注意 dB 显示只改观感不改数据本身后续提取特征仍然用原始wt。5.3 用合成信号验证脊线读数对不对先有基准再说调试参数最怕没有基准。先用已知瞬时频率的调频信号验证整条链路fs 1000; t 0:1/fs:1; trueF 50 100*t; % 真实瞬时频率 50 到 150 Hz phase 2*pi*cumsum(trueF)/fs; % 相位是频率的积分 sig sin(phase); [wt, f] cwt(sig, fs, VoicesPerOctave, 32); [fridge, ~] tfridge(wt, f); rmse sqrt(mean((fridge - trueF).^2)); fprintf(脊线提取 RMS 误差%.3f Hz\n, rmse);cumsum(trueF)/fs是对频率序列做数值积分生成调频信号的正确做法比直接写sin(2*pi*(50*t 50*t.^2))更不易算错。误差如果在几个 Hz 以内说明参数和代码链路没问题。把这段脚本存成verify_cwt.m每次调完VoicesPerOctave或FrequencyLimits先跑一遍确认误差没变大再放到真实振动数据上。本文还有配套的精品资源点击获取