1. 从一段看起来像噪声的波形说起如果你做过振动监测、语音处理或者机械故障诊断大概率遇到过这种场景采集回来一段信号时域波形上完全看不出规律幅值忽大忽小肉眼看就是一团毛线。但把它丢进频谱分析里某些频率成分突然就冒出来了再进一步做倒谱又能看到一些频谱上根本读不出来的周期性结构。频谱和倒谱这两个工具就是干这个用的。这篇内容围绕频谱分析和倒谱分析这两个经典信号处理手段展开重点讲清楚它们各自在解决什么问题、背后的数学逻辑是什么、以及在MATLAB里怎么落地实现。适合正在做信号处理相关课题的学生、需要快速上手振动或声学分析的工程师以及学过傅里叶变换但不太清楚倒谱到底有什么用的朋友。我会尽量避开纯数学推导的堆砌用实际能跑的代码和能看懂的解释把这两个概念串起来。先给一个最直观的区分频谱告诉你信号里有哪些频率成分倒谱告诉你这些频率成分之间的间隔规律。前者是有什么后者是怎么排列的。很多故障特征——比如齿轮箱的啮合频率边带、轴承的故障冲击周期——恰恰藏在排列规律里而不是单个频率的幅值上。这就是为什么光看频谱不够还得上倒谱。下面我会从频谱的工程实现讲起再过渡到倒谱的原理和代码中间穿插几个我在实际项目中踩过的坑以及MATLAB里那些文档不会明说但很关键的参数选择经验。2. 频谱分析在MATLAB里的工程化实现2.1 为什么FFT不是调个函数就完事很多人第一次用MATLAB做频谱分析代码大概长这样Y fft(x); P abs(Y); plot(P);然后发现横坐标是点数不是频率幅值也不对低频部分还对称重复了一半。这不是MATLAB的问题是FFT本身的输出格式决定的。FFT输出的是复数序列包含正负频率且不做任何归一化。要得到有物理意义的单边幅值谱至少需要三步处理取单边、归一化、构造频率轴。我见过不少论文里的频谱图幅值标的是V或g但数值明显不对——要么忘了除以N要么忘了单边谱要乘2。这个细节在答辩或评审时很容易被追问值得一开始就做对。2.2 单边幅值谱的完整构造流程假设采样率fs信号长度N做N点FFT。标准做法如下function [f, A] singleSidedSpectrum(x, fs) N length(x); X fft(x); X X(1:floor(N/2)1); % 取前半段含直流和Nyquist A abs(X) / N; % 归一化 A(2:end-1) 2 * A(2:end-1); % 单边谱非直流非Nyquist分量乘2 f (0:floor(N/2)) * fs / N; % 频率轴 end这里有几个容易出错的点。第一floor(N/2)1这个索引范围当N为偶数时包含Nyquist频率点当N为奇数时不包含代码用floor统一处理是安全的。第二乘2的操作只针对非直流、非Nyquist的分量因为这两个分量在正负频率中各占一半本身就不需要补偿。第三归一化除以N是因为FFT的求和没有除以点数物理幅值需要还原。注意如果你的信号做了窗函数比如汉宁窗归一化系数还要除以窗函数的均值否则幅值会偏小。这个坑我在做声学测量时踩过频谱整体矮了一截排查了半天才发现是窗函数没补偿。2.3 频率分辨率与采样参数的取舍频率分辨率Δf fs/N这个公式大家都知道但实际用的时候经常矛盾想要高分辨率就得加长采样时间但设备内存或实时性又不允许。我的经验是先确定你关心的最小频率间隔再反推需要的采样时长。举个例子你要区分50Hz和52Hz两个成分间隔2Hz那Δf至少要做到1Hz以下采样时长T 1/Δf ≥ 1秒。如果采样率是1000Hz那N至少1000点。但如果你只是看整体频谱趋势不需要区分这么细N可以小一些配合平均来降噪。另一个常被忽略的是频率混叠。采样率必须大于信号最高频率的2倍这是奈奎斯特定理。但实际信号往往有高频噪声即使你关心的成分在低频高频噪声混叠进来也会污染频谱。所以硬件上加抗混叠滤波器或者软件上先做低通滤波是必要的步骤。2.4 实测中的频谱泄漏与窗函数选择频谱泄漏的根源是信号截断导致非整周期采样。FFT隐含假设信号是周期的如果截断处不连续频谱就会漏到邻近频点。解决办法是加窗让信号两端平滑过渡到零。不同窗函数的取舍窗类型主瓣宽度旁瓣衰减适用场景矩形窗最窄最差(-13dB)整周期采样、瞬态信号汉宁窗中等较好(-31dB)通用连续信号汉明窗中等好(-43dB)语音、振动平顶窗最宽最好幅值精确测量凯泽窗可调可调需要灵活权衡时选窗的核心逻辑是如果你更关心频率定位精度选主瓣窄的如果更关心幅值精度选旁瓣低的。平顶窗幅值精度最高但频率分辨率最差适合校准场景。汉宁窗是通用默认选择我在大多数振动分析里都用它。加窗之后记得做幅值补偿补偿系数是窗函数序列的均值。MATLAB里可以这样算w hanning(N); xw x .* w; CG mean(w); % 相干增益 A A / CG; % 补偿3. 倒谱当频谱不够用时该看什么3.1 倒谱解决的到底是什么问题频谱能把信号分解成频率成分但有一类信息它表达不好频谱上的周期性结构。最典型的例子是齿轮箱的振动信号。齿轮啮合会产生一个高频的啮合频率同时由于轴的转动啮合频率周围会出现一系列边带边带间隔等于轴的转频。在频谱上你看到的是一堆密集的峰很难直接读出边带间隔是多少。倒谱的做法是对频谱再取一次对数然后做一次傅里叶变换或逆变换。这样频谱上的周期性结构就会在倒谱域表现为一个明显的峰峰的位置直接对应周期。用生活化的类比频谱像是一首歌的和弦列表倒谱则是告诉你这些和弦每隔多久重复一次。数学上倒谱的定义有两种常见形式实倒谱C(τ) IFFT(log(|FFT(x)|))功率倒谱C(τ) IFFT(log(|FFT(x)|²))取对数的目的是把卷积关系变成加法关系。信号中的激励源和传递路径往往是卷积关系取对数后变成相加再变换就能分离。这是倒谱能分离源和路径的根本原因。3.2 实倒谱与复倒谱的区别和选择复倒谱保留相位信息定义为IFFT(log(FFT(x)))其中log是复对数。实倒谱只取幅值丢弃相位。两者用途不同实倒谱用于分析频谱的周期性比如边带间隔、谐波结构。计算简单结果实值最常用。复倒谱用于同态滤波、解卷积比如去除混响、分离卷积信号。需要处理相位展开问题实现复杂。大多数工程场景用实倒谱就够了。如果你只是想知道频谱上有没有等间隔的峰实倒谱完全胜任。复倒谱更多出现在语音信号处理的解卷积任务里。3.3 MATLAB实现倒谱的关键步骤实倒谱的MATLAB实现并不复杂但有几个细节决定结果对不对function [q, c] realCepstrum(x, fs) N length(x); X fft(x); logMag log(abs(X) eps); % 加eps防止log(0) c ifft(logMag); c real(c); % 取实部 % 构造倒谱横轴单位秒或采样点 q (0:N-1) / fs; % 通常只看前半段 c c(1:floor(N/2)1); q q(1:floor(N/2)1); endeps的加入很关键。如果频谱某个频点恰好为零log会返回-inf整个倒谱就废了。加一个极小值可以避免这个问题同时不影响结果。倒谱横轴的单位是时间秒但它的物理含义不是时间而是频谱上周期成分的重复周期。如果倒谱在τ0.01s处有峰说明频谱上每隔1/0.01100Hz就有一个峰。这个对应关系是理解倒谱的关键。3.4 倒谱中的低通滤波倒谱编辑倒谱有一个很实用的操作叫倒谱编辑cepstral editing本质上是在倒谱域做滤波。因为倒谱域中缓慢变化的频谱成分集中在低倒谱频率区快速变化的周期性成分集中在高倒谱频率区。如果你想去除传递路径的影响保留激励源特征可以在倒谱域把低倒谱部分置零再变换回频谱。具体操作% 倒谱编辑保留高倒谱成分 c_edit c; c_edit(1:cutoff) 0; % 去除低倒谱成分 % 变换回频谱 logMag_edit fft(c_edit); X_edit exp(logMag_edit); x_edit real(ifft(X_edit));这个操作在故障诊断里很有用。传递路径比如传感器安装位置、结构传递特性通常表现为频谱上的缓慢调制对应低倒谱而故障冲击的周期性对应高倒谱。通过倒谱编辑可以把两者分开。提示倒谱编辑的cutoff选择没有固定公式需要根据倒谱图上低倒谱区和高倒谱区的自然分界来定。我的经验是先画倒谱图找到第一个明显低谷的位置作为cutoff。4. 频谱与倒谱的联合诊断实战4.1 一个模拟的齿轮箱故障信号为了把两个工具串起来我构造一个模拟信号齿轮啮合频率500Hz轴转频20Hz故障导致啮合频率被转频调制同时叠加噪声。fs 5000; T 2; t 0:1/fs:T-1/fs; fc 500; % 啮合频率 fm 20; % 转频调制频率 x (1 0.5*cos(2*pi*fm*t)) .* cos(2*pi*fc*t); x x 0.3*randn(size(t)); % 加噪声这个信号的频谱会在500Hz处有一个主峰两侧间隔20Hz有边带。倒谱则会在τ1/200.05s处出现一个峰直接指示调制周期。4.2 频谱图上能读出什么、读不出什么对上述信号做频谱你会看到500Hz主峰和一系列边带。如果边带数量少、间隔大肉眼能读出来。但实际齿轮箱信号往往有多个啮合频率、多个轴的转频边带交织在一起频谱上就是一团密集的峰。这时候你很难判断哪个间隔对应哪个故障。另外频谱的幅值受传递路径影响很大。同一个故障传感器位置换一下频谱幅值可能差好几倍。但倒谱峰的位置不受传递路径影响因为传递路径是卷积关系在倒谱域表现为加性成分不改变峰的位置。这是倒谱相对于频谱的一个核心优势。4.3 倒谱峰的位置如何反推故障周期倒谱横轴τ对应的时间其倒数就是频谱上的周期间隔。如果倒谱在τ0.05s有峰说明频谱上每隔20Hz有周期结构。如果已知某个轴的转频是20Hz那这个峰就指向该轴的故障。实际操作中我会这样做先做频谱确认主要的频率成分和可能的调制现象。做倒谱找出所有显著峰的位置。把倒谱峰位置换算成频率间隔和已知的转频、啮合频率对照。如果某个倒谱峰对应的间隔与某轴转频吻合基本可以定位到该轴。这个流程在轴承故障诊断里尤其有效。轴承故障的特征频率内圈、外圈、滚动体通常不是直接出现在频谱上而是以冲击周期的形式存在倒谱能直接把这个周期抓出来。4.4 两个工具的分工与配合逻辑总结一下两者的分工分析目标频谱倒谱有哪些频率成分强弱频率成分的幅值强不直接频谱上的周期结构弱强边带间隔需人工读自动成峰传递路径影响敏感不敏感故障周期定位间接直接我的习惯是先频谱后倒谱频谱用来建立整体认知倒谱用来提取周期特征。两者结合基本能覆盖大多数旋转机械的故障诊断需求。5. 那些文档里不会写的参数陷阱5.1 采样率选高了反而坏事很多人觉得采样率越高越好其实不然。采样率高了同样的采样时长下数据量变大计算变慢更重要的是如果信号本身带宽有限过高的采样率会把大量噪声也采进来反而降低信噪比。我的经验是采样率取关心最高频率的2.5到4倍就够了留一点余量给抗混叠滤波器。5.2 倒谱的端点效应怎么处理倒谱做IFFT时如果频谱在两端不连续会出现类似频谱泄漏的端点效应。处理方法和频谱类似可以在取对数之前对频谱做平滑或者在倒谱域加窗。但更简单的做法是保证信号长度是2的幂次这样FFT和IFFT的效率最高端点效应也相对可控。5.3 对数底数的选择影响什么倒谱里取对数底数用e还是10还是2其实不影响峰的位置只影响幅值的缩放。因为换底公式只是一个常数因子而常数因子在IFFT后表现为冲激函数的缩放。所以不用纠结底数用自然对数log就行。5.4 平均与分段降低方差的关键单次FFT或倒谱的方差很大尤其是噪声环境下。标准做法是分段平均把长信号分成若干段每段做FFT然后对功率谱平均。倒谱也可以类似处理但注意倒谱是非线性变换有log平均要在功率谱域做然后再取log做倒谱。% 分段平均功率谱 nSeg 8; segLen floor(N/nSeg); Pavg zeros(1, segLen); for k 1:nSeg idx (k-1)*segLen (1:segLen); Xk fft(x(idx)); Pavg Pavg abs(Xk).^2; end Pavg Pavg / nSeg; % 再取log做倒谱 c real(ifft(log(Pavg eps)));这个细节在论文里经常被省略但实际做的时候如果不平均倒谱图会非常毛糙峰很难辨认。6. 从代码到结论几个我踩过的坑第一个坑是忘记去趋势。信号里如果有直流偏置或线性趋势频谱上会在0Hz附近出现巨大峰值倒谱也会被这个低频成分主导。做分析前先detrend一下能省很多事。第二个坑是倒谱横轴单位搞混。有人把倒谱横轴当成频率结果峰的位置全解释错了。记住倒谱横轴是时间单位秒它的倒数是频率间隔。第三个坑是窗函数加在倒谱之前还是之后。正确顺序是原始信号加窗做FFT取log做IFFT得到倒谱。不要在倒谱域加窗那会改变峰的形态。第四个坑是把倒谱峰和频谱峰一一对应。它们不是一一对应的。倒谱的一个峰可能对应频谱上的一组峰代表的是周期性而不是单个频率。这些经验在教科书里往往一笔带过但实际做项目时每一个都可能让你多花半天时间排查。希望这篇内容能帮你少走点弯路。频谱和倒谱这对组合用熟了之后面对再复杂的振动信号你至少知道该从哪里下手。