先聊点实在的短时傅里叶变换STFT在信号处理里几乎是绕不开的标配工具尤其在语音、振动、电力、生物电信号这些领域配合MATLAB做分析和验证算是工程实践里的常见组合。很多朋友最初接触时习惯直接拖一个spectrogram函数出图看起来挺直观但真到了项目里要调参数、改算法、定位异常数据时才发现连“图为什么长这样”都说不清楚。这篇文章不打算只丢一个函数给你。我会从STFT要解决的问题讲起把公式拆开揉碎再给出MATLAB的实现代码和调参经验最后整理我实际踩过的坑。无论你是刚接触时频分析的新手还是想把手写实现和内置函数对齐的老手这篇文章应该都能给你一些参考。1. 先说清楚STFT到底解决了什么问题1.1 傅里叶变换的“盲区”全局频谱丢了时间信息傅里叶变换大家都熟就是把一个时域信号分解成不同频率的正弦波之和。但不知道你有没有想过一个问题做完FFT之后你只得到一个频率谱它告诉你“这个信号里有100Hz的分量、有5kHz的分量”却完全不告诉你这些频率是什么时候出现的。打个比方。假如把信号当成一本书FFT相当于统计整本书里“爱”这个字出现了多少次但你看完统计结果完全不知道这个词是在第几章、第几页出现的。对于平稳信号——比如稳定的正弦波——这个统计结果是足够的因为频率成分从头到尾都没变过。但现实世界的信号几乎都不是平稳的语音里一句话的基频在不停变化轴承故障时的冲击特征只在特定转速阶段出现脑电里某个节律可能只在某个状态下爆发。这就是FFT的“盲区”它是全局分析时间信息被彻底抹掉了。1.2 非平稳信号真实世界里的大多数信号举三个我实际处理过的例子你就能直观感受到这个问题有多普遍。第一个是语音信号。人说话时元音和辅音的频谱差异巨大一句话里音调还在不断起伏。你要是把整段语音直接做FFT得到的频谱是几十个音素频谱的算术平均混成一团根本看不出任何有用的结构。第二个是旋转机械的振动信号。设备启动阶段转速从0逐步上升到额定值此时转频、倍频、边带成分都在随时间移动。如果直接对整个过程做一次FFT得到的是一个“糊掉”的频谱因为不同时刻的频率成分被平均到了一起峰值变得又宽又矮。第三个是电力系统的暂态信号。电压暂降、电弧故障这些现象通常只持续几个毫秒到几个工频周期频率成分是突然出现又突然消失的。常规FFT做出来以后这些瞬态成分会被埋没在稳态50Hz的巨大分量里几乎看不到。这些场景有一个共同点信号的频率成分随时间是变化的我们不但想知道“有哪些频率”还想知道“这些频率在什么时间出现、持续了多久、强度怎么变化”。这正是STFT的用武之地。1.3 STFT的核心思想切段、加窗、逐个分析STFT的思路其实特别朴素既然直接对整个信号做FFT会丢失时间信息那就把信号切成一小段一小段的假设每一小段内部是近似平稳的然后对每一小段分别做FFT最后把结果按时间顺序排列起来。这个过程类似你看书时不满足于全书统计于是决定一页一页地查词频。每一页的统计结果代表“这一页里各种词出现了多少次”按页码顺序排下来你既能看到总体分布也能看出不同章节之间的差异。但这里有个关键细节直接“切段”不行生硬地截断会让每一段信号的边界处产生突变这个突变会在频域造成严重的频谱泄漏。所以正确做法是先给每一段信号乘以一个窗函数让段内信号在边界处平滑地衰减到0再去做FFT。这个“先加窗、再变换”的操作就是STFT的核心流程。2. STFT的数学原理与核心参数拆解2.1 公式拆解从连续到离散STFT的连续形式定义是这样的X(t, f) ∫ x(τ) w(τ - t) e^(-j2πfτ) dτ看着复杂其实含义很直接我们在时间t附近取一块信号用窗函数w(τ-t)把它框出来然后对这一块做傅里叶变换得到的是“时刻t处的频谱”。让t不断滑动就得到一张二维的时频图。工程上信号都是离散采样的所以更常用的是离散形式X(m, k) Σ_{n0}^{L-1} x(n mH) w(n) e^(-j2πkn/L)这里每个符号都有实际含义m是帧序号表示第几帧k是频率索引对应第几个频点L是窗长也就是每帧参与计算的样本点数H是帧移hop size也就是窗每次滑动的样本点数w(n)是窗函数长度为L注意这个式子里最关键的一点我们不是从n0一直积到信号末尾而是从n0积到L-1也就是只看长度为L的一段。这就是“短时”两个字的由来。代码实现时还要加一个参数NFFT也就是做FFT时实际使用的点数。NFFT可以和窗长L相等也可以比L大差别我后面细说。2.2 窗函数为什么要加窗窗怎么选前面说“直接切段会造成频谱泄漏”这里展开讲一下。假设有一段纯正弦信号你恰好截取了整数倍周期的一段FFT结果是一个干净的谱线。但如果你截取的恰好是3.7个周期段首和段尾的相位不连续FFT就会把能量“漏”到旁边的频率上表现为谱线变宽、出现旁瓣这就是频谱泄漏。加窗的作用就是让每段信号的首尾平滑归零消除这种相位不连续。但天下没有免费的午餐不同窗函数在“抑制泄漏”和“保持频率分辨率”之间有不同的取舍。窗函数选择的核心指标有两个主瓣宽度和旁瓣高度。主瓣越窄两个相近频率越容易分开旁瓣越低对弱信号的干扰越小。这两者通常是矛盾的。我用得比较多的是这几种窗窗函数主瓣宽度旁瓣水平适用场景矩形窗最窄最高-13dB频率分辨率优先、信号周期性好汉宁窗较宽较低-31dB通用分析最常用兼顾两者海明窗较宽略低-43dB与汉宁类似旁瓣衰减更平滑布莱克曼窗最宽最低-58dB强干扰背景下分析弱信号实际项目里我绝大多数情况直接用汉宁窗只有在需要精确测量幅值且信号是周期整倍数采样时才会换回矩形窗。至于MATLAB实现hann(L,periodic)返回的就是周期汉宁窗用于谱分析时应该选periodic模式而不是symmetric模式后者用于滤波器设计。2.3 三个关键参数窗长、帧移、NFFTSTFT的参数看似多核心只有三个窗长L、帧移H、FFT点数NFFT。理解了这三个参数你就理解了STFT的绝大部分行为。窗长L决定一帧信号的时间跨度。窗越长一帧里的样本越多频率分辨率越高但时间分辨率越差——因为你无法区分窗长内发生的频率变化。窗越短时间定位越准但频率分辨率越差。帧移H决定相邻两帧的重叠程度。如果H L两帧首尾相接完全不重叠这通常会导致时间轴上的信息不连贯时频图看起来“一条条”的。常见做法是取H L/4或L/2也就是重叠75%或50%。重叠的目的不仅是让图好看更是为了保证每个时域样本都被多个窗覆盖减少信息损失。spectrogram函数默认的重叠就是窗长的75%这个默认值经过了大量工程验证多数情况下是合理的。NFFT是每帧做FFT的点数。如果NFFT L需要在加窗后的数据后面补零。补零不等于提高真实频率分辨率它只是让频谱在频率轴上采样得更密看起来更细腻平滑。真实分辨率由窗长L决定。如果NFFT LMATLAB会直接报错或截断处理这在代码里容易被忽略。2.4 逃不掉的权衡时间分辨率与频率分辨率STFT有一个内在矛盾海森堡测不准原理在信号处理里的体现时间分辨率和频率分辨率的乘积存在下限。窗函数越短你对时间定位越准但频率上就越“看不清”窗函数越长频率越精确但在时间轴上就越“迟钝”。用生活类比来说这就是用显微镜看东西放大倍数越高你能看到的细节越细但视野范围越小。STFT的窗长就是放大倍数频率分辨率是细节时间分辨率是视野。这个矛盾没有完美解只能根据应用场景选择妥协方向分析语音的基频变化需要较好的时间定位窗长通常取20-30ms对应音频采样率下约256-512个点分析机械振动的窄带故障特征需要较高的频率分辨率窗长可以取几百毫秒甚至更长分析电力暂态信号瞬态只有几个毫秒必须用短窗否则瞬态成分会被平均掉参数怎么选归根结底取决于你的信号里“你想看到多快的变化”和“你想分辨多细的频率”之间哪个更重要。3. MATLAB实现函数调用与手写算法两条路3.1 三行代码上手spectrogramMATLAB内置的spectrogram函数封装了完整的STFT流程先看一个最基本的用法fs 8000; % 采样率 8kHz t 0:1/fs:1; x chirp(t, 100, 1, 2000); % 1秒内从100Hz线性扫频到2kHz win hann(256, periodic); % 256点汉宁窗 noverlap 192; % 重叠192点即帧移64点重叠75% nfft 512; % FFT点数 [s, f, t_out] spectrogram(x, win, noverlap, nfft, fs); figure; imagesc(t_out, f, 20*log10(abs(s))); axis xy; xlabel(时间 (s)); ylabel(频率 (Hz)); colorbar;这段代码里返回的s是复数谱矩阵每一列是一帧的FFT结果f是频率轴t_out是每一帧对应的时间中心。imagesc用来显示时频图axis xy是为了让频率轴从低到高正常排列不然MATLAB默认的图片坐标是反的。spectrogram的经典用法还可以简化成一句spectrogram(x, 256, 192, 512, fs, yaxis);不接收返回值时MATLAB会直接画出一张彩色时频图yaxis选项让频率轴竖直朝上。这对快速看一眼数据的特征非常方便。3.2 手写STFT核心代码逐段拆解内置函数好用但有一个问题它是一个黑盒。很多人在实际项目里需要改造算法、嵌入到自己的处理流程中这时候只靠内置函数远远不够。我建议每个做信号处理的人都至少手写一遍STFT写完以后你对加窗、分帧、补零的理解会透彻很多。手写实现分四步分帧、加窗、FFT、拼接。第一步分帧与确定帧数function [X, f, t_out] my_stft(x, win, hop, nfft, fs) L length(win); % 帧长等于窗长 N length(x); % 信号总长度 % 帧数计算确保最后一帧不越界 num_frames floor((N - L) / hop) 1; % 预分配频谱矩阵 X zeros(nfft, num_frames); % 每帧的时间中心 t_out (0:num_frames-1) * hop / fs;注意帧数的计算方式。(N - L) / hop 1的含义是第一帧从样本1开始之后每移动hop个样本产生一帧最后一帧必须完全落在信号范围内。也可以用floor((N-L)/hop)再判断边界情况但上面的公式是工程里最常用的写法。第二步逐帧加窗做FFTfor m 1:num_frames % 取当前帧的样本区间 idx_start (m - 1) * hop 1; idx_end idx_start L - 1; frame x(idx_start:idx_end); % 加窗 frame_windowed frame .* win; % 补零到nfft长度再做FFT X(:, m) fft(frame_windowed, nfft); end这里有两个容易踩的坑。第一窗函数win如果是列向量hann默认返回列向量而frame是从信号里取出来的行向量直接点乘会报维度错误所以要确保两者方向一致。我在代码里统一用win转置成行向量。第二fft(frame_windowed, nfft)这个写法会在长度不足nfft时自动补零如果超过nfft则自动截断这一点用好能省掉手动补零的代码。第三步构建频率轴f (0:nfft-1) * fs / nfft; % 双边谱频率轴 end这个频率轴对应的是完整的双边FFT结果。如果你只想显示单边谱实际工程中更常见需要把频率轴截取到f fs/2对应幅度乘以2直流分量除外。对比一下手写版本返回的X和内置spectrogram返回的s在数值上几乎完全一致只是spectrogram默认做了归一化处理单位不同。比较时需要注意缩放关系。3.3 验证与对比手写结果和内置函数怎么对齐我自己验证手写代码时有一套固定的方法用同一个扫频信号跑内置spectrogram和手写my_stft对比两者的图形和数值。数值对比时注意三个差异点spectrogram默认对窗函数做了能量归一化所以幅值比手写版本的FFT结果小手写版本没做归一化对应的是标准FFT幅值spectrogram返回的f是单边谱的频率向量默认情况下而手写版本是完整双边谱spectrogram的t_out与手写版本可能有半个帧移的偏移取决于对帧中心的定义我比较常用的验证方法是把两边的幅度谱归一化后做差值差值的量级应该在浮点精度级别否则说明某个参数理解有误。4. 参数怎么调不同信号的实战选择4.1 调参顺序先窗型、再窗长、再帧移、再NFFT遇到具体信号时我建议按固定顺序调参而不是想到哪个改哪个。第一步定窗型。默认汉宁窗可以解决95%的问题只有在强干扰下要分离很接近的频率分量时才考虑换成海明窗或布莱克曼窗。第二步定窗长。这一部最难也最影响结果。核心方法是先粗看信号特征如果你要观察的瞬态事件持续大约50ms那窗长应该明显小于这个持续时间比如取20ms左右否则瞬态会被平均掉如果你要区分两个间隔很近的频率窗长得足够长才能分清。这里给一个估算公式频率分辨率约等于fs/L也就是窗长对应的时长取倒数。比如采样率8kHz、窗长256点频率分辨率约为31.25Hz这意味着两个频率相差不到31.25Hz时在时频图上很难分开。第三步定帧移。帧移决定重叠率一般取窗长的1/4到1/2。重叠越多时频图越平滑但计算量越大重叠太少帧与帧之间跳变明显图像会看到“指纹状”条纹。工程上最常用的折中是75%重叠也就是帧移为窗长的1/4。第四步定NFFT。NFFT主要影响频率轴上的点数一般取大于等于窗长的幂次即可。NFFT越大频率曲线的显示越细腻但并不会带来真正的分辨率提升。推荐取NFFT等于下一个比窗长大的2的幂次比如窗长256直接取256窗长300可以取512。4.2 三个具体场景的参数量级参考我整理了一张表是几个典型信号场景下的起始参考参数。注意这只是“起步点”实际要根据你的信号特征微调信号类型采样率窗长帧移NFFT说明语音分析16kHz512点32ms128点8ms1024兼顾基频和音节边界轴承振动48kHz4096点85ms1024点21ms8192需要分辨边带频率电力暂态12.8kHz128点10ms32点2.5ms256快速捕捉暂态过程脑电250Hz128点0.51s64点0.256s256分辨alpha/beta节律语音那个例子多说一句32ms窗长对应的频率分辨率约31Hz这正好能区分开相邻的谐波分量又不至于把音节边界抹掉是语音分析的经典选择。4.3 频谱泄露、幅值校准与边界失真的处理参数调好之后还会有三个绕不开的事项要处理。频谱泄露的残余问题。加窗能大幅抑制频谱泄漏但不能完全消除。如果你观察到的频谱在某些频率点上有异常高的“底座”可以先检查窗函数有没有加错。一个常见乌龙是窗函数和信号方向不一致导致加窗失败结果相当于用了矩形窗旁瓣高得离谱。幅度校准。你如果要用STFT做精确幅值分析比如振动信号的包络谱提取一定要注意归一化问题。手写FFT版本直接得到的幅值是真实幅值的一半对于正弦分量能量分散在正负频率上但加了周期汉宁窗之后幅值又会被压缩约一半。实际使用时我通常在显示之前乘一个校准系数把加窗造成的衰减补回来。这个系数可以这样算对窗函数求平均增益mean(win)然后用2 / mean(win)作为幅度恢复系数。我用这个方法做过整套幅值校准误差可以控制在2%以内。边界失真。时频图最左边和最右边的几帧因为窗超出信号边界数据不完整计算出的谱往往不可靠。这不算bug是STFT的固有特性。处理方式有两种一是直接把边界区域的时频数据标记为不可信分析时跳过二是对信号做反射延拓把信号两端镜像扩展一部分再计算STFT这样边界效应会被推到延拓区域里感兴趣的有效信号段就不受影响了。5. 常见问题排查速查表5.1 一张表看懂常见问题与解决方法这几年我在社区和实际项目里帮人看过不少STFT相关问题整理成一个速查表遇到问题直接对照排查现象可能原因解决办法时频图一片模糊看不出频率峰窗太长或重叠太少缩短窗长提高重叠率高频段出现明显的横向亮带加窗失败等效矩形窗检查窗向量维度和信号维度是否一致频率轴显示为0到8kHz但信号只有4kHz显示的是双边谱截取单边谱显示时间轴和预期完全对不上没有正确传入fs或帧移单位混淆确认fs参数用秒为单位计算帧移图像左右边缘有两道异常亮区边界效应反射延拓或忽略边界帧图像有规律的“条纹”重叠率太低提高重叠到75%外行的正弦幅值明显偏小没有做窗增益校准用2/mean(win)校准幅值频谱上出现“假峰”信号中有工频或采样率整倍数干扰先做滤波预处理再STFT5.2 单独说几个绕不开的坑有几个坑虽然不是STFT本身的问题但在MATLAB环境里特别常见值得单独提醒。第一个是图像坐标轴方向。直接image(20*log10(abs(s)))显示时图默认y轴是倒置的频率看起来是“从上往下”增长的。很多人第一次出图都会懵一下怀疑自己算法写错了。解决方案就是axis xy或者set(gca,YDir,normal)顺手也把坐标刻度改成频率值而非点数。第二个是中文注释和显示乱码。新版MATLAB在较新版本中默认用UTF-8编码但旧版本脚本用GBK保存互相打开会乱码。这种问题在团队协作、换电脑继续写代码时特别常见。我现在养成的习惯是脚本里涉及中文注释时统一在文件开头加一句说明或者干脆用英文注释省得在不同版本之间来回折腾。如果你拿到一个旧版本写的脚本在自己新版里打开文件菜单里选“另存为”把编码改成UTF-8通常能解决。第三个是spectrogram和imagesc显示的幅值尺度问题。spectrogram直接调用不带输出参数时画的是对数幅值谱但你拿返回值自己用imagesc显示时默认是线性幅值视觉效果差异巨大会造成误读。我建议把自己绘图时的显示统一改成20*log10(abs(s) eps)加一个小的常数防止log(0)产生-Inf。我的一些实际体会写到这里STFT从原理到代码到调参到排错就算完整过了一遍。我个人实际工作中的工作流是先用内置函数快速出图看整体特征建立对信号的基本感觉然后根据分析需求更关注时间还是更关注频率确定窗长量级手写实现做定制化处理最后用幅值校准后的结果跟内置函数对比确认两者没有系统性偏差。这套流程用了好几年效率和准确度都比较稳定。最后分享一个小技巧调参时别只看最终图把X矩阵的数值也打印出来看几行。很多参数问题在图上看差异并不明显但数值分布能一眼暴露异常——比如某几帧的幅度特别大通常就是边界帧比如幅度整体偏小那就是校准没做。数据是检查代码逻辑最可靠的依据这是每次调试时我都会提醒自己的事。