搞信号处理的兄弟应该都有这种经历:手里抓着一串采样回来的波形数据,脑子里第一个问题永远是“这玩意儿里面到底藏着哪些频率成分?”。时域波形只能看到幅度怎么变,但频率信息完全被埋在一条起起伏伏的曲线里。这时候傅里叶变换就该登场了,而在工程落地环节,MATLAB几乎成了频谱分析的事实标准工具。这个系列写了这么多篇,从傅里叶变换的数学原理一路讲到工程实现,今天这篇实战笔记,就专门聊聊MATLAB在频谱分析里那些最常用、最容易被忽略,也最值得抠细节的操作。
要说明的是,这篇不是简单的fft函数说明书,而是把从信号生成、预处理、加窗、频谱绘制到问题排查的完整流程走一遍。适合刚接触信号处理仿真、正在做课程设计或者实际项目里需要分析振动、音频、电力谐波数据的读者,有基础的朋友也能在这儿捡到一些平时文档里不写的细节。
1. 为什么频谱分析非得用MATLAB:从原理到工具的逻辑闭环
1.1 傅里叶变换到底是干什么的
先把最底层的逻辑盘清楚。傅里叶变换做了一件什么事?它把一段时域信号拆解成一组不同频率、不同幅度、不同相位的正弦波叠加。你可以这样理解:一锅高汤端上来,你想知道里面放了哪几种调料、每种放了多少,傅里叶变换就是那个帮你把调料成分和用量分析出来的“味觉分析仪”。
数学上,连续时间信号的傅里叶变换公式是:
X(f) = ∫ x(t) * e^(-j2πft) dt
但计算机没法处理连续积分,我们手里的数据永远是一串离散采样点,所以工程上用的其实是离散傅里叶变换(DFT),MATLAB里的fft函数做的就是这个离散版本,只不过它通过蝶形运算把复杂度从O(N²)压到了O(NlogN),所以叫“快速傅里叶变换”。
直接说结论:x(t)是时域里的信号,X(f)就是频域里的复函数,它的模|X(f)|对应频率成分的幅值大小,辐角对应相位。这两个量就是我们做频谱分析最关心的东西。
1.2 MATLAB在频谱分析里的核心优势
有些人会问,频谱分析用Python不行吗?用C写不行吗?都行,但MATLAB在信号处理场景有几个难以替代的优势。
第一,交互式调试太舒服了。你可以在命令行里随手生成一个信号,画图看看波形,觉得不对立刻改参数,不需要编译,不需要“写完代码-编译-运行-关闭”的循环。这对探索性分析极其重要,因为频谱分析本身就是个反复试参数的过程。
第二,矩阵运算是底层原生能力。fft的本质是矩阵变换,MATLAB对向量化运算的优化做得非常彻底。处理几百万个采样点的数组,MATLAB的fft在大多数机器上都是毫秒级响应,写起来还不必手动管理内存。
第三,信号处理工具箱和绘图体系非常完整。从窗函数生成(hann、hamming、blackman这些一行调用),到功率谱估计(pwelch、periodogram),再到时频分析(spectrogram),MATLAB把信号处理领域的标准算法都封装成了可靠的函数,极大地减少了重复造轮子的时间。
要强调的是,工具只是手段,理解原理才是根本。MATLAB让实现变得简单,但如果你不清楚频率轴怎么构造、幅值为什么要乘2、加窗之后幅度为什么会变,代码跑出结果你也不知道对不对。所以下面这些实操细节,比记住函数名重要得多。
1.3 一次完整频谱分析的流程拆解
完整的频谱分析绝对不是一句“Y = fft(x)”就结束的。我一般在项目里按下面这套流程走:
| 步骤 | 做什么 | 关键点 |
|---|---|---|
| 1 | 数据获取 | 确认采样率fs、采样点数N、数据是否有缺失或毛刺 |
| 2 | 预处理 | 去除直流分量、去趋势、剔除异常点 |
| 3 | 加窗处理 | 根据泄漏抑制需求选择窗函数 |
| 4 | 执行FFT | 确定FFT点数(是否补零) |
| 5 | 频率轴映射 | 把FFT结果的下标换算成真实物理频率 |
| 6 | 幅值换算 | 单边谱的幅值还原与窗函数增益补偿 |
| 7 | 可视化与分析 | 幅值谱/功率谱/相位谱的绘制与特征提取 |
这里面的每一步都有不少隐藏细节,接下来我会逐个拆开讲。
2. 上手必会:MATLAB频谱分析三大核心函数
2.1 fft:快速傅里叶变换的调用姿势
先把最基础的调用写出来。假设我们有一个采样率fs = 1000Hz的信号,时长1秒,由50Hz和120Hz两个正弦波叠加而成:
fs = 1000; % 采样率,单位Hz t = 0:1/fs:1-1/fs; % 时间向量,共1000个点 x = sin(2*pi*50*t) + 0.5*sin(2*pi*120*t); % 两个频率叠加 N = length(x); % 采样点数,这里N=1000 Y = fft(x); % 直接对x做FFT,结果长度也是N跑完这行,Y就是一个长度为N的复数数组。很多新手第一次拿到Y就直接plot(abs(Y)),画出来的图横坐标是0到999的点数,完全看不出频率信息。这就是频谱分析里第一个大坑:FFT结果的下标并不是物理频率。
Y的第一个元素Y(1)对应直流分量(0Hz),Y(2)对应频率fs/N = 1Hz,Y(k)对应的频率是(k-1)*fs/N。所以要从结果里读出物理频率,必须自己构造频率轴。
顺带说一句,fft函数的第二个参数很有用,比如Y = fft(x, NFFT)可以指定做多长点数的FFT。如果NFFT大于length(x),MATLAB会自动在数据末尾补零;如果NFFT小于length(x),MATLAB会自动截断数据。补零这种操作不是随便用的,后面会专门讲。
2.2 fftshift:把零频挪到中间
直接fft得到的结果,前半段对应正频率,后半段对应负频率。对于实信号来说,负频率部分和正频率部分是共轭对称的,信息量有冗余,所以大多数分析只用前半段就够了。
但遇到某些场景,比如画频谱瀑布图、做频域滤波、分析I/Q调制信号时,你更希望看到以0Hz为中心、负频率在左正频率在右的双边谱。这时候就需要fftshift:
Y_shift = fftshift(Y); f_shift = (-N/2 : N/2-1) * fs / N; plot(f_shift, abs(Y_shift));fftshift的本质就是把Y的前半段和后半段整体对调,让零频分量从数组开头挪到正中间。理解它的原理很重要:它没有改变任何数据的数值,只是重新排列了下标顺序。
2.3 频率轴构造与幅值还原
频率轴构造这部分,我要写得非常详细,因为这个细节直接决定了你的频谱图物理意义是否准确。
频率分辨率Δf = fs / N。这是频谱中两条相邻谱线之间的频率间隔,它只取决于采样率fs和FFT点数N。如果采样1秒1000个点,Δf就是1Hz;如果采样0.5秒500个点,Δf就是2Hz。分辨率这个概念后面还会反复提到。
完整的频率轴可以这样构造:
f = (0:N-1) * fs / N; % 双边谱完整频率轴而实际画单边谱时,通常只用前半部分:
half_N = floor(N/2); f_one = f(1:half_N); Y_one = Y(1:half_N);幅值还原是最多人踩坑的地方。fft结果的绝对值|Y(k)|是什么?它和真实幅值不相等。对于正弦信号,FFT结果中某根谱线的幅度等于该正弦幅值的N/2倍。要还原真实幅值,必须做两步操作:
第一步,除以N。因为FFT的定义里没有归一化,sum(x)被摊进了每个频点,所以要除以N才能让直流分量正确。
第二步,乘以2。单边谱只保留了正频率那一半,但实信号的能量在正负频率各占一半,所以正频率谱线的幅值要乘2才能还原全能量。
直流分量(k=1)不需要乘2,因为它没有对应的负频率分量。完整代码片:
mag = abs(Y_one) / N; mag(2:end) = mag(2:end) * 2; % 除直流外,其余乘2还原幅值 plot(f_one, mag);这段逻辑必须刻在脑子里。你要是偷懒直接用plot(abs(fft(x))),画出来的纵轴数值会大得离谱,而且没有物理意义。
3. 实操案例:从信号生成到频谱图绘制全流程
3.1 构造一个多频叠加信号
纸上谈兵没意思,直接来一个完整可复现的案例。假设我们有这样一段信号:50Hz正弦幅值1,120Hz正弦幅值0.5,再叠加一部分随机噪声,采样率设成1000Hz,采样时长还是1秒。
fs = 1000; t = 0:1/fs:1-1/fs; x = sin(2*pi*50*t) + 0.5*sin(2*pi*120*t) + 0.2*randn(size(t)); N = length(x);先画时域图:
subplot(2,1,1); plot(t, x); xlabel('时间 (s)'); ylabel('幅值'); title('含噪多频信号时域波形');从时域波形上,我们能看出信号大约在50Hz和120Hz附近,但由于噪声干扰,直观判断并不精确。这就是频谱分析要解决的问题。
接着做FFT:
Y = fft(x); f = (0:N-1) * fs / N; half_N = floor(N/2); f_one = f(1:half_N); mag = abs(Y(1:half_N)) / N; mag(2:end) = mag(2:end) * 2; subplot(2,1,2); plot(f_one, mag); xlabel('频率 (Hz)'); ylabel('幅值'); title('单边幅值谱');跑完这段,你会在50Hz和120Hz处看到两个清晰的尖峰,如果噪声不大,峰值高度会非常接近真实幅值1和0.5。这里我特别强调“接近”而不是“等于”,是因为噪声的能量会漏到各个频率,以及频谱泄漏的影响,下面立刻就会碰到。
3.2 加窗与补零:频谱泄漏的处理
如果你把上面这段代码跑出来,然后仔细观察谱线根部,会发现能量并不完全集中在50Hz和120Hz两个点上,而是向两侧蔓延开去,形成一条“拖尾”。这就是频谱泄漏。
泄漏的本质是:我们对无限长信号做了截断,相当于在时域乘了一个矩形窗,而矩形窗的频谱是sinc函数形状,它在频域会展宽信号原本的谱线。尤其是当信号频率不是FFT点数整数倍时(比如fs=1000、N=1000,50Hz刚好是整数倍,但如果信号是49.8Hz呢?),能量就会显著地“漏”到相邻谱线上。
解决泄漏的标准手段是加窗。常用窗函数有汉宁窗、海明窗、布莱克曼窗等,原理都是在时域让信号两端平滑衰减到0,从而压低截断造成的旁瓣。
w = hann(N); % 生成汉宁窗,长度N xw = x .* w'; % 逐元素相乘,注意维度 Yw = fft(xw);加窗之后,泄漏的旁瓣被大幅抑制,但也要付出代价:主瓣会变宽,也就是频率分辨率变差;同时因为窗函数把信号某些位置的幅值压低了,总能量减小,谱线幅值也会变小。因此加窗之后需要对幅值做补偿——也就是接下来要讲的相干增益校正。
补零是另一个常被误解的操作。有人以为补零能提高频率分辨率,这是个经典误区。fft(x, 8192)并不会让原本1秒的500点数据突然获得更高的物理分辨率,它只是在原有频谱上插值出更多频点,让峰值看起来更平滑、更接近真实谱峰。真正的分辨率提升只有两个途径:延长采样时间(增大N),或者减小采样率(减小fs,但受奈奎斯特限制)。补零适合用来找到谱峰的更精确位置,但不解决两个很近频率成分的分离问题。
3.3 幅值谱、功率谱、相位谱一次性搞定
一个正经的频谱分析报告,一般要给出三种谱:幅值谱、功率谱和相位谱。
幅值谱就是上面画的mag,表示每个频率成分的幅值大小。功率谱是幅值平方相关的量,工程上更常用的是功率谱密度(PSD),单位是W/Hz或dB/Hz,它把功率按照频率密度分布描述,方便比较不同带宽下的能量。
最简单的功率谱实现:
power = (abs(Y(1:half_N)).^2) / (N * fs); % 单边功率谱密度近似 power(2:end) = power(2:end) * 2;也可以更省心地调用pwelch,它能自动做分段平均、加窗,得到平滑的功率谱估计,适合随机信号分析:
[pxx, f_p] = pwelch(x, hann(256), 128, 512, fs); plot(f_p, 10*log10(pxx));相位谱需要用angle:
phase = angle(Y(1:half_N)); subplot(3,1,3); plot(f_one, unwrap(phase)); % unwrap用于消除相位跳变这里有个细节:相位谱在噪声存在时会非常混乱,因为噪声在所有频率都有分量,导致很小幅值的频点相位是随机的。所以实际项目里一般先设定一个幅值阈值,只分析幅值超过阈值的谱线相位。
一个完整的绘图布局就像这样:
figure; subplot(3,1,1); plot(f_one, mag); title('幅值谱'); subplot(3,1,2); plot(f_one, 10*log10(power)); title('功率谱密度 (dB/Hz)'); subplot(3,1,3); plot(f_one, unwrap(phase)); title('相位谱');4. 常见坑与排查实录:频谱分析翻车现场
4.1 频率轴错了一半?采样定理与奈奎斯特
实战里最经典的翻车现场,就是画出来的频谱在右边出现一条“对称的”镜像带。比如你做了一段fs=1000Hz的信号,FFT之后发现150Hz处有一个峰,但在850Hz处也有一个一模一样的峰。这不是信号真的有850Hz分量,而是FFT结果的后半段对应负频率,被很多人误画出来了。
解决方法是只取前半段分析,也就是前面提到的half_N处理。但要注意,即使只取前半段,频率轴的构造也必须正确:f = (0:half_N-1) * fs / N,最后一个频率点严格小于fs/2。
这里牵出采样定理这条铁律:采样率fs至少要是信号最高频率的2倍,否则高频成分会混叠到低频区间,形成假峰。我调试过一台设备,采样率设成500Hz去采300Hz的机械振动,结果频谱图上出现了200Hz的“鬼峰”,那个200Hz其实是300Hz混叠下来的。排查半天才定位到是采样率不足。
检查这个坑的方法很粗暴:把采样率大幅提高一倍,如果频谱图上某些峰的位置变了,或者新的峰出现了,多半是混叠。真正的信号频率不会随采样率变化而漂移。
4.2 幅值对不上?窗函数的相干增益补偿
加窗导致幅值衰减,这个问题几乎人人都会遇到。你以为信号是1V,加窗后幅值谱峰变成了0.5V,下意识以为是代码错了,其实不是。
每一种窗函数都有自己的相干增益(Coherent Gain,CG),定义为窗函数的平均值mean(w)。矩形窗的CG等于1,所以不加窗时幅值天然准确。汉宁窗的CG约等于0.5,所以加汉宁窗后所有频点的幅值都打了对折。要恢复真实幅值,就得除以CG:
cg = mean(w); mag_compensated = mag_raw / cg;几种常用窗函数的相干增益如下:
| 窗函数 | 相干增益CG | 主瓣宽度 | 旁瓣衰减 |
|---|---|---|---|
| 矩形窗 | 1.0 | 最窄 | 约-13dB(较差) |
| 汉宁窗 | 0.5 | 较宽 | 约-31dB |
| 海明窗 | 0.54 | 较宽 | 约-41dB |
| 布莱克曼窗 | 0.42 | 最宽 | 约-58dB(很好) |
这个表格能解释为什么有些测量仪器里幅值切换不同窗函数时读数会变。你在做频谱分析报告时,必须注明用了什么窗,以及是否做了CG补偿,否则别人拿你的幅值数据无法复现。
我的习惯是:如果关心的是频率位置和泄漏抑制,用汉宁窗;如果关心幅值精度且信号频率恰好是整周期截断,直接用矩形窗;如果要分离两个频率很近的成分,优先保证频率分辨率,选主瓣更窄的窗。
4.3 频谱图太乱?分辨率与泄漏的取舍
“频谱图太乱”本质上逃不脱分辨率和泄漏这两个问题。当两个频率成分间隔小于频率分辨率Δf时,它们在频谱图上会融合成一个峰,你根本分辨不出是两列频率。比如fs=1000、N=1000时Δf=1Hz,两个频率如果只差0.5Hz,谱峰就会重叠。
提高分辨率的唯一硬办法是增大N,也就是采集更长的信号。但项目里信号时长往往受限,这时候可以适当减小fs,但前提是仍满足采样定理。注意补零不能提高分辨率,前面已经强调过。
泄漏会让谱峰底部长出旁瓣,干扰相邻较弱频率成分的识别。此时加窗可以压旁瓣,代价是主瓣变宽。所以窗函数选择永远是在做“分辨率”和“旁瓣抑制”的权衡。一个工程经验:如果目标频率间隔较大(比如几倍Δf),优先用汉宁或海明窗;如果信号本身就比较干净,而且你希望看到最窄的谱峰,就用矩形窗。
我还要提一个容易被忽视的细节:FFT点数不一定要等于信号长度。你可以对一段数据补零到2的整数幂(比如4096、8192),虽然不提升物理分辨率,但可以提升FFT计算速度,也让谱峰位置更圆滑,便于人工判读峰值频率。现代MATLAB里fft即使长度不是2的幂也很快,但大数据量场景我还是习惯补到2的幂。
4.4 边界效应与栅栏效应
栅栏效应这个名字听起来很玄,实际意思是:DFT只能在一系列离散频点上“张望”频谱,如果真实信号频率落在两条频点之间,你看到的峰值会低于真实幅值,像是被栅栏挡住了一部分。比如fs=1000、N=1000,频率网格是0、1Hz、2Hz……真实正弦是47.3Hz,那么它的能量会被分摊到47Hz和48Hz两根谱线上,每根都显得矮。
应付栅栏效应的办法之一就是补零插值,让离散频点更密,峰值读得更准。另一个办法是使用频谱插值算法,比如在峰值附近做抛物线拟合,估计真实峰位和峰值。这个细节在测频类仪器里特别重要,因为频率估计精度直接影响转速、振动特征提取的准确性。
边界效应说的是:我们在有限区间内做FFT,相当于截断信号后把它当成周期信号处理。如果截断处信号不连续,频谱就会出现Gibbs现象——峰根部的振荡纹波。加窗能缓解这个问题,但无法根除,因为信息丢失了就是丢了,窗函数只是让丢失的代价更可控。
5. 进阶玩法:短时傅里叶变换与实时频谱监测
5.1 非平稳信号用什么
前面所有讨论都隐含一个假设:信号的频率成分不随时间变化。但工程里很多信号是频率随时间变化的,比如变频器的启动过程、语音的发音过程、轴承故障的瞬态振动。这时候做一次FFT只能得到一个“整个时间段平均”的频谱,完全看不出频率演进过程。
短时傅里叶变换(STFT)的思路非常朴素:把信号切成一段段短窗,对每一段分别做FFT,然后把结果按时间轴排开。这样就得到一个时间-频率-幅值的三维谱图,也就是常说的语谱图或频谱瀑布图。
这段历史很有意思:短时傅里叶变换最早广泛应用在语音信号分析里,后来成为机械故障诊断、生物医学信号处理的标配工具。它解决的恰恰是普通FFT“只看全局、不看局部”的盲点。
5.2 spectrogram的基本操作
MATLAB里直接调用spectrogram就能完成STFT:
fs = 1000; t = 0:1/fs:4-1/fs; % 构造一个频率随时间从50Hz线性升到200Hz的chirp信号 f0 = 50; f1 = 200; x = chirp(t, f0, t(end), f1); window = hann(256); noverlap = 200; nfft = 512; [S, F, T, P] = spectrogram(x, window, noverlap, nfft, fs);输出中的S是复数短时傅里叶系数,F是频率轴,T是时间轴,P是功率谱密度。直接画图:
figure; imagesc(T, F, 10*log10(P)); axis xy; xlabel('时间 (s)'); ylabel('频率 (Hz)'); colormap jet; colorbar; title('Chirp信号的时频谱');你会在时频谱上看到一条从50Hz平滑爬升到200Hz的亮带,这就是频率随时间变化的最直观可视化。这里的几个参数需要解释:
window长度决定时间分辨率和频率分辨率的平衡。窗越长,频率分辨率越好,但时间分辨率越差,因为窗内信号的局部时间被“糊”在一起了。窗越短,时间分辨率越好,但频率分辨率越差。noverlap是相邻窗重叠的样本数,重叠能减少窗边缘的影响,让谱图更平滑,代价是计算量增加。我一般设成窗长的75%左右。
5.3 从离线分析到实时系统的思路
离线分析和实时系统在频谱分析上有根本差别。离线时整个数据都在内存里,你可以反复迭代、选择最优参数;实时系统里数据是流水一样进来的,你必须在一个严格的时间预算内完成一帧数据的FFT并把结果用掉。
实时频谱分析的基本思想是分块处理。假定音频采样率48000Hz,每帧处理1024个点,那么每帧时长约21.3ms。系统每收到1024个点就做一次FFT,得到频谱,然后更新显示或送往下游特征提取模块。MATLAB里,如果你是在做算法原型验证,可以用AudioToolbox的dsp.AudioFileReader配合dsp.FFT实现流式频谱分析。
一个容易踩的实时坑是时延控制。如果你用了75%的重叠窗,那么每次窗口滑动256个点,帧和帧之间更新的间隔是5.3ms,但谱图显示的“当前时间”应该对应窗口中心位置,而不是窗口末尾。不然你看到的峰值频率总是滞后半帧,在实时故障诊断里这个滞后可能直接导致误判。
我建议做实时系统的人先离线把算法验证透了,再搬到流式框架里跑。因为实时系统一旦跑起来,你没法轻易停下来调试中间状态,而MATLAB的离线调试体验要好太多,很多“玄学”问题往往在离线环境下一眼就能看出来。
最后分享一点个人体会
这个系列写到这里,我最大的体会是:MATLAB的fft函数永远只是三行代码的事,真正的工程难度全在理解信号本身的物理特性,以及正确地选择分析参数。我调试过不少振动、音频、生物电信号,每次频谱图“看着不对劲”,最后追根溯源,要么是采样率设置不对,要么是窗函数没补偿,要么是频率轴没换算对,几乎都不是fft本身算错了。
如果你看完这篇想立刻动手,我建议你拿一段自己手里的真实数据跑一遍文中的完整流程:先画时域波形,再做单边幅值谱,再试试加窗和pwelch,最后用spectrogram看看频谱有没有随时间变化。实践中踩过的坑,比背一百个函数名管用得多。做信号处理这行,永远是“时域看不出来的,拉到频域看;全局看不出来的,拉到时间-频率平面看”。把这几板斧练熟,你手里的数据会开口说话。