ARTICLE DETAIL

资讯详情

深耕网站视觉设计与运营推广的一线实战洞察。

STFT短时傅里叶变换与spectrogram语谱图详解:从原理到Python实战

STFT短时傅里叶变换与spectrogram语谱图详解:从原理到Python实战 简介本资源是一套面向信号处理初学者与MATLAB进阶用户的STFT短时傅里叶变换实战教学材料聚焦spectrogram函数的原理理解、交互式操作与编程实现解决非平稳信号时频分析中“看不懂参数、画不出准确图、调不好分辨率”的典型痛点。压缩包含3个核心文件1个实测EMG生理信号数据.mat、1个可直接运行的spectrogram解析脚本.m以及1个85分钟高清教学视频.mp4总大小205.12MB结构精炼、即下即用。已有5263人学习下载视频内容系统拆解为六大模块——从函数本质与适用场景到分步推演STFT计算流程从GUI界面操作要点到代码调用中窗口长度、重叠率、FFT点数等关键参数的精准设置最后涵盖时频图着色、坐标标注、动态范围优化等可视化技巧。配套脚本与实测数据支持读者立即复现、对比验证真正实现“学完即用、调参不盲”。1. STFT是什么为什么需要它在信号处理里频谱分析是绕不开的老朋友。我们手头有一段时间信号想看看它的频率成分最直接的工具就是傅里叶变换。傅里叶变换把一个完整时域信号扔到频域里得到全局的频谱这个频谱告诉我们“信号里有哪些频率”但它回答不了一个问题这些频率是什么时候出现的比如一段音乐里前两秒是低音中间两秒是高音最后两秒又是低音。你用傅里叶变换看频谱只会看到低频和高频成分都存在具体哪段旋律对应哪个时刻完全无从分辨。这就是傅里叶变换的核心局限它把时间信息完全抹平了你在频域里找不到任何和“时间顺序”有关的东西。STFTShort-Time Fourier Transform短时傅里叶变换就是冲着这个问题来的。核心思路说白了很简单把长信号切成一帧一帧的短片段假设每个片段内部信号是平稳的然后对每一帧单独做傅里叶变换再把各帧的频谱按时间顺序排开。这样一来输出结果就是一个二维矩阵——横轴是时间纵轴是频率颜色/数值代表该时刻该频率成分的能量或幅度。这种可视化结果就是我们常说的语谱图spectrogram。如果把语音信号丢进STFT语谱图上能看到明显的声纹结构——横条纹是共振峰竖条纹是不同音节的边界斜条纹是音调上扬或者下降。这是语音识别、音乐分析、故障诊断、振动检测等领域的基本功工具。我自己最初接触STFT是在做振动信号故障诊断的时候当时要在时域波形里找轴承故障特征频率波形看不出门道一跑spectrogram马上清楚了故障对应的频率成分在哪段时间出现、能量多强一目了然。对于刚开始接触信号处理的同学来说STFT是你从频域思维进阶到时频域思维的必经台阶。不会STFT你看信号只能“二选一”——要么看时域要么看频域会了STFT你可以同时看两个维度很多工程问题的答案就藏在这两个维度的交叉里。2. STFT的基本原理与核心参数2.1 分帧、加窗与傅里叶变换三步走STFT的过程拆开看就三步第一步分帧。把连续信号按固定长度切成一帧一帧帧长一般叫frame_size或n_fft单位是采样点。比如采样率 16000 Hz帧长取 400 个采样点那么一帧的时间就是 25 毫秒。这个时间长度不是随便定的语音处理里常见的帧长在 20-30 毫秒之间因为人说话时声道特征在这个时间尺度内近似稳定。第二步加窗。直接对信号片段做傅里叶变换会产生频谱泄漏问题。帧的边界是硬截断的在时域上表现为矩形窗矩形窗在频域有很高的旁瓣会让频谱出现严重的“拖尾”现象本来只有一个频率成分结果附近频率也出现较高的能量。解决办法是给每帧信号乘一个窗函数比如汉宁窗Hanning、汉明窗Hamming、布莱克曼窗Blackman等让帧的两端平滑衰减到接近零减小边界处的不连续性。第三步对每一帧做FFT。快速傅里叶变换是数字信号处理的基石算法复杂度从直接计算的O(N^2)降到了O(N log N)。对每帧做完FFT后去掉镜像部分实信号FFT结果是对称的只用前一半就行取幅度平方或对数幅度作为该帧的频谱然后按时间顺序堆叠起来就得到spectrogram矩阵。2.2 帧移、帧长与时间分辨率/频率分辨率的博弈STFT里面最考验直觉的一对矛盾是时间分辨率和频率分辨率想要时间上精细频率上就变粗想要频率上精细时间上就变粗。这个关系可以用数学式说清楚。时间分辨率取决于帧移hop_length也就是相邻两帧起点之间的距离频率分辨率取决于n_fft也就是FFT点数频率分辨率 采样率 / n_fft。以采样率 16000 Hz、n_fft400为例频率分辨率是 16000 / 400 40 Hz。如果你把n_fft提到 1024频率分辨率变成约 15.6 Hz但每帧时间覆盖从 25 ms 变成 64 ms这64毫秒内信号如果变化很快就会在语谱图上出现时间方向的“糊”感。换一个更直白的比喻手机拍照的时候要么拍得清楚但视野窄要么拍得广但每块地方都模糊。STFT的窗口就是这个镜头短窗口看到的是“瞬间的变化”长窗口看到的是“稳定的频率结构”。实际项目里怎么选我的经验是先从应用场景倒推。2.3 为什么stft的频域结果是复数很多新手第一次拿到STFT结果时会愣住为什么输出的矩阵里全是复数幅值怎么取这里要澄清一点傅里叶变换本质上是在做“信号和一系列复指数基函数的相似度计算”输出结果天然包含实部和虚部分别对应余弦分量和正弦分量的权重。实信号的FFT结果在正负频率上共轭对称所以只保留正频率部分信息不会损失。工程中使用STFT通常有两种取数的姿势取幅度谱abs(spectrogram)看各频率成分的能量变化取功率谱abs(spectrogram)**2按能量分布来分析。这两种都能画出熟悉的语谱图。需要注意的是如果后续要做信号重构比如变调不变速、时域拉伸必须保留完整的复数结果和相位信息丢掉相位直接做逆变换重构出来的时域信号基本就是噪音。2.4 频谱泄漏与窗函数选择前面说了边界截断会产生频谱泄漏这里把窗函数的选择逻辑再讲透一点。窗函数有许多种常用的有矩形窗旁瓣最高但主瓣最窄适合瞬时频率精确测量工程上用的少因为泄漏太严重汉宁窗主瓣稍宽旁瓣明显衰减是通用首选汉明窗和汉宁窗很像但两端的衰减更缓在语音识别里非常常见布莱克曼窗旁瓣更低频率分辨能力更弱适合对幅度精度要求高的场景。我个人的习惯是先上汉宁窗跑一版效果不够再试旁瓣更低的窗函数。大多数情况下汉宁窗都能给出足够清晰的语谱图没有必要把时间花在纠结窗函数上。真正影响项目成败的是n_fft和hop_length这两个参数。还需要提一个细节加窗后每帧能量会低于原始信号的帧能量因为窗函数把两端压低了。如果需要做准确的能量分析可以在计算功率谱后乘以一个能量校正系数或者使用已归一化的窗函数。大部分库函数比如librosa、scipy默认会做归一化但自己手写STFT时一定不要忘了这个坑。3. spectrogram函数用法拆解3.1 不同库的spectrogram函数对比现在说到本文的主角spectrogram函数。需要先提醒一点spectrogram这个名字在几个主流的Python库中都存在参数和返回值略有差异很多人都在这里踩过坑。我遇到过不止一个同事把scipy的spectrogram结果直接拿到matplotlib的specgram里对照发现数值对不上就开始怀疑人生。最常用的三个版本库函数返回内容scipy.signalspectrogram(x, fs, window, nperseg, noverlap, nfft)f, t, Sxxmatplotlib.pyplotspecgram(x, NFFT, Fs, noverlap)spectrum, freqs, t, imlibrosastft(y, n_fft, hop_length, win_length, window)复数矩阵Z其中scipy是最接近“纯函数”的封装适合把频谱数据拿去做后续分析matplotlib的specgram直接在绘图接口里封装了所有逻辑画图方便但要拿到中间计算结果稍微绕一点librosa是音乐和语音分析的标配它的stft返回复数矩阵需要用abs取模再转成dB值画图。除Python外MATLAB里也有spectrogram函数语法类似但默认参数不太一样比如MATLAB默认返回短时功率谱密度而scipy默认返回短时功率谱幅度。如果项目需要跨语言复现同一套算法参数的默认值一定要逐个核对我在这里吃过一次亏后面细讲。3.2 scipy.signal.spectrogram 参数逐项精讲用scipy的spectrogram函数作为主线来拆因为它的参数最全、最接近STFT的理论定义。函数签名如下scipy.signal.spectrogram(x, fs1.0, window(tukey, 0.25), npersegNone, noverlapNone, nfftNone, detrendconstant, return_onesidedTrue, scalingdensity, axis-1, modepsd)逐个说关键的x输入信号一维数组。超过一维要指定axis参数默认是最后一维。fs采样率默认 1.0。这个参数只影响返回的时间数组t和频率数组f的单位不影响Sxx的数值本身。所以哪怕你忘了填fsSxx照样是对的只是横纵轴的坐标刻度的物理意义没了。window窗函数可以是字符串、元组或数组。字符串如hann、hamming、blackmanharris元组如(tukey, 0.25)表示Tukey窗第二项是锥形部分的占比。如果传入一个数组则该数组被直接用作窗函数数组长度需要和nperseg一致。这里有个容易迷惑的点window默认是(tukey, 0.25)而不是很多人以为的汉宁窗。Tukey窗中间有一段是平的类似于矩形窗和汉宁窗的折中。如果习惯性地认为默认是汉宁窗结果频谱旁瓣偏高排查半天发现是窗的问题。实际项目里我一般显式指定windowhann不依赖默认值。nperseg每帧的长度单位是采样点。如果没指定默认是信号长度的八分之一且取256的上限也就是如果信号长度大于2048默认就是256。这个默认值在很多应用里偏小了尤其是低频分析。建议手动指定。noverlap相邻帧重叠的采样点数量。如果不指定默认是nperseg // 2即50%重叠率。重叠率决定语谱图在时间方向的平滑程度也决定计算量。50%是比较常规的起点高分辨率需求可以用75%甚至更高。注意noverlap必须小于nperseg否则报错。nfftFFT点数。如果nfft大于nperseg函数会对帧做零填充相当于在频率方向上插值让频谱看起来“更光滑”。零填充不会增加真实频率分辨率纯粹是视觉上的平滑效果不要被它迷惑。return_onesided对于实数信号返回单边频谱True频率范围从0到奈奎斯特频率对于复数信号必须设为False返回双边频谱频率范围从负奈奎斯特到正奈奎斯特。scaling这个是很多人忽略的关键参数。density表示返回的是功率谱密度PSD单位是 V^2/Hzspectrum返回的是幅度平方谱单位是 V^2。做随机振动分析或噪声分析时用density做确定性信号分析时用spectrum。数值上两者差一个频率分辨率的倍数不看文档直接取Sxx做比较很容易得出错误结论。mode返回值形式。psd返回功率谱密度magnitude返回幅度谱angle返回相位谱complex返回复数谱。按需选择默认psd。3.3 librosa.stft 与 scipy 的差异语音和音乐方向的朋友大概率会用到librosa。librosa的stft和scipy的spectrogram有几点明显不同librosa.stft返回的是复数矩阵形状为(1 n_fft/2, 1 len(y)/hop_length)第一维是频率第二维是时间帧。librosa.stft没有noverlap参数而是用hop_length直接控制帧移配合win_length控制窗长。librosa默认窗是汉宁窗windowhann和scipy默认不同这一点对初学者来说反而更友好。librosa默认对信号做居中填充centerTrue会在信号前后各补n_fft // 2个零保证每一帧的中心对齐到时间网格上。两个库的返回数值也存在差异。scipy的spectrogram在scalingspectrum、modemagnitude时返回的是幅度谱而librosa的stft返回的复数矩阵直接取模得到的是幅度谱。两者在数值上会差一个跟窗函数有关的增益系数严格对比时需要用窗函数均值做归一化。如果只是做可视化这个差异无所谓如果要做特征提取并跨库验证就得小心了。3.4 从spectrogram矩阵到可视化语谱图拿到Sxx矩阵之后画图也有讲究。很多人直接imshow(Sxx)完事结果图像上一片漆黑或一片惨白这是因为频谱的动态范围太大了——能量高的频率分量和能量低的分量可以相差好几个数量级线性色标根本显示不出细节。常规做法是转成对数刻度或分贝dB刻度。分贝的计算公式是dB 10 * log10(Sxx / ref)其中ref是参考值通常取Sxx的最大值或取1表示相对于单位幅度的分贝值。librosa.display.specshow封装了这个过程内部会转成 dB 显示。如果自己用 matplotlib 画可以用pcolormesh(f, t, 10 * np.log10(Sxx 1e-10), shadingauto)来避免大动态范围下的细节丢失加一个极小的常数是为了防止log10(0)产生 -inf。4. 实操用Python完成一次完整时频分析4.1 准备测试信号与可视化环境纸上谈兵到这里直接上一个完整实例。假设我们要分析一段模拟信号2秒时长的采样率8000 Hz信号由两个频率分量组成前1秒是400 Hz正弦波后1秒是1500 Hz正弦波并叠加一点高斯白噪声。先构建环境import numpy as np import matplotlib.pyplot as plt from scipy.signal import spectrogram fs 8000 duration 2.0 t_total np.linspace(0, duration, int(fs * duration), endpointFalse) # 前1秒 400 Hz后1秒 1500 Hz加10%噪声 seg1 0.8 * np.sin(2 * np.pi * 400 * t_total[:fs]) seg2 0.8 * np.sin(2 * np.pi * 1500 * t_total[fs:]) noise 0.1 * np.random.randn(len(t_total)) x np.concatenate([seg1, seg2]) noise4.2 设置合理的STFT参数按照上文思路来确定参数fs 8000信号主要能量集中在1500 Hz所以频率分辨率要足够区分400 Hz和1500 Hz取nperseg 256频率分辨率 8000 / 256 31.25 Hz完全够用noverlap nperseg // 2 12850%重叠时间方向上帧间隔为128 / 8000 16 ms足够观察前1秒到后1秒的频率切换window hann通用首选scaling spectrum因为这里是确定性信号看幅度平方谱合适mode psd或magnitude都可以这里选psd来绘制功率谱。核心代码f, t, Sxx spectrogram( x, fsfs, windowhann, nperseg256, noverlap128, scalingspectrum, modepsd ) print(频率轴形状:, f.shape) print(时间轴形状:, t.shape) print(语谱矩阵形状:, Sxx.shape)4.3 绘制语谱图并解读结果画图时注意转成dB刻度plt.figure(figsize(10, 6)) plt.pcolormesh(t, f, 10 * np.log10(Sxx 1e-12), shadingauto, cmapmagma) plt.xlabel(时间 (s)) plt.ylabel(频率 (Hz)) plt.title(STFT Spectrogram) plt.colorbar(label功率 (dB)) plt.ylim(0, 2500) plt.tight_layout() plt.show()运行结果应该能直观看到两段时间的频率差异前1秒在400 Hz附近有一道亮纹后1秒在1500 Hz附近出现另一道亮纹。因为加了噪声背景有均匀的细碎噪点这是正常现象。如果nperseg放大到1024你会看到频率方向上更细腻的纹理但400 Hz→1500 Hz的切换时刻会变得模糊切换点附近会出现一段过渡的斜坡——这就是时间分辨率下降的表现。有一个容易被忽略的细节matplotlib的pcolormesh默认以像素中心对齐坐标而t和f数组表示的是每个bin的边界还是中心不同库的约定不同。scipy.signal.spectrogram返回的t和f是每个bin的中心坐标可以直接用作pcolormesh的坐标。如果发现图像整体偏移了半个格子多半是坐标约定没对齐。4.4 用librosa复现同样的分析为了对照再给出librosa版本语音或者音乐方向的同学可以直接抄import librosa import librosa.display n_fft 256 hop_length 128 Z librosa.stft(x, n_fftn_fft, hop_lengthhop_length, windowhann) S_db librosa.amplitude_to_db(np.abs(Z), refnp.max) plt.figure(figsize(10, 6)) librosa.display.specshow(S_db, srfs, hop_lengthhop_length, x_axistime, y_axishz, cmapmagma) plt.colorbar(format%2.0f dB) plt.title(librosa STFT spectrogram) plt.ylim(0, 2500) plt.tight_layout() plt.show()librosa.display.specshow会自动处理坐标轴、刻度标签以及频率轴从Hz到kHz的转换非常省事。librosa.stft返回的帧数比scipy版本多因为它默认做了居中填充。如果要在两个库之间做特征对齐比对建议在调用librosa.stft时显式指定centerFalse并手动处理信号边界。5. 实战中的坑与经验STFT参数选择与问题排查5.1 参数选择的经验法则我在实际项目里总结了一套选参流程可以当成checklist来用先定频率分辨率需求。比如要区分相距20 Hz的两个频率分量采样率16000 Hz那么nperseg至少要 16000/20 800向上取整到2的幂就是1024。再定时间分辨率需求。如果要用语谱图观察事件开始和结束的时刻帧移决定了时间精度。比如需要10 ms的时间精度hop_length fs * 0.01 160那么noverlap nperseg - hop_length 864重叠率约84%计算量会变大但效果能保证。最后做一次可视化验证。参数选完之后用包含已知频率切换的测试信号跑一遍看切换沿是否清晰、频率是否准确再决定是否需要微调。这个方法比盯着文档啃有效得多。5.2 常见报错与排查记录问题一ValueError: noverlap must be less than nperseg这个报错的原因很直接noverlap大于等于nperseg了。如果想让帧移为0即完全不重叠noverlap应该设为0而不是nperseg。检查代码里noverlap和nperseg的赋值逻辑尤其注意当nperseg被动态计算时noverlap是否同步更新。问题二频率轴数据与预期不符峰值频率偏了常见原因有三类第一采样率fs传错了第二输入的信号是复数但没有设置return_onesidedFalse导致频谱折叠第三信号里有直流分量0 Hz或高频噪声动态范围太大看起来峰值被“压”到了边缘。排查顺序建议先从检查输入信号长度和fs开始再用一个已知频率的纯正弦波做标定一般能快速定位。问题三语谱图全是黑的或花的语谱图颜色分布异常大多不是STFT本身错了而是显示端的问题。线性色标遇到大动态范围信号就会这样。解决办法是转成dB刻度并合理设置色标的上下限。另外Sxx里存在全零帧比如静音段时log10(0)会产生警告加一个小epsilon值就能避免例如10 * np.log10(Sxx 1e-12)。问题四librosa和scipy结果对不上前面提到了两者在数值定义上有差异。librosa返回幅度谱且没有做PSD归一化scipy在scalingdensity时返回归一化的功率谱密度。如果要把两者对齐可以用scipy.signal.spectrogram(x, fs, windowhann, nperseg256, noverlap128, scalingspectrum, modemagnitude)拿到幅度谱再除以窗函数的均值和abs(librosa.stft(x, n_fft256, hop_length128, centerFalse))对比。两者形状和数值基本一致差异只在于归一化方式和边界填充策略。5.3 逆变换与信号重构的注意事项STFT更进阶的玩法是逆变换istft。当需要做时域拉伸、变调、去噪时通常流程是分帧加窗→STFT→在时频域做处理如掩码去噪→ISTFT还原成时域信号。这里有两个非常容易被忽略的点窗函数与重构的匹配不是所有窗函数都适合做完美重构。Hann窗配合50%重叠率满足COLAConstant Overlap-Add条件可以做到精确重构。如果你换了一种窗或者改了重叠率重构出来的信号可能带幅度调制噪声。相位信息必须保留很多人做去噪时只处理幅度谱保留原始相位然后做ISTFT。这不是最优的但至少能听。如果相位也被修改了重构信号会出现明显的金属感或破碎感。要真正提高质量就得用相位恢复算法这个超出本文范围但值得在实战前大概了解。5.4 性能优化大数据量下的STFT怎么跑得快STFT的计算量主要受nfft、帧数和信号长度影响。我们做一个简单估算假设采样率44100 Hz60秒音频信号长度约264.6万点n_fft2048、hop_length512帧数约5168帧。每一帧做一次2048点FFT总FFT次数5168次单次2048点FFT大约几十微秒量级总耗时在0.2秒左右。普通电脑都能轻松处理。但如果hop_length缩到128帧数变成约20672帧计算量上升4倍对于长时间实时流式处理就要小心了。优化手段从上到下依次排把n_fft固定为2的幂利用FFT的基-2算法加速在可接受范围内减少重叠率将信号分块后用多线程并行处理或者用GPU加速库如cuSignal对实时场景用流式STFT每次只处理一个块而不是等完整信号到位后再处理。5.5 几个容易记混的细节最后整理几个我踩过坑之后才记住的细节。spectrogram返回的t数组第一个值不一定从0开始这和边界填充策略有关。scipy默认不填充第一帧从位置0开始librosa默认填充第一帧的窗口中心在n_fft // 2位置所以t[0]是一个很小但不等于0的值。画图的时候横轴要对上参考时刻这个偏移要心里有数。nperseg和nfft是两个不同的参数前者是实际参与加窗的信号长度后者是FFT运算点数。nfft nperseg时才发生零填充相当于频域插值。有些库把这两个参数合并了比如MATLAB的nfft同时控制两者迁移代码时尤其容易搞混。STFT的帧数计算公式是1 (len(x) - nperseg) // hop_lengthscipy不填充的情况或1 len(x) // hop_lengthlibrosa填充的情况。如果和你拿到的矩阵形状对不上多检查一步边界填充的默认设置。6. 延伸从spectrogram到更高级的时频分析工具STFT是时频分析的地基但它在实际使用中有一个绕不开的天生短板窗口长度固定分辨率的博弈无法调和。短窗口看时间细节长窗口看频率细节一旦选定了nperseg整个分析过程中的时频分辨率就锁死了。对于包含瞬态冲击和长时稳定成分的混合信号比如机械设备故障信号、生物电信号固定窗口的STFT往往顾此失彼。对症的武器是小波变换Wavelet Transform。小波变换用不同尺度的小波基函数去匹配信号高频段用窄窗口、低频段用宽窗口天然具备“变焦”能力。在Python里用pywt.cwt可以快速得到小波尺度图它和spectrogram在表现形式上很像但分辨率特性完全不同。如果信号里既有明显的瞬态冲击比如轴承故障的脉冲又有稳定的周期性分量小波尺度图往往比STFT语谱图更容易看清。还有一种常见变体是梅尔频谱Mel-spectrogram它在STFT结果基础上把频率轴按照梅尔刻度重映射模拟人耳对频率的非线性感知。语音识别、音乐流派分类等任务里梅尔频谱就是标配输入特征。librosa.feature.melspectrogram一行代码就能完成从波形到梅尔频谱的转换但它的底层就是先算STFT再做梅尔滤波器组积分理解STFT是理解梅尔频谱的必要前提。对这些扩展工具我的建议是先把STFT吃透能把nperseg、hop_length、window、scaling这些参数对结果的影响在心中模拟出一个大概的图像再上手小波和梅尔频谱你会发现它们的原理就容易很多了。很多初学者直接跳到高级工具结果参数调不动就是因为在STFT这一段缺了“手感”。最后再分享一个小技巧在调试STFT参数时我习惯构造一个“频率跳变纯音脉冲噪声”的合成信号作为测试样例每次修改nperseg或hop_length后跑一遍语谱图用眼睛观察频率跳变沿的陡峭程度和频率峰的宽窄。这比盯着参数文档空想要直观得多也能帮你快速建立对时频分辨率博弈的直觉。这个习惯我一直保留到现在遇到不确定的参数组合就先跑一版合成信号验证省了不少排查时间。本文还有配套的精品资源点击获取
返回列表