★ FEATURED ARTICLE 搞信号处理的兄弟应该都有这种经历手里抓着一串采样回来的波形数据脑子里第一个问题永远是“这玩意儿里面到底藏着哪些频率成分”。时域波形只能看到幅度怎么变但频率信息完全被埋在一条起起伏伏的曲线里。这时候傅里叶变换就该登场了而在工程落地环节MATLAB几乎成了频谱分析的事实标准工具。这个系列写了这么多篇从傅里叶变换的数学原理一路讲到工程实现今天这篇实战笔记就专门聊聊MATLAB在频谱分析里那些最常用、最容易被忽略也最值得抠细节的操作。要说明的是这篇不是简单的fft函数说明书而是把从信号生成、预处理、加窗、频谱绘制到问题排查的完整流程走一遍。适合刚接触信号处理仿真、正在做课程设计或者实际项目里需要分析振动、音频、电力谐波数据的读者有基础的朋友也能在这儿捡到一些平时文档里不写的细节。1. 为什么频谱分析非得用MATLAB从原理到工具的逻辑闭环1.1 傅里叶变换到底是干什么的先把最底层的逻辑盘清楚。傅里叶变换做了一件什么事它把一段时域信号拆解成一组不同频率、不同幅度、不同相位的正弦波叠加。你可以这样理解一锅高汤端上来你想知道里面放了哪几种调料、每种放了多少傅里叶变换就是那个帮你把调料成分和用量分析出来的“味觉分析仪”。数学上连续时间信号的傅里叶变换公式是X(f) ∫ x(t) * e^(-j2πft) dt但计算机没法处理连续积分我们手里的数据永远是一串离散采样点所以工程上用的其实是离散傅里叶变换DFTMATLAB里的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再到时频分析spectrogramMATLAB把信号处理领域的标准算法都封装成了可靠的函数极大地减少了重复造轮子的时间。要强调的是工具只是手段理解原理才是根本。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); % 采样点数这里N1000 Y fft(x); % 直接对x做FFT结果长度也是N跑完这行Y就是一个长度为N的复数数组。很多新手第一次拿到Y就直接plot(abs(Y))画出来的图横坐标是0到999的点数完全看不出频率信息。这就是频谱分析里第一个大坑FFT结果的下标并不是物理频率。Y的第一个元素Y(1)对应直流分量0HzY(2)对应频率fs/N 1HzY(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为中心、负频率在左正频率在右的双边谱。这时候就需要fftshiftY_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才能还原全能量。直流分量k1不需要乘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正弦幅值1120Hz正弦幅值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附近但由于噪声干扰直观判断并不精确。这就是频谱分析要解决的问题。接着做FFTY 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点数整数倍时比如fs1000、N100050Hz刚好是整数倍但如果信号是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));相位谱需要用anglephase 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 频率轴错了一半采样定理与奈奎斯特实战里最经典的翻车现场就是画出来的频谱在右边出现一条“对称的”镜像带。比如你做了一段fs1000Hz的信号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 GainCG定义为窗函数的平均值mean(w)。矩形窗的CG等于1所以不加窗时幅值天然准确。汉宁窗的CG约等于0.5所以加汉宁窗后所有频点的幅值都打了对折。要恢复真实幅值就得除以CGcg mean(w); mag_compensated mag_raw / cg;几种常用窗函数的相干增益如下窗函数相干增益CG主瓣宽度旁瓣衰减矩形窗1.0最窄约-13dB较差汉宁窗0.5较宽约-31dB海明窗0.54较宽约-41dB布莱克曼窗0.42最宽约-58dB很好这个表格能解释为什么有些测量仪器里幅值切换不同窗函数时读数会变。你在做频谱分析报告时必须注明用了什么窗以及是否做了CG补偿否则别人拿你的幅值数据无法复现。我的习惯是如果关心的是频率位置和泄漏抑制用汉宁窗如果关心幅值精度且信号频率恰好是整周期截断直接用矩形窗如果要分离两个频率很近的成分优先保证频率分辨率选主瓣更窄的窗。4.3 频谱图太乱分辨率与泄漏的取舍“频谱图太乱”本质上逃不脱分辨率和泄漏这两个问题。当两个频率成分间隔小于频率分辨率Δf时它们在频谱图上会融合成一个峰你根本分辨不出是两列频率。比如fs1000、N1000时Δf1Hz两个频率如果只差0.5Hz谱峰就会重叠。提高分辨率的唯一硬办法是增大N也就是采集更长的信号。但项目里信号时长往往受限这时候可以适当减小fs但前提是仍满足采样定理。注意补零不能提高分辨率前面已经强调过。泄漏会让谱峰底部长出旁瓣干扰相邻较弱频率成分的识别。此时加窗可以压旁瓣代价是主瓣变宽。所以窗函数选择永远是在做“分辨率”和“旁瓣抑制”的权衡。一个工程经验如果目标频率间隔较大比如几倍Δf优先用汉宁或海明窗如果信号本身就比较干净而且你希望看到最窄的谱峰就用矩形窗。我还要提一个容易被忽视的细节FFT点数不一定要等于信号长度。你可以对一段数据补零到2的整数幂比如4096、8192虽然不提升物理分辨率但可以提升FFT计算速度也让谱峰位置更圆滑便于人工判读峰值频率。现代MATLAB里fft即使长度不是2的幂也很快但大数据量场景我还是习惯补到2的幂。4.4 边界效应与栅栏效应栅栏效应这个名字听起来很玄实际意思是DFT只能在一系列离散频点上“张望”频谱如果真实信号频率落在两条频点之间你看到的峰值会低于真实幅值像是被栅栏挡住了一部分。比如fs1000、N1000频率网格是0、1Hz、2Hz……真实正弦是47.3Hz那么它的能量会被分摊到47Hz和48Hz两根谱线上每根都显得矮。应付栅栏效应的办法之一就是补零插值让离散频点更密峰值读得更准。另一个办法是使用频谱插值算法比如在峰值附近做抛物线拟合估计真实峰位和峰值。这个细节在测频类仪器里特别重要因为频率估计精度直接影响转速、振动特征提取的准确性。边界效应说的是我们在有限区间内做FFT相当于截断信号后把它当成周期信号处理。如果截断处信号不连续频谱就会出现Gibbs现象——峰根部的振荡纹波。加窗能缓解这个问题但无法根除因为信息丢失了就是丢了窗函数只是让丢失的代价更可控。5. 进阶玩法短时傅里叶变换与实时频谱监测5.1 非平稳信号用什么前面所有讨论都隐含一个假设信号的频率成分不随时间变化。但工程里很多信号是频率随时间变化的比如变频器的启动过程、语音的发音过程、轴承故障的瞬态振动。这时候做一次FFT只能得到一个“整个时间段平均”的频谱完全看不出频率演进过程。短时傅里叶变换STFT的思路非常朴素把信号切成一段段短窗对每一段分别做FFT然后把结果按时间轴排开。这样就得到一个时间-频率-幅值的三维谱图也就是常说的语谱图或频谱瀑布图。这段历史很有意思短时傅里叶变换最早广泛应用在语音信号分析里后来成为机械故障诊断、生物医学信号处理的标配工具。它解决的恰恰是普通FFT“只看全局、不看局部”的盲点。5.2 spectrogram的基本操作MATLAB里直接调用spectrogram就能完成STFTfs 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看看频谱有没有随时间变化。实践中踩过的坑比背一百个函数名管用得多。做信号处理这行永远是“时域看不出来的拉到频域看全局看不出来的拉到时间-频率平面看”。把这几板斧练熟你手里的数据会开口说话。
阅读完成 · 觉得有帮助？
