做信号处理的人早晚要碰时频分析。STFT短时傅里叶变换是其中最基础、最常用也最适合入门的一种。我刚工作那会儿处理齿轮箱振动信号FFT频谱上明明有一堆峰值但设备处于启停变速工况频率一直在漂整段信号做傅里叶变换等于把时间信息全平均掉了怎么看都对应不上实际故障时刻。后来换用STFT把时变频率拉到时间-频率平面上一看问题瞬间清楚。这篇文章就从原理讲起再用Matlab从底层手写一套STFT代码不调用spectrogram这些现成API顺带聊聊窗长、重叠、补零这些绕不开的坑。适合刚开始学信号处理、做故障诊断或者语音分析的朋友想搞明白时频图到底怎么来的照着代码敲一遍就会了。1. 为什么需要时频分析从傅里叶变换的两大局限说起1.1 傅里叶变换把时间轴“抹掉”了先看传统傅里叶变换在干什么。对连续信号(x(t))傅里叶变换是[ X(f)\int_{-\infty}^{\infty}x(t)e^{-j2\pi ft}dt ]这个积分把时间变量(t)从负无穷积到正无穷输出(X(f))只和频率有关和时刻无关。写成离散形式也一样做FFT时默认把整段信号当成一个周期信号来分解。只要信号是平稳的也就是统计特性不随时间变化这种做法没问题。可工程上遇到的信号大量是非平稳的语音里的元音和辅音、雷达回波里的多普勒突变、旋转机械启停过程中的转频爬升这些信号的频率成分是随时间变化的。整段做FFT得到的结果是所有时刻频谱的叠加平均。它告诉你“这个频段有能量”却没办法告诉你“这个频率是出现在第几秒”。我做齿轮箱诊断时遇到的情况特别典型。齿轮箱在升速阶段齿面出现点蚀故障特征频率是转频的整数倍但转频从10Hz一路升到50Hz对应的特征频率也跟着变。整段FFT把这种线性调频成分变成了一段宽宽的隆起峰值扁平根本无法和理论值对齐。而实际上故障是在特定转速区间才激发出明显冲击的我想要的是“哪个时刻、转速到多少、这个故障频率出现了”传统傅里叶变换做不到。1.2 非平稳信号的时频联合需求为了解决这个问题最简单的思路就是既然全局傅里叶变换丢失时间信息那把信号切成很多小段对每一小段分别做傅里叶变换再把结果拼起来。这样一来横轴有了时间纵轴是频率颜色深浅表示能量大小就得到一张二维的时频图。这就是“短时”两个字的核心含义只在一小段时间窗口内做频谱分析窗口移动覆盖整个信号。时频分析领域其实有不止一种工具STFT、小波变换、Wigner-Ville分布、Hilbert-Huang变换各有利弊。STFT最大的优点是直观、线性、快速没有交叉项的干扰算法最成熟工程上用得也最多。虽然它受制于窗函数时间分辨率和频率分辨率不能兼得但只要你理解了这个限制绝大多数应用场景下STFT都够用了。这篇文章完整实现STFT就是给后面理解小波、搞懂谱图打基础的。2. STFT核心思想把长信号切成小段再傅里叶2.1 从连续定义到离散公式短时傅里叶变换的连续形式是[ X(t,f)\int_{-\infty}^{\infty}x(\tau)w(\tau-t)e^{-j2\pi f\tau}d\tau ](w(\tau-t))是一个以当前时刻(t)为中心、长度为有限支撑的窗函数它在(t)附近取1远离开来逐渐衰减到0把信号截断成一小段。把这个窗沿着时间轴滑动每滑动一个位置就对窗内信号做一次傅里叶变换。所以STFT本质上是“加了窗的傅里叶变换族”。实际计算机处理的是离散信号。若采样率为(f_s)信号为(x[n])窗长为(N)帧移为(H)FFT点数为(N_{FFT})那么离散STFT可以写成[ X[m,k]\sum_{n0}^{N-1}x[nmH]w[n]e^{-j2\pi kn/N_{FFT}} ]这里(m)是帧编号对应时间轴第(m)个窗(k)是频率序号对应频率(k f_s / N_{FFT})。把(m)从0数到(M-1)就得到一列频谱向量最后堆成一个大小为((N_{FFT}/21)\times M)的复数矩阵。这个矩阵就是STFT的核心结果通常我们只看它的幅值再取对数转成dB画成伪彩图。要注意的是(N_{FFT})和窗长(N)不一定要相等。窗长决定实际参与加权的采样点数(N_{FFT})决定了FFT的运算长度。(N_{FFT})大于窗长时本质上是在给截断后的信号补零补零能让频谱曲线更平滑但不会提高真实的物理分辨率。后面我会专门讲这个。2.2 窗函数怎么选汉宁窗、汉明窗、矩形窗在STFT里窗函数的作用不只是“截断”还会影响频谱的质量。直接截断相当于乘矩形窗矩形窗在频域的主瓣很窄但旁瓣很高是-13dB左右会让频谱产生严重的泄漏强频率成分旁边出现一堆假峰。窗函数的核心指标是主瓣宽度和旁瓣衰减主瓣越窄分辨相近频率的能力越强旁瓣越低抑制频谱泄漏的能力越强。两者很难兼得。常见窗函数大致如下窗类型主瓣宽度归一化最大旁瓣幅度适用场景矩形窗最窄-13dB瞬态信号同时定位时间起点汉宁窗中等-31dB通用分析平衡分辨率和泄漏我用得最多汉明窗中等-43dB语音分析旁瓣衰减略好于汉宁布莱克曼窗较宽-58dB需要强抑制泄漏但不关心时间定位时高斯窗可调可调时间-频率分辨率可调适合非平稳生活化一点理解矩形窗像拿一把剪刀硬生生剪断信号断点处的突变在频域里就是一串高频旁瓣汉宁窗像用一块软布在信号边缘逐步过渡让截断变得平滑频谱自然也干净得多。实际做STFT我一般默认用汉宁窗只有在分析雷达脉冲这类本就短暂、需要精确时间起点的信号时才考虑矩形窗或更短的高斯窗。2.3 窗长、步长、补零与分辨率的关系STFT分辨率的核心约束来自海森堡-盖博不等式时间分辨率(\Delta t)和频率分辨率(\Delta f)的乘积存在下限不可能同时无限提高。窗长越短时间定位越准但频率分辨率越差窗长越长频率分辨越细但时间定位越模糊。这就像用相机拍奔跑的人快门时间短能凝固瞬间但画面可能糊到看不清细节快门时间长能看清轮廓但人物拖影严重。具体到离散实现如果不做补零频率分辨率约等于[ \Delta f \approx \frac{f_s}{N} ]其中(N)是窗长。采样率1kHz时窗长100点对应10Hz分辨率窗长1000点对应1Hz。如果想让频谱上两个相差5Hz的频率成分分开窗长至少得让(\Delta f)小于5Hz也就是(N f_s / 5)。时间分辨率则由窗长和帧移共同决定窗长约等于“看这一小段信号用了多长的时间”帧移(H)决定了时间轴上的采样间隔时间点约为每(H/f_s)秒一个。帧移小比如重叠75%能让时频图更平滑但计算量也更大帧移大比如完全不重叠会漏掉信号变化细节时频图容易出现格子感。补零的情况要单独说。很多人以为把(N_{FFT})设得很大分辨率就变高了这是个误区。补零只是在FFT的末尾加零改变的是频谱抽样的密度让FFT结果看起来曲线更细致但两个相邻频率峰值的真实可分辨性没有提升因为物理分辨由窗长决定。举个例子100点窗补零到1024点频谱上的“点数”多了但如果信号里有两个相隔5Hz的成分100点窗照样分不开。补零的好处是能让谱峰位置更精确在绘图时更美观所以实操中我常把(N_{FFT})设为大于窗长的最近的2的幂。3. 手动实现STFT一行行写清楚不调Matlab API3.1 构造测试信号线性调频信号加瞬态冲击写代码之前先准备一个测试信号。最好的测试信号是线性调频chirp信号它的频率随时间线性变化在时频图上应该是一条清晰斜线方便直观判断STFT是否正确。我不调用chirp函数直接用正弦表达式生成fs 1000; % 采样率 1000 Hz t_total 2; % 总时长 2 秒 t 0:1/fs:t_total-1/fs; % 时间序列 N_total length(t); % 频率从 50 Hz 线性增长到 300 Hz f0 50; f1 300; phase 2*pi*(f0*t (f1-f0)/(2*t_total)*t.^2); x_chirp sin(phase); % 在 0.8 秒处叠加一个 500 Hz 的瞬态冲击 impulse_t0 0.8; impulse_tau 0.02; x_impulse 0.6 * exp(-((t-impulse_t0).^2)/(2*impulse_tau^2)) .* sin(2*pi*500*t); x x_chirp x_impulse;这段代码不依赖任何工具箱信号里有两个关键特征一是50~300Hz的持续调频斜线二是0.8秒处的一个短暂高频冲击。STFT做出来后应该能在时频图上同时看到斜线和冲击对应的竖直亮条。如果代码写错了这两个特征会以各种诡异的方式变形。3.2 核心函数编写分帧、加窗、FFT、存储现在写STFT的核心函数。思路很简单把信号按窗长分帧每帧乘窗做FFT把结果存到矩阵列里。为了避免调用Matlab的spectrogram或stft我们全程用循环实现逻辑一目了然。function [S, f, t_center] my_stft(x, win, hop, nfft, fs) % 手动实现短时傅里叶变换 % 输入: % x - 输入信号行向量 % win - 窗函数向量长度等于窗长 % hop - 帧移相邻窗起点间隔的采样点数 % nfft - FFT点数建议 窗长 % fs - 采样率用于生成坐标 % 输出: % S - 单边幅值谱矩阵大小为 (nfft/21) x numFrames % f - 频率坐标向量 % t_center- 每帧对应的时间中心坐标 if size(x,1) size(x,2) x x.; % 统一转成行向量 end xlen length(x); winLen length(win); step hop; % 计算帧数丢弃末尾不足一窗的零头 numFrames floor((xlen - winLen) / step) 1; S zeros(floor(nfft/2) 1, numFrames); for m 1:numFrames startIdx (m-1) * step 1; idx startIdx : startIdx winLen - 1; segment x(idx) .* win; % 加窗 spec fft(segment, nfft); % nfft点FFT if mod(nfft, 2) 0 S(:, m) spec(1:nfft/2 1); % 单边谱 else S(:, m) spec(1:(nfft1)/2); % nfft为奇数的情况 end end % 频率坐标单位 Hz f (0:size(S,1)-1) * fs / nfft; % 时间坐标用窗中心时刻单位 s % 第m帧起点是 (m-1)*step中心是起点 winLen/2 t_center ((0:numFrames-1) * step winLen/2) / fs; end有几个细节必须注意。第一分帧时用了floor信号末尾不足一个窗长的部分直接丢弃这是最常见的做法避免边界补零造成虚假能量。第二取单边谱时对于偶数nfft频点包括直流和奈奎斯特频率所以长度是nfft/2 1对于奇数我单独处理了一下不过实际使用中通常设成2的幂很少遇到奇数。第三时间坐标用窗中心而不是窗起点这样每一帧的频谱就对应窗中心那个时刻画图时不会向右偏移半个窗长这个细节很多人容易忽略。3.3 画时频图从复数矩阵到伪彩图调用上面这个函数画时频图。我们需要把复数的模取出来转成dB单位否则动态范围太大小能量成分在图上根本看不清。% 参数设置 winLen 256; % 窗长 win hann(winLen, periodic); % 汉宁窗 hop winLen/4; % 帧移75%重叠 nfft 1024; % FFT点数大于窗长做补零 % 计算STFT [S, f, t_center] my_stft(x, win, hop, nfft, fs); % 画图 figure; imagesc(t_center, f, 20*log10(abs(S) eps)); axis xy; % 让纵轴频率从小到大从下往上显示 xlabel(时间 / s); ylabel(频率 / Hz); colormap(jet); clim([-80 20]); % 动态范围具体值可根据情况调整 colorbar; title(STFT时频图手动实现);这里20*log10(abs(S)eps)非常关键。原始幅值谱动态范围可能从几百到0.001直接用imagesc会是一片蓝色海洋什么都看不见。转成dB后把颜色轴限制在-80dB到20dB就能同时看清强分量和弱分量。clim的取值要看具体情况比如噪声底在-100dB左右可以放宽到-100~20。另外axis xy是为了纠正imagesc默认的纵轴方向不写的话频率轴会倒过来这种低级错误排查起来特别耗神。3.4 验证与解读时频图跑完上述代码你会在时频图上看到一条从50Hz延伸到300Hz的清晰斜线这是chirp信号的频率变化轨迹在0.8秒附近、500Hz上下会出现一条竖直的亮色条纹这是瞬态冲击的宽频特征。两个特征都能对上就可以确认STFT代码写对了。我做了一个简单验证在采样率1000Hz、窗长256、重叠75%的条件下代码运行耗时为几十毫秒级别处理2秒长度的信号绰绰有余。再把同一段信号用Matlab自带的spectrogram结果对比峰值位置完全一致幅值差异在浮点误差范围内。这说明自己写的循环实现完全没有问题后续你想在硬件平台或Python里复刻逻辑也是完全相同的。4. 参数选择与踩坑记录4.1 窗长选择的实战心得窗长是STFT里最需要调的核心参数。我刚开始做时频分析时习惯所有信号都用1024点窗结果分析振动信号时把所有瞬态冲击都抹平了时间分辨率严重不足。后来总结出一条实用经验先明确你更关心“什么频率间隔需要分开”还是“什么时间点需要对准”。如果需要区分相距2Hz的两个频率成分采样率2000Hz时窗长至少要1000点如果需要定位毫秒级的冲击采样率2000Hz时窗长最多只能取几十个点。这两者不可兼得只能折中。很多老工程师的做法是“多次尝试”。先用短窗看整体形态再逐步加长窗看频谱细节是否会变清晰直到出现明显的拖影再往回退。这比死记硬背参数表更可靠。实际操作中重叠率我一般设成50%~75%。重叠率太高计算量增大但收益递减重叠率太低时频图会有明显的栅格感尤其在非平稳信号上看起来像二维码。75%是我的习惯默认值。4.2 常见问题速查表下面整理STFT实操里最容易踩的坑和对应排查方法都是我自己遇到过的。现象可能原因解决办法时频图“糊成一片”看不出明显曲线窗长太长时间分辨率不够减小窗长比如从1024降到256时频图上有很多水平横纹能量溢出窗函数旁瓣太高频谱泄漏严重换用汉宁窗或布莱克曼窗某个频率分量忽明忽暗、有空洞帧移过大时间采样太稀增加重叠率hop设为窗长的1/8~1/4图的横轴时间整体偏移了半个窗长时间坐标用窗起点而非窗中心按窗中心计算时间坐标频率轴倒置没有使用axis xy加一行axis xy颜色全是蓝色看不到细节未转dB或颜色轴范围不对用20*log10(abs(S)eps)调整clim频谱分辨率不够两个峰粘在一起窗长太短补零也无法解决增大窗长注意时间分辨率会牺牲信号首尾有异常的频率分量边缘窗覆盖不完整丢弃首尾不完整的帧或边缘补零并舍弃边缘段频谱上有直流附近的大块亮斑信号有直流偏置或加窗不平滑先去均值再检查窗形态4.3 与Matlab自带spectrogram的结果对拍为了确认手动实现没写错我建议你做一个盲测用自带spectrogram算一次与自己写的函数算一次然后对比差值。比如[Sm, fm, tm] spectrogram(x, win, winLen - hop, nfft, fs); max(max(abs(S - Sm)))如果正确这个差值应该在1e-12级别。我跑下来是1.4211e-13。这里有一个容易踩的点spectrogram默认的时间轴返回的是经过窗中心修正的我的实现也做了同样处理所以对齐得很自然。如果不处理横向时间会差四分之一个窗长左右画在同一张图上看不太出来但对比数据就会发现系统性偏移。对比的另一个价值是当你的代码结果和工具包出现差异时先排查参数比如重叠率、FFT点数、单边谱取法是不是一致。很多“看起来不对”其实都是参数没对齐造成的而不是逻辑问题。5. 经验总结与扩展方向5.1 调参之外容易被忽略的细节STFT上手不难但要把时频图用好有几个细节我强调过太多次了。一个是先看时频图的全貌再决定要不要做谱平滑或增加动态范围不要一上来就套参数。另一个是对于非平稳信号最好把STFT的帧长和物理过程挂钩。比如分析启停工况下的设备整段频率变化跨度很大这时可以用自适应窗长的STFT或多尺度分析但这就超出STFT的范围了。另外写代码时养成好的习惯函数输入参数尽量显式传入不要写死文件路径或采样率对信号先做去均值和去趋势否则直流分量会让时频图底部出现一条亮带影响观察其他弱分量。预处理这个步骤很多人一开始不在意等到调完图才回头看数据白白浪费时间。5.2 从STFT走向更复杂的时频工具理解STFT之后再去看小波变换就不难了。STFT的窗固定不变相当于在整个时频平面上用同一把尺子去量小波变换的尺度随频率变化高频处时间窗自动缩短低频处频率窗自动展宽更适合宽频带信号。这两者不是替代关系STFT仍然在语音、振动、雷达的很多场景里是首选因为它结果直接、相位信息完整、逆变换实现简单。如果要进一步做特征提取STFT矩阵还能用来算谱熵、能谱密度、边际谱或者作为二维输入送给卷积神经网络做故障分类。这些扩展方向都基于对STFT矩阵的深入理解。我建议你把自己的函数改造成支持输入信号流、边采边算的形式这在实时监测里很实用核心就是维护一个滑动缓冲区每次只对新进来的数据块做一帧STFT并不复杂。最后再分享一个个人经验在实际项目中不要迷信某一个窗长参数最好把窗长、重叠率、FFT点数写成交互式变量先用粗参数快速预览再根据信号特性细调。我做现场数据时常准备一块二次查看的工具区把“短窗看事件、长窗看边带”当作两套固定配置来回切换效果比追求一个“万能参数”好得多。做时频分析说到底是对时间与频率分辨率权衡的理解工具只是实现手段代码越底层你对这个权衡的体会就越深。