1. 从噪声中提取信号FFT实战入门第一次用Matlab做FFT分析时我对着屏幕上那堆乱七八糟的频谱线完全懵了。直到后来才发现原来真实世界的信号从来不会像教科书里那样干净。让我们从一个更接近实际的案例开始假设你从传感器采集到一组数据里面混合了50Hz工频干扰、100Hz的设备振动信号还有随机噪声。这时候FFT就像给信号做CT扫描能清晰看到各个频率成分的病灶。先来生成这个模拟信号Fs 1000; % 采样率1000Hz t 0:1/Fs:1-1/Fs; % 1秒时间轴 signal 0.7*sin(2*pi*100*t) sin(2*pi*50*t); % 主信号 noise 0.5*randn(size(t)); % 高斯白噪声 x signal noise; % 含噪信号画图看看时域波形figure; subplot(2,1,1); plot(t(1:200),x(1:200)); % 显示前200个点 title(含噪声的时域信号); xlabel(时间(s)); ylabel(幅值);现在运行FFT分析N length(x); X fft(x); f (0:N-1)*(Fs/N); % 频率轴 subplot(2,1,2); plot(f(1:N/2), abs(X(1:N/2))); % 单边谱 title(FFT频谱); xlabel(频率(Hz)); ylabel(幅值); grid on;你会看到频谱图上除了50Hz和100Hz的尖峰还有底部起伏的噪声基底。这就是典型的实际信号特征——确定性的周期信号叠加随机噪声。我经常用这个例子教新手理解频谱泄露现象如果信号长度不是信号周期的整数倍能量就会泄露到相邻频段形成虚假的频谱分量。2. 频谱分析的三大核心技巧2.1 采样参数设置的艺术采样率Fs和采样点数N直接影响频谱分辨率。去年调试电机振动时我就因为参数设置不当漏掉了一个关键谐波。这里有个实用公式频率分辨率Δf Fs/N最高分析频率Fmax Fs/2比如要分析500Hz以内的信号采样率至少1000Hz。如果想区分间隔1Hz的两个频率成分需要至少1秒的信号长度。举个例子Fs 2000; % 采样率提升到2000Hz t 0:1/Fs:2-1/Fs; % 2秒时长 x sin(2*pi*100*t) 0.5*sin(2*pi*101*t); % 两个相近频率 X fft(x); f (0:length(x)-1)*Fs/length(x); plot(f(1:200), abs(X(1:200))); % 局部放大2.2 加窗函数抑制频谱泄露直接截断信号相当于乘矩形窗会导致严重的频谱泄露。工程中常用汉宁窗(Hanning)或汉明窗(Hamming)来缓解win hann(N); % 生成汉宁窗 x_windowed x .* win; % 加窗 X_windowed fft(x_windowed);比较加窗前后的频谱figure; subplot(2,1,1); plot(f(1:N/2), abs(X(1:N/2))); title(不加窗频谱); subplot(2,1,2); plot(f(1:N/2), abs(X_windowed(1:N/2))); title(加汉宁窗频谱);加窗后主瓣会变宽但旁瓣泄漏显著降低。我在做轴承故障诊断时加窗使信噪比提升了15dB以上。2.3 零填充提高频谱可视性零填充不会增加真实分辨率但可以让频谱曲线更光滑N_pad 4*N; % 4倍零填充 X_pad fft(x, N_pad); f_pad (0:N_pad-1)*Fs/N_pad; plot(f_pad(1:N_pad/2), abs(X_pad(1:N_pad/2)));3. 功率谱估计从理论到实践3.1 周期图法直接计算PSD最简单的功率谱估计方法Pxx abs(fft(x)).^2 / (Fs*N); % 周期图 f (0:N/2)*Fs/N; plot(f, 10*log10(Pxx(1:N/21))); % 转换为dB title(周期图功率谱); xlabel(频率(Hz)); ylabel(功率谱密度(dB/Hz));但这种方法方差大我在分析风电齿轮箱信号时发现数据波动太大难以判断故障特征。3.2 Welch平均法稳定可靠的PSD估计Matlab的pwelch函数实现了Welch方法[pxx, f] pwelch(x, hann(256), 128, 1024, Fs); plot(f, 10*log10(pxx));关键参数说明256点汉宁窗128点重叠1024点FFT 实测这种设置对平稳信号效果很好但分析冲击信号时需要减小窗长。3.3 多窗谱估计应对复杂场景对于非平稳信号我推荐用pmtm函数多窗谱估计[pxx, f] pmtm(x, 3.5, 1024, Fs); % 3.5个DPSS窗 plot(f, 10*log10(pxx));4. 工业级应用案例解析4.1 电机振动分析实战去年帮某工厂分析电机异响采集到以下数据Fs 5000; % 5kHz采样 t 0:1/Fs:10-1/Fs; % 10秒数据 x 0.8*sin(2*pi*25*t) 0.3*sin(2*pi*137*t)... 0.1*randn(size(t)); % 25Hz转频137Hz轴承故障用Welch方法计算PSD[pxx, f] pwelch(x, hann(4096), 2048, 8192, Fs); plot(f(f200), 10*log10(pxx(f200))); % 观察200Hz以内发现137Hz处明显峰值对应轴承外圈故障特征频率。工厂更换轴承后振动值下降70%。4.2 音频信号处理技巧分析语音信号时我常用短时傅里叶变换(STFT)[s, Fs] audioread(speech.wav); spectrogram(s, hann(256), 128, 1024, Fs, yaxis);这个案例教会我对于时变信号单纯的FFT会丢失时间信息需要联合时频分析。4.3 嵌入式系统实现要点在DSP上部署FFT算法时要注意使用定点数运算减少资源占用预计算旋转因子节省计算时间合理选择FFT点数通常取2的整数幂一个C语言实现模板#include arm_math.h arm_cfft_instance_f32 S; arm_cfft_init_f32(S, 1024); // 初始化1024点FFT arm_cfft_f32(S, input, 0, 1); // 执行FFT