尧图网站设计 尧图网站设计YAOTU DESIGN
ARTICLE DETAIL

资讯详情

深耕网站设计与一线实操的经验洞察。

MATLAB时频分析工具箱实战:从短时傅里叶到Wigner-Ville

MATLAB时频分析工具箱实战:从短时傅里叶到Wigner-Ville 简介面向MATLAB信号处理学习与研究人员的时频分析工具箱压缩包涵盖短时傅里叶变换、Wigner-Ville分布等常用时频分析方法的函数实现可用于信号特征提取、时频图绘制与算法对比验证。包体共146个文件以139个m脚本为主要代码主体辅以少量mat数据文件、txt说明文档及htm使用帮助压缩包约961KB结构清晰紧凑。内容包含从TFDEMO基础演示到TFRVIEW、TFRQVIEW等交互式时频显示工具配套的CONTENTS文件与新手必看帮助文档能够引导用户快速上手理解各函数的输入输出与适用场景。目前已有743人学习下载借助这些脚本和示例研究者可直接修改参数、替换测试信号将时频分析能力迁移至语音、振动或生物医学信号处理项目中是算法学习与工程验证的实用参考。 做信号处理的人早晚会碰到一个问题傅里叶变换能告诉你信号里有哪些频率却说不清这些频率在什么时候出现。这时候就得请出时频分析。我最早接触MATLAB时频分析工具箱是因为一组轴承振动数据怎么都看不出故障特征后来发现不是数据问题是分析工具没选对。这套基于经典时间-频率分析理论的工具箱把我从重复造轮子的状态里解放出来也是我当时入门时频分析最顺手的工具。无论你是在校学生还是做设备状态监测、语音识别、雷达信号处理的工程师只要跟非平稳信号打交道这个工具箱都值得花半小时装上试试。1. 这个工具箱到底是什么解决什么问题1.1 傅里叶变换的“盲区”和时频分析的登场傅里叶变换是信号处理的地基但它有一个绕不过去的默认前提信号是平稳的或者至少在统计特性上不随时间变化。可现实中的信号大多数不听话。轴承磨损产生的冲击、语音中的音节、雷达回波里的多普勒变化全是频率随时间变化的结果。拿傅里叶变换去做这类信号只能得到一个平均光谱时间维被彻底抹掉了。比如一段持续1秒、频率从100Hz线性扫到600Hz的调频信号频谱上只有一个横跨100到600Hz的宽包络根本看不出扫频方向更看不出每个时刻对应的频率。时频分析就是为了解决这个矛盾出现的它把一维时间信号映射到二维的时间-频率平面上让频率随时间的变化轨迹直接“画”出来。1.2 工具箱里都有什么它能用在哪些地方MATLAB时频分析工具箱很多资料里叫TFTBTime-Frequency Toolbox不是MATLAB自带的那个信号处理工具箱而是一个面向科研和教学的开源函数集合。它的代码最早由法国一批信号处理实验室发布核心覆盖了三大类算法第一类是线性时频表示比如短时傅里叶变换STFT、连续小波变换第二类是双线性时频分布也就是Cohen类中的Wigner-Ville分布、Choi-Williams分布、伪平滑Wigner分布等第三类是时频后处理工具比如脊线提取、瞬时频率估计、时频重排。我上学时拿到这个压缩包第一反应是函数名字怎么都这么怪比如tfrstft、tfrwv、tfrcw后来翻源码才明白这些封装得很好的函数调用前甚至不用理解全部公式就可以先看到效果。这套工具箱特别适合三类人一是刚接触时频分析、想快速看到算法效果的学生二是做故障诊断或生物医学信号处理的工程师需要先验证某个时频分布是否适合数据三是想研究算法内部实现、读源码做论文复现的科研人员。因为它提供的是MATLAB源码在编辑器里按CtrlD就能看到每一步计算比纯看公式直观得多。用上之后很多需要从零写的算法直接变成了三五行代码。2. 拿到压缩包之后安装、路径与版本兼容2.1 解压与路径设置压缩包解压后通常会有一个主目录和若干子目录比如tftb、tftb/refbooks之类。有人直接双击进去点tfrstft.m运行结果MATLAB报“未定义函数或变量”十有八九是路径没设。我的建议是先把整个目录解压到一个没有中文和空格的路径下比如D:\Tools\TFTB然后在MATLAB命令行执行一句addpath(genpath(D:\Tools\TFTB)); savepath;addpath负责添加路径genpath会把该目录下所有子文件夹一并加进来savepath把路径保存到永久配置里下次启动不用重复设置。这里有个细节如果只addpath(D:\Tools\TFTB)子文件夹里的函数照样找不到所以一定要用genpath。用图形界面的朋友可以在“设置路径”对话框里点“添加并包含子文件夹”效果一样。需要说明的是这个压缩包并不是官方MathWorks工具箱因此不会出现在“附加功能”里必须手动加路径。2.2 验证安装是否成功安装完别急着跑大数据先用最简单的命令验证。在命令行输入which tfrstft如果返回的是D:\Tools\TFTB\tfrstft.m这样的路径说明路径已经生效。再敲help tfrstft能看到详细的帮助文本有些版本帮助文本是法语或半英文不影响理解核心参数。然后我习惯跑一个demo验证整体环境。比如生成一段简短信号t (0:127) / 100; x sin(2*pi*30*t); [tfr, tIdx, fIdx] tfrstft(x, 1:length(x), 64, hamming(31)); imagesc(tIdx/100, fIdx*100, abs(tfr)); axis xy;如果图形窗口里出现明显的水平亮带说明工具箱工作正常。这一步很快能避免后续分析时因为环境问题产生不必要的怀疑。2.3 不同MATLAB版本的兼容性注意这个工具箱诞生得早不少函数还是上世纪九十年代的风格在MATLAB R2016之后使用基本没问题但偶尔会遇到个别绘图函数因为MATLAB绘图体系更新而报警告。印象最深的是tfrdraw这个内部绘图函数旧版能自动生成带坐标轴的图形新版可能提示某些图形属性被移除。遇到这类情况我的办法是不依赖它内部绘图自己用imagesc或contourf手动画结果更可控。另一个兼容性问题是数据类型老代码默认输入是double的列向量如果喂进去一个逻辑数组或单精度数组某些版本的矩阵运算会报维度错误。所以建议在处理前统一用double(x(:))转换。还有一点这类开源工具箱不会随MATLAB版本自动更新如果你电脑上装了R2021a又装了R2023a记得在每个版本里分别执行addpath或者直接写一个startup.m启动脚本里面放好addpath(genpath(...))省得每次换版本都折腾。3. 核心算法与使用要点不只会上手还要知道为什么3.1 短时傅里叶变换参数选择的门道短时傅里叶变换是最容易上手的时频分析工具核心思想就是把信号切成小段每段加一个窗函数后再做FFT最后把所有段的谱按时间拼起来。工具箱里对应tfrstft函数。但真正用好它关键是懂窗长。窗长越长频率分辨率越高但时间分辨率越差窗长越短越能捕捉瞬态变化但频率会糊成一片。这就像用相机拍一个正在移动的人快门慢了背景清楚但人模糊快门快了人停住了但噪声变多。实际选窗长没有一个放之四海皆准的值我调试轴承信号时常从127点起步如果希望看清两个相距很近的频率峰就加长到255或511如果希望定位冲击发生的时刻就缩短到31或63。重叠率也很重要相邻窗之间设50%到75%的重叠能让时频图更平滑避免出现明显的栅栏感。代码里对应的参数就是窗函数h其他参数如N是FFT点数可以取为2的幂次注意N只影响频域插值密度不能改变真实分辨率。3.2 Wigner-Ville分布与交叉项问题如果只用STFT很多人会觉得时频图上频率轨迹“胖”不够精细。于是有人转向Wigner-Ville分布工具箱里对应tfrwv。它的理论时频聚集性非常好单分量线性调频信号会变成一条很细的斜线看起来惊艳。但Wigner-Ville是双线性变换对多分量信号会产生交叉项这是它的“阿喀琉斯之踵”。什么叫交叉项简单说如果信号是频率为100Hz和300Hz的两个正弦叠加WVD时频图上除了100Hz和300Hz两条真实线在200Hz附近会出现一条虚假的强度振荡亮带强度甚至可能接近真实分量容易误导人。我一开始就在这上面栽过跟头拿一段含两个频率成分的仿真信号跑出漂亮的图却把交叉项当成了真实故障频率差点闹出笑话。解决交叉项的办法是使用平滑类分布比如Choi-Williams分布tfrcw或平滑伪Wigner分布tfrspwv它们在时频聚集性和交叉项抑制之间做了折中。工具箱里这些函数都封装好了代价是计算量变大、参数比如核函数宽度需要额外调节。如果你只是做快速探索可以先跑STFT确定信号大概有几个分量再决定要不要上双线性分布。3.3 瞬时频率估计与特征提取时频分析不只是画图好看更重要的用途是估计瞬时频率。工具箱里提供了instfreq等函数可以直接从信号本身提取瞬时频率也可以先计算时频分布再用脊线检测策略找到能量最大值对应的频率轨迹。在旋转机械故障诊断中轴承某个部件出现局部缺陷时振动信号的瞬时频率会出现周期性波动通过时频图上的脊线能把这种波动拉出来。具体做法是对时频矩阵沿频率方向找最大值得到每个时刻的瞬时频率再用中值滤波去掉离群点。这个过程在本文第4部分我会给出示例代码。需要提醒的是瞬时频率估计在信号瞬时带宽较窄时才可靠如果某个时刻同时存在多个强分量最大值法会在这几个分量之间跳变视觉上表现为脊线断裂或穿插这时要考虑多分量分解或先滤波分离后再逐一估计。4. 实操从信号构造到时频图制作4.1 构造一个非平稳测试信号为了演示我自己搭一个仿真信号模拟现实中的“频率斜坡短时脉冲”组合。采样率设2000Hz时长1秒前0.8秒频率从100Hz线性扫到600Hz0.5秒附近叠加一段0.1秒的300Hz正弦脉冲。这样的信号既包含频率随时间连续变化又包含短时瞬态很适合用来对比不同时频分布的表现。fs 2000; t (0:1999) / fs; x chirp(t, 100, 0.8, 600, linear); pulseIdx 1001:1200; % 0.5s到0.6s x(pulseIdx) x(pulseIdx) 1.2 * sin(2*pi*300*t(pulseIdx)); x x(:);生成后可以用plot先看一眼但plot看到的是混叠在一起的时域波形很难确认频率变化规律这正是需要时频图的地方。4.2 调用工具箱生成时频图先调用tfrstft生成STFT时频图再调用tfrwv生成WVD时频图并排对比。figure(Color, w); subplot(2,1,1); [tfr, tIdx, fIdx] tfrstft(x, 1:length(x), 512, hamming(255)); imagesc(tIdx/fs, fIdx*fs, abs(tfr)); axis xy; xlabel(Time (s)); ylabel(Frequency (Hz)); title(STFT); colorbar; subplot(2,1,2); [tfrWv, ~, fWv] tfrwv(x, 1:length(x), 512); imagesc(tIdx/fs, fWv*fs, abs(tfrWv)); axis xy; xlabel(Time (s)); ylabel(Frequency (Hz)); title(WVD); colorbar;运行后STFT图上能看到一条从100Hz延伸到600Hz的斜亮带中间脉冲处300Hz附近多出一小段水平亮线。WVD图上同样有这两部分但斜亮带内外会出现一些干涉条纹尤其曲线和水平线相交的地方这正是交叉项。第一次看的时候别被这些条纹吓到这是双线性时频分布的正常现象。这里要特别说明坐标换算工具箱里返回的fIdx一般是归一化频率范围在0到0.5之间实际频率需要乘以采样率fs。如果画出来的图频率范围不对多半是忘了乘fs或者使用了length(x)而不是FFT点数。另外tfrwv计算量比tfrstft大信号较长时可能跑得很慢我的习惯是先截取一段数据测试参数再跑全段。4.3 从时频分布中提取瞬时频率有时候光看图还不够需要把时频脊线以数值形式取出来。我常用的一种简便方法是在指定频率区间内对STFT时频矩阵沿频率轴找极大值得到每个时间点的瞬时频率。以刚才的信号为例限制在80Hz到650Hz之间直接观察整体脊线freqRange fIdx * fs; maskRange (freqRange 80) (freqRange 650); [~, imax] max(abs(tfr(maskRange,:)), [], 1); subFreqs freqRange(maskRange); fInst subFreqs(imax); figure; imagesc(tIdx/fs, subFreqs, abs(tfr(maskRange,:))); axis xy; hold on; plot(tIdx/fs, fInst, r., MarkerSize, 4); xlabel(Time (s)); ylabel(Frequency (Hz)); title(STFT with ridge);运行后红色点会大致贴合STFT图上那条从低到高的斜亮带。需要说明的是这种最大值法实现最简单但遇到多分量时容易跳变。更稳健的做法是用动态规划或惩罚路径算法搜索连续脊线工具箱里有些后处理函数也有第三方代码可参考但基础的max法已经能覆盖很多单分量场景。5. 常见问题与排查技巧实录5.1 三个最常见的报错报错信息常见原因解决办法未定义函数或变量 tfrstft工具箱路径未添加addpath(genpath(...))重启前savepath输入矩阵维度不一致x不是列向量或t长度不匹配统一用x(:)t (0:length(x)-1)图形窗口空白或报字体错误老绘图函数与新版MATLAB不兼容不用tfrdraw改用imagesc/contourf遇到过太多人卡在第一个报错上其实工具箱本身没坏就是MATLAB找不到它。路径设置好之后关掉MATLAB再重开一次确认savepath生效基本就能解决。维度问题则多半是把行向量直接传给了函数TFTB老版本对输入格式很敏感加一个(:)转置能省掉一半麻烦。图形问题属于版本历史债建议绕开内部绘图函数自己控制坐标轴反而更灵活。5.2 工具箱与MATLAB自带函数的选型MATLAB自带的Signal Processing Toolbox里有spectrogram、pspectrumDSP System Toolbox里有wvd也能做时频图。很多读者会问我那还有必要用这个老工具箱吗我的看法是分场景。自带函数接口规范、计算效率高、文档全适合直接在生产代码里用而这个第三方工具箱的优势是算法齐全、源码开放适合学习原理和快速验证论文方法。比如我想对比STFT、WVD、Choi-Williams三种分布的效果在TFTB里只需要换一个函数名在自带函数里就要倒腾各种参数对象。反过来如果只是画个语谱图我肯定用spectrogram因为它自动处理了归一化和坐标少踩很多坑。场景推荐快速查看语谱图、做常规时频显示MATLAB自带spectrogram/pspectrum研究多种时频分布、复现论文算法TFTBtfrstft/tfrwv/tfrcw等需要集成到Simulink或大规模工程MATLAB自带函数优先第三方工具箱只做预研5.3 分析结果的可视化调整时频图默认的colormap在不同显示器下可能差异很大我建议统一用parula或jet并在图旁加colorbar标明强度归一化单位。导出图片时用exportgraphics(gcf, filename.png, Resolution, 300)可以保证清晰不要用老式print -dpng容易四周留白。还有个小技巧当信号动态范围很大时abs(tfr)直接画图会把弱分量淹没这时候可以对幅值取对数比如imagesc(tIdx/fs, subFreqs, 20*log10(abs(tfr(maskRange,:))eps))弱脊线立刻清晰可见。不过要注意log会同时放大噪声参数调节要适度。我在实际项目里用这个工具箱最大的感受是它像一个时频算法的“试验台”算法一多反而逼着你去理解每类分布的本质区别。踩过交叉项的坑之后我现在拿到信号会先问自己要抓的是稳态频率还是瞬态冲击信号是单分量还是多分量对这些问题的回答直接决定了选STFT、WVD还是平滑Wigner。最后再分享一个很实用的小习惯跑任何时频分析前先对信号做一次零均值化和去趋势否则低频大能量会把整个时频图压扁很多细节根本看不见。祝你们也能用这个工具箱挖出信号里隐藏的规律。本文还有配套的精品资源点击获取
返回列表