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)/2+1); % 单侧频谱 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)/2+1); 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)对于CWT:
- Morse小波(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这些工具在实际项目中帮团队节省了大量重复工作。特别是智能参数选择功能,让新人也能快速获得可用的分析结果,而不用反复调整窗口参数。