简介：MATLAB短时傅里叶变换与Morlet小波变换实现包，面向信号处理方向的学生和科研人员，帮助解决非平稳信号时频分析中的常见问题。包内共2个m文件，压缩包大小仅816B，是轻量级的学习演示脚本，虽然体积小，但代码注释清晰，可直接在MATLAB中运行观察输出。目前已有815人学习下载，适合入门阶段对比理解STFT与Morlet小波变换的适用差异。脚本覆盖短时傅里叶变换的窗函数选取、分段加窗与时频图绘制流程，以及Morlet小波变换的复小波基、尺度参数调节和连续小波变换函数cwt用法，通过实例展示两种方法在时间分辨率与频率分辨率上的取舍。读者按注释修改参数即可迁移至语音、地震或心电等信号分析场景，有助于从原理过渡到实践。
1. 为什么非平稳信号分析要同时准备短时傅里叶变换和 Morlet 小波
普通傅里叶变换只能告诉我们一段信号里有哪些频率，却说不清这些频率出现在什么时候。对于语音、机械振动、脑电这类非平稳信号，时间信息恰恰是重点。短时傅里叶变换（STFT）把信号切成小段，逐段做 FFT，得到时间-频率二维分布；Morlet 小波则用一个带高斯包络的复指数小波去匹配信号局部特征，低频长窗、高频短窗，分辨率更灵活。这套资料里的 WORK3.m 和 work4.m 正好演示了两种方法在 MATLAB 里的落地方式，适合正在做信号分析、故障诊断或生物电处理的人。这里先把两者放在一个工程环境里拆开讲，然后再回到脚本看参数和排错。
2. 短时傅里叶变换在 MATLAB 里的窗口、步长与 spectrogram 参数
2.1 从 FFT 到 STFT：窗口长度的物理含义
假设一列信号 x[n] 的采样率 fs 为 1000 Hz，整体持续 10 秒，直接做 FFT 只能得到 0～500 Hz 的频谱，峰值可能出现在 100 Hz，但丢失了这个 100 Hz 成分是出现在第 2 秒还是第 8 秒。STFT 的处理方式是先选一个长度为 N 的窗函数，比如汉明窗，在时间轴上滑动，每滑动一步就把窗内的信号截取出来做一次 FFT，得到一个列频谱。所有列频谱按时间顺序排起来，就形成二维时频矩阵，横轴时间、纵轴频率、颜色代表幅值。
这里最关键的参数是窗长 N。令 fs 为采样率，STFT 的频率分辨率由 Δf = fs / N 决定。N 越大，频率分辨得越细；但时间窗越宽，一个时刻附近的信号被平均得越厉害，时间分辨率就越差。这与小波分析里“低频看频、高频看时”的思路形成鲜明对比。实际工程里我一般先确定最低关心频率。如果想区分 10 Hz 和 12 Hz 两个相邻谱峰，至少需要 Δf ≤ 2 Hz，那么在 fs=1000 Hz 时，N 不要小于 500，取 512 点比较稳妥；如果只关心 50 Hz 以上的工频及谐波，256 点就够了。
窗函数的形状决定了频谱泄漏和旁瓣。矩形窗主瓣最窄，但旁瓣衰减慢，容易把强频率成分泄漏到相邻频点；汉宁窗、汉明窗旁瓣衰减快，适合大多数工程场景；布莱克曼窗衰减更快，但主瓣更宽，谱峰显得钝。下表是常见窗的选择依据。
| 窗类型 | 主瓣宽度（相对矩形窗） | 旁瓣峰值 | 典型用途 |
|---|---|---|---|
| 矩形窗 | 1.0 | -13 dB | 瞬态明显、窗内近似平稳 |
| 汉宁窗 | 2.0 | -31 dB | 通用信号分析，默认首选 |
| 汉明窗 | 2.0 | -41 dB | 语音和窄带信号 |
| 布莱克曼窗 | 3.0 | -58 dB | 需要极低旁瓣的功率谱 |
第一次做时频分析时，不要迷信复杂窗。先用汉宁窗，把时间分辨率调出来，观察谱线是否泄漏，再换窗。换窗带来的是旁瓣形状变化，而窗长带来的则是时间和频率分辨率取舍，后者影响大得多。
2.2 用 spectrogram 画出一张可用的时频图
MATLAB 里最常用的是spectrogram，它把分帧、加窗、FFT、重叠全部封装好。假设工作区已有信号t1，典型调用如下：
% WORK3.m 的常见结构：t1 是含噪非平稳信号 fs = 1000; % 采样率，单位 Hz t = (0:length(t1)-1)/fs; % 时间轴 win = hamming(256); % 256 点窗，约 0.256 秒 noverlap = 220; % 相邻窗口重叠 220 点，步长 36 点 nfft = 512; % 频域补零到 512 点 [S, f, tp] = spectrogram(t1, win, noverlap, nfft, fs); imagesc(tp, f, 20*log10(abs(S)+eps)); % 用 dB 显示动态范围 axis xy; xlabel('Time (s)'); ylabel('Frequency (Hz)'); colorbar;代码里win是窗函数向量，长度就是窗长 N；noverlap是重叠点数，窗口滑动步长等于length(win) - noverlap，本例是 36 点。重叠越多，时间轴越平滑，但相邻帧相关性也越高，计算量上升。nfft是 FFT 点数，一般取大于窗长的 2 的幂，这里 512 比窗长 256 大，相当于在频谱上插值加密，并不会改变真实频率分辨率。显示时用 dB 转换，是因为原始abs(S)的动态范围可能跨越几个数量级，强峰会把弱成分压得看不见。
如果使用较新版本的 Signal Processing Toolbox，可以换成stft函数：
[S2, f2, tp2] = stft(t1, fs, 'Window', hamming(256), ... 'OverlapLength', 220, 'FFTLength', 512);stft和spectrogram的命名差异容易让人困惑：spectrogram返回的时间向量对应每段窗口的中心时刻，而stft也默认返回中心时刻，两者画出的图基本一致。主要区别在于输出顺序和默认参数。用help spectrogram可以看到，老版本里如果省略返回值，函数会直接画三维瀑布图；新版本推荐使用带返回值的写法，再自己用imagesc或surf显示。
2.3 判断窗口参数是否合理的三条标准
拿到 STFT 图后，先不要急着调颜色。通常有三条判断标准。第一，时频图里不应出现所有频率同时被加粗的横向亮带，如果有，说明窗口太短，频率泄漏严重。第二，瞬态冲击在时间轴上应该表现为竖直的窄亮线，如果亮线横向拖尾，说明时间分辨率不足。第三，对平稳正弦信号，谱峰宽度应接近 Δf=fs/N，比如 fs=1000、N=256 时，理论分辨率约 3.9 Hz，谱峰宽度若明显大于这个值，就要检查窗类型和补零设置。用这三条标准检查 WORK3.m 的输出，比反复改色标更有效。
3. Morlet 小波变换的 cwt 实现与尺度-频率换算
3.1 Morlet 小波的时频特性：复指数、高斯包络和中心频率
Morlet 小波由带高斯包络的复指数函数构成，连续形式可以写成 ψ(t) = (1/sqrt(pi)) * exp(i·ω0·t) * exp(-t²)。高斯包络让小波在时间上有很好的局域性，复指数又让它带有明确的频率偏向。和 STFT 固定窗长不同，小波变换通过伸缩母小波改变分析带宽：尺度越大，小波在时间上越被拉宽，对应中心频率越低，频率分辨率越细；尺度越小，小波越窄，时间分辨率越好。所以整个时频平面在不同频段上的分辨率并不均匀，这对低频持续长、高频持续短的自然信号非常友好。
在较老的 MATLAB 资料里，常能看到'morl'这种写法，它代表默认中心频率的 Morlet 小波。centfrq('morl')大约返回 0.8125，这是尺度转频率的关键参数。需要特别注意的是，这个中心频率是小波本身的中心频率，不是信号中的物理频率。物理频率还与采样率有关：同样的尺度，在 fs=1000 和 fs=2000 两种采样率下，对应的实际频率差一倍。很多刚接触小波的人会把尺度轴直接当成频率轴，这是最容易犯的错误。
3.2 cwtfilterbank 与老式 cwt：两种调用方式的差异
新版 MATLAB 建议用cwtfilterbank构造滤波器组，再调用wt计算系数。这样做的原因是它把尺度序列、归一化、边界处理都封装好，直接返回物理频率轴，省去手工换算。典型结构如下：
% work4.m 的常见结构：用解析 Morlet 小波处理 t1 fs = 1000; fb = cwtfilterbank('SignalLength', length(t1), ... 'SamplingFrequency', fs, ... 'Wavelet', 'amor'); % 'amor' 是解析 Morlet 小波 [cfs, f] = wt(fb, t1); % cfs 复数矩阵，f 为物理频率 imagesc(t, f, abs(cfs)); set(gca, 'YScale', 'log'); % 小波分析习惯用对数频率轴 xlabel('Time (s)'); ylabel('Frequency (Hz)'); colorbar;这里Wavelet指定为'amor'，也就是解析 Morlet 小波。返回的f已经是一组物理频率向量，行数由滤波器组内部的尺度数以及VoicesPerOctave决定。这个参数默认是 10，表示每个倍频程内放置 10 个尺度。把VoicesPerOctave提高到 16，时频图会明显更平滑，但矩阵行数、计算时间和内存占用都同步上升；对长信号，我会先保持默认 10，确认频率范围正确后再加密。
老代码里常见格式是cwt(t1, scales, 'morl')，新版已经用新语法替代。但旧脚本仍然需要看懂，因为 WORK3.m、work4.m 这类资料很可能来自几年前。手动换算伪频率的方式是：
fc = centfrq('morl'); % 中心频率约 0.8125 dt = 1 / fs; % 采样周期 scales = 1:128; % 尺度序列 freqs = fc ./ (scales * dt); % 伪频率向量这行代码的逻辑是：尺度为 a 时，小波被拉伸 a 倍，其中心频率变为 fc/(a·dt)。换算结果称作“伪频率”，因为 Morlet 小波有带宽，不是单频，实际匹配到的频率会在该值附近波动。用scal2frq(scales, 'morl', dt)可以得到同样结果。如果直接用cwtfilterbank，这些换算已经被封装在内部，但不代表可以忽视它的存在，因为修改采样率后，f向量会自动变化，而老式脚本里的scales不会。
3.3 尺度范围如何选：先定最低频率，再定最大尺度
选尺度序列时，先要明确关心频率范围。最大尺度 amax 由最低关心频率 fmin 决定：amax = fc / (fmin · dt)。例如 fs=1000，fmin=1 Hz，dt=0.001，fc=0.8125，那么 amax 约为 812。直接用 1:64 这种固定范围，很可能覆盖不到低频段。常见做法是把尺度按等比数列分布，而不是等间隔，因为小波在低频段每个尺度对应的频率间隔很小，等间隔会造成低频分辨率过剩、高频分辨率不足。
cwtfilterbank自动生成的尺度序列更适合工程使用，老式cwt则需要自己构造等比数列，一般用2.^(1:0.2:8)这类形式。除了尺度，还要注意边界效应。信号两端的小波系数，特别是低频成分，受边界截断影响严重，表现为时频图两端出现异常亮纹。判断边界区的一个简单方法是把一小段零值拼在信号前后，重新做变换，对比前后两次结果中系数开始明显变化的区域，这个区域就是不可信边界。我一般只保留时频图中间 90% 区域作为有效分析区。
4. 同一段 t1 信号上 STFT 与 Morlet 小波的时频图像差异
4.1 固定窗 vs 恒 Q 窗：分辨率在时频平面上的分布
STFT 对整段信号使用同一个窗长，所以在时频图上，每个频率处的分辨率相同，表现为均匀栅格。Morlet 小波则改变分析窗宽：频率越低窗越长，频率分辨率更高；频率越高窗越短，时间分辨率更高。这种特性称为恒 Q 分析，Q 值近似等于中心频率与带宽的比值。工程上不需要死记定义，只需要看两个现象：低频段小波频谱比 STFT 更细，高频段小波时间定位比 STFT 更敏锐。
除了幅值，小波系数还是复数，包含相位信息。STFT 得到的S也是复数，但大多数人只看abs(S)，忽略相位。小波系数的相位可以用于瞬时频率估计或相位同步分析，在脑电、语音和机械故障诊断里都有用处。实际处理时，abs(cfs)反映能量，angle(cfs)反映相位，如果只画幅值图，等于丢弃了一半信息。
4.2 从 WORK3.m 和 work4.m 的输出反推脚本结构
两个脚本的命名基本暗示了作业或实验顺序：WORK3.m 是第一版 STFT 实现，work4.m 是小波版本。按常见写法，两者前几行会共用数据加载部分，比如load t1，然后各自调用工具箱函数。区别集中在中间几行：
% WORK3.m 里类似： [S, f, tp] = spectrogram(t1, hamming(256), 224, 512, fs); imagesc(tp, f, 20*log10(abs(S))); % work4.m 里类似： [cfs, f] = cwt(t1, 'amor', fs); imagesc(t, f, abs(cfs));从这段对比可以看出，STFT 的代码更直白，所有参数暴露在函数调用里，适合精细化调节；小波代码更简洁，但内部尺度序列由SignalLength和SamplingFrequency隐式决定。如果两个脚本针对同一段 t1，时间轴必须统一，否则对比时频图会错位。我一般会写一个小函数把两者输出对齐到相同时间网格：
tp_common = linspace(0, (length(t1)-1)/fs, 500); S_interp = interp1(tp, abs(S)', tp_common)'; cfs_interp = interp1(t, abs(cfs)', tp_common)';interp1把不同时间精度的谱线统一重采样到同一组时间点，这样后面做相关性比较时矩阵维度一致，不会因为列数不同报错。需要留意的是，重采样会引入插值误差，用于比较峰值位置是可行的，比较绝对幅值则不建议。另一个细节是归一化：新版cwtfilterbank默认对每个尺度做 L1 归一化，不同尺度的系数幅值可以直接比较；而旧版手工构造小波时容易忽略这一步，导致低频系数天然偏大，这不是信号本身的特征。
4.3 什么时候用 STFT，什么时候用 Morlet
下面这种分法可以应对大部分分析场景。
| 因素 | STFT | Morlet 小波 |
|---|---|---|
| 频率分辨率 | 各频段恒定 | 低频好、高频差 |
| 时间分辨率 | 各频段恒定 | 高频好、低频差 |
| 突变定位 | 依赖窗长，容易模糊 | 高尺度下更准 |
| 计算量 | 小，适合在线 | 大，适合离线 |
| 参数可调性 | 窗长、重叠、FFT 点数 | 尺度数和中心频率 |
| 典型应用 | 振动监测、语音端点检测 | 脑电、地震、故障诊断 |
工程中常把两者结合：先用 STFT 快速浏览整个频带的能量分布，圈出问题频段，再用 Morlet 小波对圈定频段做精细时频分析。这不是绕路，而是用小波的变分辨率特性去替代反复调 STFT 窗长的成本。参数层面，STFT 里的nfft和窗长决定谱线平滑度，小波里的尺度序列决定分析频段覆盖范围，两者没有一一对应关系，但都要满足奈奎斯特约束，最高分析频率不超过 fs/2。
5. 从 WORK3.m 到 work4.m：调参顺序、伪频率验证与常见报错
5.1 小波频率轴对不对，用单频正弦验证
拿到 work4.m 后，不要直接分析 t1。先构造一个已知频率的正弦信号验证频率轴是否正确：
fs = 1000; t = (0:1000)/fs; x = sin(2*pi*50*t); % 50 Hz 单频 [cfs, f] = cwt(x, 'amor', fs); [~, idx] = max(abs(cfs(:))); % 矩阵展平后最大位置 [row, ~] = ind2sub(size(cfs), idx); fprintf('探测到峰值频率: %.2f Hz\n', f(row));如果输出接近 50 Hz，说明小波频率轴正确；如果差很多，先查fs是否传给了滤波器组，再看有没有手动改过VoicesPerOctave。这里用max(abs(cfs(:)))只适合单频信号，实际信号需要按行扫描峰值。
5.2 两个脚本共用的排错清单
STFT 脚本里最常见的报错是noverlap大于等于窗长，或者length(win)大于信号长度，这时 MATLAB 会提示窗口超出数据范围。小波脚本里常见的报错是输入信号为整数类型，比如int16，部分版本会提示数据类型不支持，解决办法是double(t1)。另外，cwtfilterbank的SignalLength必须与实际信号长度一致，修改过加载数据后忘了同步修改，会出现维度不匹配的红色报错。检查顺序是：采样率、信号类型、信号长度、工具包路径。
5.3 一张图里同时放两种结果的排版技巧
用subplot(2,1,1)和subplot(2,1,2)并排显示时，色标范围要统一，否则视觉上无法比较。STFT 图的频率轴线性显示即可，小波图建议用set(gca, 'YScale', 'log')对数轴，这样低频段的细节才能展开。调参时先固定时间轴，再改窗长或尺度，观察同一时刻谱峰的变化。完成这些验证后，再回头处理 t1 的噪声和瞬态成分，你就能知道自己真正需要的是时间分辨率还是频率分辨率了。
本文还有配套的精品资源，点击获取
