从频谱分析到小波变换：MATLAB实战指南（附完整代码实现）

1. 从时间域到频率域信号分析的起点第一次接触信号处理时我最困惑的就是为什么要做频谱分析。直到有次用麦克风录下一段钢琴曲看着示波器上跳动的波形却完全听不出旋律才明白时间域波形的局限性。傅里叶变换就像给声音装上了频率眼镜让我们能直接看到构成复杂信号的各个频率成分。在MATLAB中实现傅里叶变换简单得惊人。假设我们有个包含50Hz和120Hz正弦波的混合信号Fs 1000; % 采样频率 t 0:1/Fs:1-1/Fs; % 时间向量 x 0.7*sin(2*pi*50*t) sin(2*pi*120*t); % 混合信号 Y fft(x); % 快速傅里叶变换 P2 abs(Y/length(x)); % 双侧频谱 P1 P2(1:length(x)/21); % 单侧频谱 f Fs*(0:(length(x)/2))/length(x); % 频率轴 plot(f,P1) % 绘制频谱图运行这段代码你会清晰地看到50Hz和120Hz两个峰。但当我用这段代码分析鸟鸣录音时发现了一个严重问题频谱只能告诉我存在哪些频率却无法知道这些频率何时出现。就像知道菜里有盐和糖但不知道何时加盐、何时加糖——这就是传统傅里叶变换对非平稳信号的无力感。2. 短时傅里叶变换给频谱加上时间维度为了解决这个问题我尝试了短时傅里叶变换(STFT)。原理很简单把长信号切成小段每段分别做傅里叶变换。这就像用手机拍视频——每秒24帧的静态照片连续播放就形成了动态画面。MATLAB的spectrogram函数让STFT实现变得直观% 生成频率变化的扫频信号 x chirp(t,20,t(end),100); window hann(128); % 汉宁窗 noverlap 120; % 重叠样本数 nfft 256; % FFT点数 spectrogram(x,window,noverlap,nfft,Fs,yaxis)这里有几个关键参数需要把握窗口长度我常用128-256个样本。太短会导致频率分辨率差就像近视眼太长会模糊快速变化像快门太慢拍运动物体重叠样本通常取窗口长度的75%-90%我习惯用120窗函数汉宁窗(hann)适合大多数情况矩形窗适合瞬态信号实际调试时我发现窗口选择对结果影响巨大。有次分析发动机振动信号用默认参数完全看不出异常把窗口从256减到128后瞬间在频谱图上看到了明显的冲击特征——这正是轴承故障的典型表现。3. 小波变换多分辨率分析的利器STFT虽然实用但固定窗口尺寸始终是个硬伤。记得分析地震信号时低频部分需要大窗口看细节高频部分又需要小窗口定位精确时间STFT无论如何调整参数都难以兼顾。这时**连续小波变换(CWT)**就像量身定制的解决方案。MATLAB的cwt函数内置了多种小波load mtlb; % 载入语音信号 [cfs,frq] cwt(mtlb,Fs,amor); % 使用Morlet小波 contour(t,frq,abs(cfs)) % 绘制时频等高线图 colorbar小波变换最神奇的特性是自适应分辨率低频区用拉伸的小波高频率分辨率高频区用压缩的小波高时间分辨率。这就像用可调焦显微镜观察信号——需要看整体结构时用低倍镜查看细节时切换到高倍镜。我在EEG信号分析中发现alpha波(8-13Hz)和gamma波(30-80Hz)需要不同的观察尺度。用STFT要么看不清alpha的精确频率要么抓不住gamma的突发时间而CWT可以同时清晰展示两者特征。4. 实战对比三种方法的综合应用去年参与一个工业设备监测项目时完整经历了从频谱分析到小波变换的技术升级。初期用FFT发现某轴承存在12kHz共振但无法确定是固有特性还是故障改用STFT后发现该频率会周期性出现暗示可能存在剥落缺陷最后通过CWT精确锁定了冲击发生时刻定位到具体旋转周期对应的机械位置。这是完整的MATLAB对比代码% 生成模拟故障信号 t 0:1/1e4:1; x sin(2*pi*50*t); % 基频振动 x(5000:500:end) x(5000:500:end) ... % 周期性冲击 0.5*exp(-1000*(t(5000:500:end)-t(5000)).^2).*sin(2*pi*12e3*t(5000:500:end)); % FFT分析 Y fft(x); P2 abs(Y/length(x)); P1 P2(1:length(x)/21); f 1e4*(0:(length(x)/2))/length(x); subplot(3,1,1) plot(f,P1) title(FFT频谱) % STFT分析 subplot(3,1,2) spectrogram(x,256,250,256,1e4,yaxis) title(STFT时频图) % CWT分析 subplot(3,1,3) cwt(x,1e4,amor) title(小波变换)从输出结果可以明显看出FFT只能看到12kHz成分存在STFT能发现周期性而CWT还能清晰显示每次冲击的精确时间位置和频率衰减过程。5. 参数调优与可视化技巧经过多个项目实践我总结出一些实用技巧对于STFT语音信号常用256点汉明窗重叠75%机械振动建议用1024点以上矩形窗重叠50%使用surf函数创建三维时频图更直观[s,f,t] spectrogram(x,window,noverlap,f,Fs); surf(t,f,10*log10(abs(s)),EdgeColor,none) view(0,90)对于CWTMorse小波(default)适合大多数场景冲击信号推荐使用Bump小波调整尺度向量可以聚焦关键频段freqrange [100 2000]; % 重点观察100-2000Hz scales centfrq(morse)*Fs./freqrange; cwt(x,scales,morse,1/Fs)在论文绘图时我习惯添加以下美化代码colormap(jet) % 使用jet色图 set(gca,YScale,log) % 对数频率轴 xlabel(Time (s)) ylabel(Frequency (Hz)) title(时频分析结果) colorbar(southoutside) % 横向色标6. 完整工具箱的实现为了方便日常使用我将常用功能封装成三个函数1. 智能STFT分析函数function [s,f,t] smartSpectrogram(x,Fs,varargin) % 自动选择最佳窗口参数 if isempty(varargin) N min(256,length(x)/10); window hann(N); noverlap floor(0.75*N); nfft max(256,2^nextpow2(N)); else % 参数输入处理... end [s,f,t] spectrogram(x,window,noverlap,nfft,Fs); end2. 交互式CWT查看器function cwtViewer(x,Fs) figure(Name,CWT分析工具) uicontrol(Style,popup,String,{Morse,Morlet,Bump},... Callback,updatePlot); % 更多交互控件... function updatePlot(src,~) waveletType src.String{src.Value}; cwt(x,Fs,waveletType) end end3. 时频对比分析工具function tfCompare(x,Fs) subplot(1,3,1) plot((0:length(x)-1)/Fs,x) title(时域信号) subplot(1,3,2) smartSpectrogram(x,Fs) title(STFT分析) subplot(1,3,3) cwt(x,Fs,amor) title(小波分析) end这些工具在实际项目中帮团队节省了大量重复工作。特别是智能参数选择功能让新人也能快速获得可用的分析结果而不用反复调整窗口参数。