简介：这份MATLAB时频分析程序包面向信号处理初学者与工程师，覆盖短时傅里叶变换、小波变换、Wigner-Ville分布及EMD/EEMD等常见方法，配套大量带exa编号的示例脚本，可系统学习时频分析原理与实现。压缩包共40个文件，以37个.m源程序为主，辅以2个.mat数据和1个txt说明，整体仅50KB，便于快速部署与运行。已有521人学习下载，内容包含多组可执行的仿真实验，如变化检测、多普勒信号分析、EEG/ECG数据处理等，代码注释与数据文件相互对应，既适合入门者按例程理解算法细节，也可供进阶用户直接调用或改造。通过学习这些程序，可提升MATLAB编程能力并加深对时频域分析技术的掌握。
1. 时频分析程序在 matlab 里到底解决什么问题
一条语音、一段振动波，直接做 FFT，得到的是整个时间段内频率含量的平均值。信号在某一瞬间发生的频率变化，会被摊平、糊成一片。时频分析程序要解决的就是把幅度、频率、时间三个维度同时呈现出来，让工程师能看到频率随时间演化的轨迹。在 matlab 里实现时频分析，最常见路线是先 spectrogram 建立基线，再上连续小波变换补低频细节；信号呈现强非平稳、非线性特征时，走 EMD 分解加希尔伯特谱。三条路线没有绝对的先进与落后，参数选型决定效果。后面的章节按这条路径推进，每段都给出能直接跑通的最小程序和容易踩的参数坑。
2. 用 spectrogram 搭出第一版时频分析程序：窗长与分辨率怎么权衡
STFT 的思路一句话：把信号切成一帧一帧，每帧加窗做 FFT，把谱向量按时间顺序排成二维矩阵。matlab 的spectrogram把这个流程压成了一个函数，省事，但参数封装得深，理解不到位容易把坐标轴物理意义弄错。我习惯先合成一个有明确频率结构的测试信号，再拿它去检验程序。
2.1 先造一个带调频和单频分量的测试信号
fs = 1024; % 采样率，单位 Hz ts = 0:1/fs:2-1/fs; % 2 秒时长，1024 点每秒 x = chirp(ts, 40, 2, 220, 'quadratic', [], 'convex'); % 频率 40 到 220 Hz 的二次调频 x = x + 0.15 * sin(2*pi*320*ts); % 叠加一个稳定的 320 Hz 正弦分量 x = x + 0.03 * randn(size(ts)); % 加入弱噪声，模拟工程信号这个测试信号里有两个值得时频分析注意的结构：二次调频分量的瞬时频率随时间连续变化，用来检验算法对非平稳频率的追踪能力；320 Hz 单频分量用来检查频率定位是否准确。小噪声存在则是为了观察底噪会不会在时频图上形成伪峰，这是实际采集信号最常见的干扰源。
2.2 spectrogram 最小可用程序
nwin = 256; % 窗长，单位样本数 noverlap = nwin - 32; % 相邻窗重叠的样本数 nfft = 1024; % FFT 点数，可大于窗长 [ps, f, tp] = spectrogram(x, hann(nwin), noverlap, nfft, fs, 'yaxis'); figure; imagesc(tp, f, 10*log10(abs(ps) + eps)); axis xy; colormap('jet'); colorbar; xlabel('时间 / s'); ylabel('频率 / Hz');ps是复数谱矩阵，每个元素对应一对时间窗中心与频率网格；tp是各帧窗中心所在时刻，f是频率轴。显示时用10*log10(abs(ps)+eps)转成 dB 值，不加eps会在全零频点上留下空白块，局部幅值差距大的程序里特别明显。
三个容易混淆的参数单独说明。noverlap在 matlab 里是样本数，不是百分比，设成nwin-32意味着每次窗向后滑动 32 个样本，相邻窗重叠 87.5%。nfft大于nwin时，FFT 不足部分自动补零，但这只是让频谱曲线更平滑，不会把两根本来要挨在一起的谱线分开。fs必须填真实采样率，漏填或填错，频率轴会整体偏移，这个错误在后文还会再出现一次。
2.3 窗长、重叠率与频率分辨率的工程设定
先说结论：STFT 的频率分辨率由窗长决定，时间分辨率也由窗长决定，两者是同一次切割行为的两个侧面。采样率 1024 Hz、窗长 256 样本时，一窗只有 0.25 秒，能分辨的两个频率分量至少要差约 4 Hz；如果想看清 40 Hz 与 41 Hz 这两个很接近的分量，窗长要拉到 1 秒量级。nfft解决不了这件事，这不是数值精度问题，而是窗的有限时宽决定了它能分辨的最小频率间隔。
窗函数的选择同样重要，它影响主瓣宽度和旁瓣高度。下面是工程里常用的四类窗：
| 窗函数 | 主瓣宽度 | 旁瓣抑制 | 适用场景 |
|---|---|---|---|
| hann | 中等 | 约 -31 dB | 默认首选，时间与幅值精度均衡 |
| hamming | 略窄 | 约 -43 dB | 旁瓣要求更高，频率泄漏要求更严时 |
| blackman | 较宽 | 约 -58 dB | 大动态范围，弱信号被强分量掩盖时 |
| kaiser | 可调 | 可调 | 需要精确控制主瓣与旁瓣折中，beta 参数据场景调 |
提示：调整分辨率时，优先改窗长而不是 nfft。短窗看全局演变，长窗看精细频率结构，noverlap 保持 75% 到 90% 可以缓解时间轴上的阶梯感。
3. 用连续小波变换做时频分析：cwt 的频率换算与 voicesperoctave
STFT 的窗长固定，是它的根本矛盾。低频分量周期长，需要长窗才能积累足够周期数，否则频率分辨率很差；高频瞬态又需要短窗才能定位到具体时刻。固定窗长无论怎么选，都只能在一个频段上表现良好。
连续小波变换用“尺度”代替“窗长”。尺度小时，小波在时间上被压窄，适合抓高频细节；尺度大时，小波在时间上展宽，天然适合作低频分析。matlab 的cwt函数封装了这套过程，输出直接是“时间-频率”二维结果，而不需要像老代码那样先算尺度再手工换算。
3.1 为什么 STFT 在低频段看不清
用第 2 章的参数计算一下：fs = 1024，nwin = 256，频率分辨率约 4 Hz。一个 40 Hz 分量的周期是 25 ms，一个 256 样本窗里有 10 个周期；而 320 Hz 分量在一个窗里有 80 个周期。周期数量决定了谱估计的稳定程度，低频分量在窗内“信息量”天然不足，固定窗长下无论怎么调 nfft 都无法改善。
CWT 的处理方式是不再给所有频率配同一个时间窗。低频用大尺度，时间窗自动加宽，周期数够了；高频用小尺度，时间窗自动压窄，瞬态定位更准。这个自适应特性不是一种额外优化，而是时频分析方法里处理非平稳信号的基本需求。
3.2 用 cwt 画最小可用时频图
[cfs, frq] = cwt(x, fs, 'amor'); % 解析 Morlet 小波，x 沿用第 2 章的测试信号 figure; imagesc(ts, frq, abs(cfs)); set(gca, 'YDir', 'normal'); colormap('turbo'); colorbar; xlabel('时间 / s'); ylabel('频率 / Hz');cfs是复值小波系数矩阵，行对应频率轴frq，列对应时间轴。画图时要处理三个点：取abs(cfs)看幅值，取平方则容易被误读为能量谱；ts是原始信号时间向量，cfs列数与信号长度相同，所以imagesc可以直接对齐；set(gca, 'YDir', 'normal')把 y 轴翻回数学坐标系，否则频率会自上而下递减。
同样的测试信号，在 spectrogram 里 40 Hz 起点附近的调频线会有轻微发糊，换成 cwt 后低频端的曲线更细、更连续。代价是频率轴不再等间隔，低频段网格更密，图像占用的存储空间也更大。
3.3 尺度、伪频率与 voicesperoctave 怎么配合
新式cwt的第二个返回值frq已经是 Hz 单位，不需要再做尺度换算。老项目里如果用[cfs, scales] = cwt(x, scales)的写法，则需要：
f = scal2frq(scales, 'amor', fs);scal2frq给出的是伪频率，小波有一定带宽，它表示尺度中心对应的频率，不是精确的单频。新旧两种写法在数值上一致，问题出在混用：把新写法的frq当尺度去乘、去画，坐标轴会明显漂移，代码拷进新版本时最容易犯。
voicesperoctave是另一个需要主动设置的参数。它表示每个倍频程内分割的尺度个数，默认值是 10，画图看趋势够用；做脊线提取、特征分类或者要输出较平滑的频率曲线时，建议提到 16 或 24。频率轴细分档位翻倍，计算时间基本线性增长，对一般时长的工程信号可以接受。示例：
[cfs, frq] = cwt(x, fs, 'voicesperoctave', 24);3.4 STFT 与 CWT 的选型对照
| 对比维度 | spectrogram（STFT） | cwt（连续小波） |
|---|---|---|
| 分析窗 | 固定窗长 | 尺度自适应展缩 |
| 低频频率分辨率 | 受窗长限制 | 大尺度下更好 |
| 高频时间定位 | 受窗长限制 | 小尺度下更锐利 |
| 主要参数 | nwin、noverlap、nfft | 小波类型、voicesperoctave |
| 计算开销 | 低 | 中等 |
| 多分量混合信号 | 直接观察，直观 | 受小波旁瓣影响，可能出现干扰 |
工程经验是：先跑一次 spectrogram 粗看全局，再用 cwt 看低频。上面合成信号里的 320 Hz 单频分量，两种方法都能清楚识别，差别主要在调频起始段，频率变化快时 cwt 的线条更连续。
4. 用 EMD 与希尔伯特谱处理强非平稳信号：emd 和 hht 的配合方式
瞬时频率的定义很直观：实信号做 Hilbert 变换后形成解析信号，相位对时间求导就是这个时刻的频率。但这个定义只对单分量信号有意义。实际采集的振动、语音大多是多个分量叠加，直接求瞬时频率会出现负频率、相位跳变一类无物理意义的结果。
EMD 的作用就是把多分量信号拆成若干固有模态函数（IMF）加一个趋势项，之后再对每个 IMF 求瞬时频率。注意，EMD 本身不是时频分析，时频信息产生于后续的 Hilbert 变换；两者在 matlab 里由emd和hht两个函数配合完成。三种路线的定位差异可以这样看：
| 时频分析方法 | 适用信号 | 主要输出 |
|---|---|---|
| STFT | 准平稳或缓变信号 | 幅度矩阵 |
| CWT | 非平稳信号，多分辨率需求 | 小波系数矩阵 |
| EMD + HHT | 强非平稳、模态成分复杂 | 瞬时频率曲线集合 |
4.1 用 emd 把信号拆成 IMF
imf = emd(x, 'Display', 'on'); % 打印每次 sifting 迭代次数imf是一个矩阵，每一行是一个模态分量，最后一行是残差趋势项，行数由信号复杂程度自动决定。拆完之后先看前几个 IMF 是否符合物理直觉：
figure; for k = 1:min(4, size(imf, 1) - 1) subplot(4, 1, k); plot(ts, imf(k, :)); title(sprintf('IMF %d', k)); endDisplay开启后会在分解过程中打印 sifting 的迭代次数，帮助判断是否出现了过分解。如果某个 IMF 的迭代次数异常大，往往意味着信号频率成分过于接近，分解不稳定。
4.2 Sifting 停止条件与模态混叠
emd有两个参数值得关注。SiftMaxIterations控制每个 IMF 的最大筛选迭代次数，默认是 100，工程上遇到持续振荡无收敛迹象时，把它调小反而能避免分解出虚假分量。MaxNumIMF则限制最大模态数量，适合事先知道信号有几类主要成分的场景。
模态混叠是 EMD 最头疼的问题：一个 IMF 里混进了不同时间尺度的成分，或同一分量被拆进相邻两个 IMF。常见对策是集合经验模态分解思想——在信号里多次叠加不同白噪声，分别分解后取平均，让噪声在不同分解中相互抵消。matlab 不内置该函数，自己写一个循环大概几十行。代价是计算量成倍增加，噪声幅度的设定也需要针对信号幅值试。
4.3 hht 输出希尔伯特谱
[hs, fhs, ths] = hht(imf(1:end-1, :), fs, 'FrequencyResolution', 0.5); figure; imagesc(ths, fhs, abs(hs)); axis xy; colormap('turbo'); colorbar; xlabel('时间 / s'); ylabel('频率 / Hz');传入hht的必须是去掉残差的 IMF 矩阵，即imf(1:end-1, :)。残差是趋势项，直接带入会让谱图最底端出现一条贯穿全域的伪频带，干扰后续判断。FrequencyResolution的单位是 Hz，值越小频率轴划分越细，计算时间和内存占用也越大，0.5 作为初值比较合适。
提示：如果
imf只有一行，说明信号没有分解出有效模态，此时hht会报维度错误。先回到 cwt 确认信号是否真的存在明显多分量，再决定要不要继续走 HHT 路线。
5. 验证时频分析结果的三种手段与最常被忽略的 3 个参数坑
时频分析程序跑完，第一件事不是调色板，而是验证结果可不可信。我最常用的验证方法有三个，全部基于已知结构的合成信号。
5.1 用无噪线性 chirp 检验三种方法的频率追踪
fs = 1024; t = 0:1/fs:2-1/fs; x = chirp(t, 40, 2, 220, 'linear'); % 瞬时频率是严格的线性曲线该信号的理论瞬时频率是40 + 90 * t。把 spectrogram、cwt、hht 三种方法的时频峰值画在同一个图里，和这根理论直线对比，偏得越多说明参数越不合适。前面用的带噪二次调频信号适合做定性展示，验证算法正确性时换成这种“答案已知”的信号更高效。
5.2 用 tfridge 自动提取脊线，避免肉眼比色
[s, f, t] = spectrogram(x, hann(256), 128, 1024, fs, 'yaxis'); [fridge, iridge] = tfridge(abs(s).^2, f); plot(t(1:length(fridge)), fridge, 'r', 'LineWidth', 1.5);tfridge的第二个输入必须是 spectrogram 输出的频率向量，fridge返回每个时刻的主脊频率，iridge返回对应索引；对多分量信号加上'NumRidges', 2可以同时提取两条脊线。无论哪种方法，结果里凡是在理论频率附近出现第二条稳定亮线，先检查参数，别急着解释成新物理现象。
5.3 WVD 的交叉项防误读
[wv, f, t] = wvd(x, fs); imagesc(t, f, abs(wv)); axis xy;WVD 聚集性最好，但对多分量信号会产生交叉项——两个真实分量中间出现第三条虚拟分量。工程上遇到多分量时，我一般先 EMD 分解，再对每个 IMF 单独做 WVD 或 Hilbert 谱，这样能避开大部分交叉项。
5.4 三个最容易被忽略的参数坑
第一个是fs漏传。spectrogram、cwt、hht三个函数不传fs时默认按 1 Hz 处理，频率轴上所有数值都会错位，而且这种错位很难通过观察发现。第二个是把nfft当频率分辨率。点数再大也只是补零插值，分辨相邻频率靠加长窗长；两者在代码上差一个参数，在物理上差一个数量级。第三个是hht输入带了残差。残差是趋势项，瞬时频率本身无意义，直接进入hht后会在图上产生一条贯穿全域的伪频带。
整套程序建议封装成统一入参(x, fs, params)、统一出参(time, freq, amp)的三个独立函数，分别对应 STFT、CWT、EMD 三套流程。改成这种结构后，批量跑数据、用 codex 这类工具按格式改参数、甚至把时频结果当作 bilstm 的谱图输入，都只需要接一个接口，不需要再改动核心算法代码。
本文还有配套的精品资源，点击获取
