简介面向非平稳信号时频分析需求的STFT MATLAB实现适合信号处理学习者、工程师与科研人员快速上手。资源以短时傅里叶变换为核心通过加窗分段处理克服传统傅里叶变换对时变信号分析的局限可应用于声学、生物医学及电力系统瞬态信号等场景。压缩包共3个文件均为.m脚本整体仅1KB代码精简便于用户按需调整窗函数类型、窗口长度与重叠比例并结合fft绘制时频图。资源提供直接可运行的MATLAB脚本结构清晰、参数集中便于二次开发与教学演示。已有261人学习下载适合希望理解STFT原理并直接运行验证的入门及进阶用户。通过运行脚本可直观观察非平稳信号频谱随时间的变化规律为后续小波分析、匹配追踪等时频方法的学习奠定基础。1. 一张图丢掉时间轴之后非平稳信号只能靠STFT找回来看到stft.rar这种命名大部分情况是从课程网站或代码分享站拖下来的整套 MATLAB 脚本。里面的核心几乎都是同一个东西把一段信号切成小段逐段做傅里叶变换再拼成一张横轴是时间、纵轴是频率的图。这套流程本身不复杂但它是非平稳信号分析里最绕不开的起点。旋转机械的振动、语音的基频漂移、PWM 波的占空比突变、高速信号接口上的瞬态毛刺——这些信号的频率成分随时间变化一次 FFT 只能给你一张全局频谱时间信息被平均掉了8QAM 或 64QAM 的载波漂移、DTMB 信号里的突发干扰在这种图里全看不出门道。STFT 做的就是退一步看不完整段就分段看用时间分辨率换回频率随时间的演化轨迹。这篇按“从跑通到调参再到提取特征”的顺序把 STFT 在 MATLAB 里的做法讲透适合刚接触时频分析的人照步骤复现也适合已经画过频谱图但对参数取舍不笃定的熟手。2. MATLAB 里跑通 STFT 的三种写法spectrogram 命令、手写 FFT 循环、滤波器组视角STFT 在 MATLAB 里最少只需要一行命令但只调用命令不搞清内部在做什么遇到参数不知道怎么改、图形不知道怎么解释。这一章先给最直接的跑通路径再给出一个手写版本用来对照原理最后补充滤波器组视角帮你在面对泄漏、带宽这些概念时不发懵。2.1 spectrogram 一句话出图先看时频全貌先造一段非平稳信号模拟一个从 100 Hz 线性扫到 800 Hz 的啁啾中间再叠一个短时突发分量然后用spectrogram直接画频谱图。fs 4000; % 采样率 4 kHz t (0:1/fs:2-1/fs); % 2 秒时长 f0 100; f1 800; x chirp(t, f0, 2, f1); % 线性啁啾 100→800 Hz % 位置 1.2 s 附近叠一个 300 Hz 短时正弦0.2 s 时长 idx round(1.2*fs):round(1.4*fs); x(idx) x(idx) 1.5*sin(2*pi*300*t(idx)); % spectrogram: 窗长 256 点, 重叠 75%, FFT 点数 512 [ S, F, T ] spectrogram(x, hann(256), 192, 512, fs); imagesc(T, F, 20*log10(abs(S)eps)); axis xy; xlabel(时间 (s)); ylabel(频率 (Hz)); colorbar;spectrogram最核心的四个参数是窗、重叠、FFT 点数、采样率。hann(256)表示一次只看 256 个采样点对应 64 ms这个窗长决定了“这一段内信号近似平稳”的时间尺度192是相邻两窗的重叠采样数256 减 192 等于 64表示每 64 个点滑动一次窗对应时间步进 16 ms重叠的目的只是让时间轴画出来更连续不增加频率分辨率512是 FFT 长度比窗长多出来的部分做零填充得到更密的频率网格。图像用20*log10(abs(S))转成 dB 显示不然幅度小的频率分量会被强的完全压掉。注意这里colormap默认是parula深色代表能量低。你会看到一条从左下角向右上角延伸的亮线那就是 100 Hz 到 800 Hz 的扫频轨迹在 1.2 秒位置有一条水平短线是叠加的 300 Hz 突发分量。两个特征在同一张图里清晰分离这就是时频分析相对普通频谱的核心优势。2.2 手写 STFT把原理变成你随时能改的代码spectrogram封装得太多初学者容易把它当黑盒。手写一个短循环每一步都对应一个明确操作后续做脊线提取、自定义时频重排时就能直接改这块代码。function [ S, F, T ] my_stft(x, win, hop, nfft, fs) % x : 输入信号 % win : 窗函数序列如 hann(256) % hop : 帧移即每次跳过的采样点数 % nfft : FFT 点数通常 length(win) % fs : 采样率 Nx length(x); Nw length(win); nf floor((Nx - Nw) / hop) 1; % 总帧数 S zeros(nfft, nf); for n 1:nf idx (n-1)*hop 1 : (n-1)*hop Nw; seg x(idx) .* win; % 加窗抑制边缘突变 S(:,n) fft(seg, nfft); % 逐帧 FFT列对应时间帧 end F (0:nfft/2-1) * fs / nfft; % 频率轴取单边 S S(1:nfft/2, :); % 去掉镜像频率 T (0:nf-1) * hop / fs; % 时间轴对应每帧起始时刻这个函数里去均值和高通滤波都没做实际数据先x x - mean(x)更稳妥。x(idx)取当前帧的采样点乘窗后直接fft第 n 列存的是一整帧的复数频谱。对加窗的帧做 FFT 时由于窗的旁瓣泄漏频谱的每根谱线都不再是纯粹的单一频率而是主瓣叠加旁瓣的形态这一点到第 3 章详细展开。帧移hop在这里是显式参数它和spectrogram里的noverlap的关系是noverlap Nw - hop。如果想手动控制时间分辨率改hop比改noverlap更直觉。nfft取不取Nw的整数次幂影响的是计算效率和频率网格密度不影响真实频率分辨率真实分辨率只由窗长决定。2.3 滤波器组视角每条频率线都是一次带通滤波除了逐帧切片STFT 还有另一个等价的解释把信号通过一组中心频率不同的带通滤波器每个滤波器的输出幅度包络就是时频图上对应频率行的强度变化。这个视角下窗函数其实就是滤波器的冲击响应主瓣宽度决定滤波器的通带宽度旁瓣衰减决定对相邻频率的抑制能力。这个视角在排查频谱图上的“横纹”时非常有用。当你发现整个时频图出现平行的固定频率亮线即 50 Hz 工频及其谐波时这不是 STFT 方法的问题而是滤波器的旁瓣过高把强干扰泄漏到了邻近频率换旁瓣衰减更大的窗或先做陷波滤波比调窗口长度更对症。理解这个等价关系后后续遇到泄漏、频率打架时你不会下意识地乱调参数。3. 窗长、重叠率与 FFT 点数时频分辨率的三角博弈STFT 的成像质量几乎全部由窗长、重叠率和 FFT 点数决定但三者的作用常被混淆。窗长改变的是真实时频分辨率重叠率只影响时间轴的平滑度FFT 点数只影响频率网格密度。分开理解比死记公式有用得多。3.1 窗长定的是“信号段近似平稳”的时间尺度窗长的选择是时频分析里第一个要做的决定。窗越长每帧包含的采样点越多频率分辨率越高——两个靠近的频率分量能在频谱图上分开窗越短每帧的时间跨度越短频率随时间的变化越不容易被抹平。这组矛盾是海森堡不确定原理在信号处理中的体现时域窗越窄频域主瓣越宽频率分辨率越差时域窗越宽频域主瓣越窄但在时间轴上变化剧烈的频率会被一帧平均掉。举例来说一个信号在 100 ms 内从 400 Hz 跳到 410 Hz。用 256 点窗在 2 kHz 采样率下看窗长 128 ms一次 FFT 同时覆盖了频率跳变前后的两个分量频谱图上会同时出现 400 Hz 和 410 Hz 两条亮线看起来像两个频率同时存在。但如果用 32 点窗窗长 16 ms跳变前后被分到不同帧就能看到亮线在时间轴上从 400 Hz 移到 410 Hz。反过来两个相距 5 Hz 的频率分量用短窗看因为频率分辨率只有 1/16ms ≈ 62.5 Hz两条亮线直接糊成一条。具体拿信号完整性测试举例看高速信号接口上的反射毛刺毛刺持续时间通常只有几纳秒窗长必须短到能包含毛刺但不长于一整段码流通常选码元周期的 2 到 4 倍而分析 PWM 信号的占空比渐变占空比变化率很慢可以把窗拉长到几十个周期来获得更干净的频率分辨。3.2 重叠率平滑时间轴但与频率分辨率无关重叠率是初学者最容易调的参数因为调它时图像变化立竿见影。提高重叠率后时间轴上的帧间距变小频率分量在帧间的幅度和频率变化变得更连续但本质只是插值——没有引入任何新的频率信息。spectrogram默认是 75% 重叠即窗长 256 时noverlap为 192。为什么常用 75%当窗函数是汉宁窗时相邻两窗以 75% 重叠叠加窗函数的平方累加近似一条平坦的直线。这意味着每一帧在拼回时频图时能量贡献基本均匀不会被窗函数本身调制出周期性的明暗条纹。低于 50% 重叠时两个窗之间的区域能量贡献明显减弱时频图会出现垂直于时间轴的浅色条纹。我的习惯是第一遍先用 75% 重叠看全貌确认时频结构后要精确提取瞬时频率时改到 87.5% 或 90% 获得更平滑的脊线但不要在重叠率上花太多时间——如果图模糊先回去改窗长。3.3 FFT 补零只是细化网格不带来真实分辨率spectrogram里的nfft参数常被误解为“提高分辨率的手段”。实际上当nfft Nw时多出来的部分是对信号补零补零只是在 FFT 输出时插入额外的频率网格点让谱线看起来更密、峰值位置读起来更精确但两个相邻频率分量能不能分开仍然由窗长决定。一个很直观的验证对 400 Hz 和 405 Hz 两个正弦的叠加信号用 256 点窗频率分辨率约 1/64ms ≈ 15.6 Hz两个分量距离只有 5 Hz无论nfft取 1024 还是 4096图上都是一个宽峰改用 1024 点窗频率分辨率变成约 3.9 Hz即使nfft仍取 1024两个峰会明显分离开。所以调nfft前先想清楚是网格不够密还是真实分辨率不够3.4 窗函数选型与泄漏用一张表结束纠结常见的窗函数按旁瓣衰减和主瓣宽度分几类选择时就是在这两个指标间取舍。旁瓣衰减大则泄漏少但主瓣宽度大则频率选择性差主瓣窄则但旁瓣衰减有限。选窗时先确认你的首要需求是幅度精度还是频率分离度。窗函数主瓣宽度旁瓣衰减典型用途rect最窄-13 dB瞬态检测优先时间分辨率hann较窄-31 dB通用默认兼顾旁瓣抑制与带宽hamming与 hann 接近-43 dB对幅度精度要求略高的测量blackman较宽-58 dB强干扰下需要大幅压低旁瓣时kaiser参数可调可自定义需要精确控制主瓣宽度时工程上最常用的组合是 hann 窗加 75% 重叠。hamming 和 hann 的差异在于前者降低了第一个旁瓣、但远端旁瓣衰减反不及 hann所以如果关注的是离主瓣很远的干扰hann 反而更干净。blackman 适合在频谱动态范围很大时使用比如监控信号里混有强窄带干扰需要看弱分量时用 blackman 可以压低干扰的泄漏底噪代价是主瓣更宽短时里两个频率很近的分量不容易分开。rect 窗几乎只在看瞬态冲击时用因为它完全不抑制边缘突变时域分辨最好但频谱泄漏最严重。提示如果不确定当前配置的实际分辨率可用2*fs/Nw估算 3 dB 带宽0.89*fs/Nw估算主瓣宽度这两个数是你解释时频图时判断“两条线能不能分开”的依据。4. 拿到 stft.rar 旧代码时函数差异、单位换算与常见报错从stft.rar这类资源包里解压出来的代码大多是多年前的版本第一次运行通常会踩几个固定的坑。这一章把最常见的环境适配、单位换算和报错整理成直接能对照的清单。4.1 specgram 已经废弃新代码一律用 spectrogram老版本 MATLAB 里的specgram函数在 R2016a 之前的帮助文档里还能找到之后的版本不再推荐R2017b 以后直接移除。作为替代spectrogram的参数顺序完全不同specgram(x, nfft, fs, win)而spectrogram是spectrogram(x, win, noverlap, nfft, fs)。解压旧代码后第一步是全局搜索specgram把调用顺序换成spectrogram的规范格式。如果旧代码里用了specgram(x, 256, fs, hanning(256))含义是 FFT 点数 256、采样率 fs、窗函数 256 点汉宁窗。对应的新写法是spectrogram(x, hann(256), 192, 256, fs)重叠率 75% 为默认推荐。4.2 坐标轴朝向、行列顺序与单位换算spectrogram在不同调用方式下返回结果的行列语义不同这也是旧代码最常见的“图看起来不对”的原因。当返回三个输出[S, F, T]时S 是复数矩阵行对应频率、列对应时间帧当只有两个输出[S, F, T]配合imagesc画图时没问题但如果直接用plot(abs(S))或者把 S 存成 CSV 再处理行列颠倒会直接导致频率轴错误。spectrogram的imagesc显示还有个微妙的细节spectrogram(x, ..., yaxis)默认认为S的列是频率与imagesc(T, F, S)的行列语义不同直接用imagesc(T, F, S)会得到一个转置的图。我的习惯是始终用[S, F, T] spectrogram(...)三输出格式然后imagesc(T, F, abs(S))再axis xy翻转 y 轴方向这样行列语义最明确。幅度单位也有坑。spectrogram返回的 S 是复数谱幅度不是功率谱密度20*log10(abs(S))得到的是相对 dB不是 dBm。如果需要绝对功率要另做校准乘以窗的等效噪声带宽再转 dBm。旧代码里最常见的错误是直接10*log10(abs(S))——因为 S 是幅度要以功率显示需要20*log10混用之后图里所有峰值相对关系都会失真。4.3 三个高频报错与对应改法第一类Error using spectrogram. The input signal must be a vector.多发生在从文件直接读入数据后忘记取单通道或者读入的是按列存储的矩阵。改法x x(:,1)或x x(:)。第二类窗长大于信号长度。当信号只有几百个采样点却用了 1024 点的窗spectrogram直接报错或输出空矩阵。改法Nw min(256, length(x))或者先对信号做截断。对于短信号窗长过长不仅报错即使勉强画出来也会让时间分辨率低到没有任何分析价值。第三类采样率参数与时间轴不匹配。spectrogram(x, win, noverlap, nfft, fs)中的fs只影响F和T的刻度值不影响 S 的数值。如果 fs 填错图上的频率轴会成比例偏移。比如 4 kHz 采样率却填了 400横轴频率数值全部变成原来的十分之一但形状看起来完全正常。提示遇到这类问题时先用已知信号验证——扫频信号画出来后频率轴的起点和终点是否与你设置的一致这是最快定位单位换算错误的方法。5. 从频谱图反推瞬时频率最大能量脊线提取画出一张漂亮的时频图不算结束工程上更常需要的是从图中自动读出频率随时间的变化曲线比如变频器的输出频率斜坡、PWM 载波的频率漂移、高速信号接口上时钟的瞬态偏差。STFT 的时频矩阵本质上是一张二维能量图每个时间点上能量最大的频率位置就是该时刻的瞬时频率估计这就是脊线提取。5.1 脊线是时间轴上“能量最大的频率线”对于单分量信号比如纯调频的啁啾时频矩阵每一列的最大值对应的频率就是该帧信号的频率。对于多分量信号比如同时存在基波和强谐波的振动信号直接取最大值只会追踪到能量最强的分量。这时有两种做法限定频率搜索范围比如只搜 20 Hz 到 200 Hz 的基波频带或者对时频矩阵做局部峰值查找把每个时间点的多个峰值都找出来再按频带分组。5.2 MATLAB 提取与叠加代码% 沿用第 2 章的 x 信号做 STFT fs 4000; [ S, F, T ] spectrogram(x, hann(256), 192, 512, fs); P abs(S); % 幅度谱 [ ~, idx ] max(P, [], 1); % 每个时间帧取最大幅度的频率索引 f_ridge F(idx); % 映射到实际频率值 figure; imagesc(T, F, 20*log10(P eps)); axis xy; colormap(jet); hold on; plot(T, f_ridge, w-, LineWidth, 1.5); hold off; xlabel(时间 (s)); ylabel(频率 (Hz)); title(STFT 频谱图与最大能量脊线);max(P, [], 1)的[]表示在第一个维度频率上取最大值1表示沿时间方向逐列输出。F(idx)把索引转成频率。这一步如果 S 是复数矩阵直接取abs(S)再找最大值因为最大值对应的可能是负频率镜像所以先取单边频谱即S(1:nfft/2, :)才能保证F(idx)对应正确频率。脊线叠加到图上后如果信号是线性啁啾白线应该与亮带的中心位置重合。检核方法对 100 Hz 到 800 Hz、2 秒的啁啾理论瞬时频率是f(t) 100 350*t脊线拟合的斜率应该接近 350 Hz/s偏差在 5% 以内算正常。如果偏差大大概率是窗长太长导致跳变被平滑或是噪声干扰下最大值被噪声峰值抢占。5.3 跳变信号下的脊线断层与改进对非平稳度更高的信号比如 PWM 信号占空比突变或跳频通信信号直接取最大值会把跳变瞬间的脊线切成一截一截频率在突变处不连续。改进做法是加一个“频率连续性约束”的跟踪算法每个时间点不只取全局最大而是落在上一帧脊线频率邻域内找最大峰值邻域宽度一般设为频率分辨率的 2 到 3 倍。这种跟踪在处理强噪声场景时比逐帧独立取最大鲁棒得多。对多分量信号脊线提取前要先用imagesc看清有多少条亮带再按频带分别搜索。常用的替代实现是tfridge函数它支持指定脊线数量但返回的脊线质量受时频分辨率影响很大窗长设得不好时同样会跟踪到噪声上。拿 8QAM 信号做解调前的载波频偏估计时先用 75% 重叠和较大 NFFT 提取到平滑脊线再对脊线做滑动平均得到的频偏估计精度足以支撑后续同步环路工作。脊线提取是 STFT 从“看图”到“量化”的转折点也是把时频分析从仪器的叠加测量里真正分离出来的关键一步。本文还有配套的精品资源点击获取