简介这份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 的频率换算与 voicesperoctaveSTFT 的窗长固定是它的根本矛盾。低频分量周期长需要长窗才能积累足够周期数否则频率分辨率很差高频瞬态又需要短窗才能定位到具体时刻。固定窗长无论怎么选都只能在一个频段上表现良好。连续小波变换用“尺度”代替“窗长”。尺度小时小波在时间上被压窄适合抓高频细节尺度大时小波在时间上展宽天然适合作低频分析。matlab 的cwt函数封装了这套过程输出直接是“时间-频率”二维结果而不需要像老代码那样先算尺度再手工换算。3.1 为什么 STFT 在低频段看不清用第 2 章的参数计算一下fs 1024nwin 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 的选型对照对比维度spectrogramSTFTcwt连续小波分析窗固定窗长尺度自适应展缩低频频率分辨率受窗长限制大尺度下更好高频时间定位受窗长限制小尺度下更锐利主要参数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 把信号拆成 IMFimf 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 的谱图输入都只需要接一个接口不需要再改动核心算法代码。本文还有配套的精品资源点击获取