当年学随机过程的时候很多同学都会有同一个疑问确定性信号好好做傅里叶变换就行了为什么随机过程要单开一章讲频域分析功率谱这个名词又是什么意思等到真正做工程、处理噪声和振动信号时才发现这个工具不是“可选”的而是绕不开的。今天这篇就把随机过程的频域分析这条线完整梳理一遍重点讲清楚功率谱是怎么来的、维纳-辛钦定理为什么是核心以及实际使用中那些容易踩的坑。这篇内容适合正在学随机过程、信号与系统或者刚开始做信号处理相关课题的同学参考。我会把定义、推导、性质、应用和实操估计都串起来尽量讲得直白点毕竟当初我也是被“功率谱密度”这个拗口名字折磨过的人。在大多数高校的随机过程课程安排里相关函数和功率谱这部分内容一般都在排队论这类应用模块之前因为后面很多分析都要用到频域统计特性这一块不吃透后面会很吃力。1. 为什么随机过程要专门做频域分析1.1 确定性信号那套傅里叶工具套不上随机过程先回忆一下确定性信号的处理方式。一个能量有限信号 (x(t))只要满足绝对可积或平方可积条件就能做傅里叶变换得到频谱 (X(f))然后可以看幅频特性、相频特性还能做滤波、卷积等等。这一套在通信、控制、电路里用得非常顺手。但随机过程一上来就遇到问题。随机过程 (\xi(t)) 的每一个样本函数 (x_i(t)) 都不同而且一般持续无限长时间多数样本函数既不绝对可积也不平方可积傅里叶变换根本不存在。就算你强行对某个样本做截断再求变换不同样本得到的频谱差别也很大那这个“频谱”到底代表谁没有统计意义。还有一个更本质的问题随机信号的相位是随机的而功率、能量这些物理量只看幅值不看相位。所以直觉上我们关心的不是每个样本的具体频谱长什么样而是“平均意义上这个随机过程的能量或者功率在频域是如何分布的”。频率成分在统计意义上有多强这才是随机过程频域分析要回答的问题。1.2 换个思路从能量谱转向功率谱确定性信号如果能量有限可以用能量谱密度 (|X(f)|^2) 来描述能量在频域的分布。但随机过程通常是功率有限而能量无限就像家用电器一样你问一台空调“总能量是多少”没有意义因为它一直开着能量一直在累积但问“平均功率是多少”就很有意义。类比一下就清楚了。热水器的功率是2000W意思是它单位时间消耗的能量是2000焦耳每秒。至于这2000W里面有多少是加热元件贡献的、多少是控制电路贡献的这是两种不同频率机制的功率分配问题。功率谱密度想表达的正是这个整个随机过程的平均功率按频率拆开每个频率附近贡献多少。对于平稳随机过程来说统计特性不随时间平移变化平均功率是一个稳定的值不会因为你取哪一段样本而变化。这个特性保证了“平均功率在频域的分布”是一个可以定义的、稳定的量。于是问题就转化为怎么从随机过程本身出发严谨地把这个频域分布构造出来。2. 功率谱密度到底是怎么定义的2.1 从帕斯瓦尔定理出发先看确定性信号的帕斯瓦尔定理。对能量有限信号 (x(t))它的傅里叶变换是 (X(f))那么信号总能量可以用时域和频域两种方式计算[ E \int_{-\infty}^{\infty} |x(t)|^2 dt \int_{-\infty}^{\infty} |X(f)|^2 df ]这里 (|X(f)|^2) 就是能量谱密度。它的物理含义很直接在频率 (f) 附近无限小频带内包含的能量除以带宽就是该点的能量密度。随机过程没法直接套用这个公式因为样本函数不是能量有限的。但是我们可以把样本函数“截断”成一段有限长度的信号来处理就像做实验时只记录一段时间内的数据一样。截断后的信号能量有限可以做傅里叶变换然后计算这段信号的平均功率最后让截断时间趋于无穷大看平均功率是否收敛到一个稳定值。2.2 截断、取极限、再取期望三步定义平均功率谱具体操作是这样的。取随机过程的一个样本函数 (x(t))把它截断在区间 ([-T, T]) 内得到[ x_T(t) x(t) \cdot \text{rect}\left(\frac{t}{2T}\right) ]对这个截断信号做傅里叶变换得到 (X_T(f))。根据帕斯瓦尔定理截断信号的总能量为[ \int_{-\infty}^{\infty} |x_T(t)|^2 dt \int_{-\infty}^{\infty} |X_T(f)|^2 df ]左边等于 (\int_{-T}^{T} |x(t)|^2 dt)。把总能量除以时间长度 (2T)就得到这段样本在 ([-T,T]) 上的平均功率[ P_T \frac{1}{2T} \int_{-T}^{T} |x(t)|^2 dt \int_{-\infty}^{\infty} \frac{|X_T(f)|^2}{2T} df ]现在让 (T \to \infty)左边如果收敛就是该样本函数的平均功率。右边被积函数在极限下的值就叫做样本的功率谱密度。但这里有个问题单个样本函数的极限仍然可能带有随机性不同样本结果会不同。所以真正有统计意义的定义需要先把 (|X_T(f)|^2) 对随机过程取期望再取极限[ S_{\xi}(f) \lim_{T \to \infty} \frac{E\left[ |X_T(f)|^2 \right]}{2T} ]这就是功率谱密度的定义式。它描述的是随机过程的平均功率在频率轴上的分布密度单位是“功率/Hz”。注意分母里的 (2T) 是截断时间长度不是别的什么东西它把能量密度转换成功率密度这一步是整个定义的灵魂。注意定义里要求极限存在。对宽平稳随机过程来说只要相关函数满足一定条件后面会讲这个极限就是存在的。这也就是为什么平稳性这么重要——没有平稳性平均功率本身都不稳定频域分析就更无从谈起了。3. 维纳-辛钦定理相关函数与功率谱的桥梁3.1 定理内容与证明思路定义式虽然严谨但真让你按这个式子去计算某个过程的功率谱你会崩溃的因为每个样本都要做傅里叶变换还要取期望和极限计算量太大了。还好有维纳-辛钦定理这个终极武器。维纳-辛钦定理说的是对于宽平稳随机过程 (\xi(t))如果自相关函数 (R(\tau)) 绝对可积那么功率谱密度和自相关函数构成一对傅里叶变换对[ S_{\xi}(f) \int_{-\infty}^{\infty} R(\tau) e^{-j2\pi f \tau} d\tau ][ R(\tau) \int_{-\infty}^{\infty} S_{\xi}(f) e^{j2\pi f \tau} df ]换句话说时域里算自相关函数频域里就有了功率谱两者一一对应。这个定理的价值怎么强调都不过分。没有它你得吭哧吭哧地做大量样本统计有了它你只要把相关函数算出来做一次傅里叶变换就完事。证明的思路也不复杂把定义式里的期望和积分交换顺序就行。回顾一下[ E\left[ |X_T(f)|^2 \right] E\left[ \int_{-T}^{T} \int_{-T}^{T} x(u) x(v) e^{-j2\pi f (u-v)} dudv \right] ]交换期望和积分注意到 (E[x(u)x(v)] R(u-v))然后换元 (\tau u - v)整理之后会得到[ S(f) \lim_{T \to \infty} \int_{-2T}^{2T} \left(1 - \frac{|\tau|}{2T}\right) R(\tau) e^{-j2\pi f \tau} d\tau ]当 (R(\tau)) 绝对可积时(\frac{|\tau|}{2T}) 这一项在极限下趋近于0加权系数趋近于1积分范围变成整个实轴于是自然得到维纳-辛钦定理。这个过程其实就是傅里叶变换里的Fejér核逼近思想从工程角度记住结论就够了但理解推导有助于你明白定理成立的前提条件。3.2 均值非零、有周期分量怎么办实际工程中经常遇到均值不为零的过程比如传感器的直流偏置、旋转机械的基频振动分量。这时候自相关函数里会包含常数项或者周期分量导致 (R(\tau)) 不绝对可积维纳-辛钦定理形式上会遇到困难。解决办法是引入冲激函数Delta函数来扩展理论框架。均值 (\mu) 会在频率零处产生一个冲激 (\mu^2 \delta(f))周期性分量会在对应频率处产生冲激。这样处理之后维纳-辛钦定理仍然成立只是功率谱中包含了离散的谱线成分。在实际操作中我们通常会把这类分量先分离掉或者用高通滤波滤除直流后再分析。你要是直接拿着含直流的信号算功率谱零频附近会有一个巨大的尖峰把其他频率成分的细节都压得看不见。我在处理振动数据时一般会先做去均值detrend这个习惯能帮你省掉很多麻烦。4. 功率谱的性质与典型随机过程4.1 三条基本性质功率谱密度作为“密度”有几个非常重要的基本性质这些性质决定了它的使用方式。第一非负性(S(f) \ge 0)。这个看起来理所应当但其实是很强的约束。(\frac{E[|X_T(f)|^2]}{2T}) 是非负的取极限也非负所以功率谱永远不可能出现负数。这跟一般傅里叶频谱不一样频谱的实部和虚部可正可负但功率谱是能量密度的统计平均绝不可能是负的。第二偶对称性实随机过程的功率谱满足 (S(f) S(-f))。这是因为自相关函数是偶函数而偶函数的傅里叶变换是实偶函数。这带来一个重要结论正负频率贡献的功率是相等的所以工程上经常只看正频率部分然后把谱值乘以2得到“单边谱”。第三积分等于平均功率[ P_{\text{avg}} \int_{-\infty}^{\infty} S(f) df ]同时维纳-辛钦定理给出 (R(0) \int S(f) df)而 (R(0) E[|\xi(t)|^2])两个式子完全一致。这说明自相关函数在零点取值就是平均功率功率谱曲线下的总面积只要等于这个值你的谱计算就基本没做错。4.2 白噪声、有色噪声、窄带过程掌握了功率谱之后就可以用频域视角重新认识一些典型的随机过程了。这些典型模型在通信和信号处理里随处可见。白噪声是最理想化的模型。它的自相关函数是 (R(\tau) \frac{N_0}{2}\delta(\tau))功率谱密度在整个频率轴上平坦等于 (\frac{N_0}{2})。就像白光包含所有可见光频率一样白噪声在所有频率上均匀贡献功率。但要注意理想白噪声的平均功率是无限的这在实际物理系统里不可能存在。真实的热噪声、散粒噪声只是在很宽的频带内近似平坦带宽远大于系统通带工程上就可以当作白噪声处理。有色噪声的功率谱不平坦比如 (1/f) 噪声在低频处功率特别大高频处衰减。电子器件里的闪烁噪声就是典型。还有一种常见模型是一阶低通噪声功率谱形如 (S(f) \frac{A}{1 (f/f_c)^2})这种谱对应的是相关函数指数衰减的过程 (R(\tau) \sigma^2 e^{-|\tau|/\tau_c})其中 (\tau_c) 是相关时间(f_c) 是3dB截止频率两者满足 (f_c 1/(2\pi\tau_c))。窄带过程的功率谱集中在一个中心频率 (f_0) 附近的窄频带内频带宽度远小于中心频率。通信中的带通信号、结构振动中的某一阶模态响应都属于窄带过程。它的自相关函数表现为快变振荡乘以慢变包络的形式这正是窄带信号的时域特征。过程类型功率谱特征自相关函数特征典型来源理想白噪声全频段平坦(\delta(\tau)) 冲激理论模型热噪声近似低通有色噪声低频平坦高频衰减指数衰减RC电路噪声、湍流窄带过程集中在 (f_0) 附近振荡衰减包络通信载波、结构共振含周期分量过程平坦背景上叠加谱线含不衰减周期项旋转机械振动噪声4.3 工程中更常用单边谱因为实随机过程功率谱关于零频偶对称正负频率各贡献一半功率工程仪器里几乎不会显示负频率所以单边谱 (P(f)) 定义为[ P(f) 2S(f), \quad f \ge 0 ]直流分量所在频率 (f 0) 处不乘2因为零频只有一个点不存在正负两个方向的对称副本。这个换算看着简单但特别容易出错。很多人在用频谱仪或者软件里的PSD功能时对不上数值往往就是单双边谱搞混了。我的习惯是在正式处理数据前先算一个已知功率的正弦信号做标定确认整个链路里使用的是单边还是双边约定这样能避免后续结果的单位误差。5. 在工程和实验里怎么用功率谱5.1 从噪声里找周期信号功率谱最典型的应用场景就是在强噪声背景里发现微弱周期信号。单个周期信号在时域里可能完全淹没在噪声中肉眼根本看不出规律。但是周期信号的全部功率集中在单一频率上而白噪声的功率均匀散布在整个频带上。想象一下总功率为1的混合信号里如果周期信号只占0.1的功率噪声占0.9那么信噪比只有-9.5dB左右。在时域上看就是一条毛糙的曲线看不出任何周期。但做功率谱之后这0.1的功率全部压在一个频率分辨率带宽 (\Delta f) 内谱峰值会远远高出噪声基底。频带越窄积累效果越明显谱峰和噪声基底之间的比值越大。这就是频谱分析在故障诊断、语音识别等领域如此通用的原因。我做轴承故障检测的时候就经常用这个方法。正常轴承的振动谱是相对平坦的背景加少量谐波一旦外圈出现裂纹会激发某个固有频率附近的调制边带功率谱上会出现特征频率及其倍频。很多早期故障在时域波形上完全看不出来但功率谱上已经有非常明显的谱线变化。5.2 用FFT估计功率谱的正确姿势数字信号处理里估计功率谱最常用的方法是Welch法也就是“分段加窗 FFT 平均”。先说为什么不能直接对整个序列做一次FFT然后取模平方。单次FFT得到的周期图方差非常大而且这个方差不会随着数据长度增加而下降这是个反直觉的现象——你以为数据越长越准其实周期图还是那么毛糙。Welch法的思路很简单把长度为 (N) 的数据分成若干重叠段每段做FFT取模平方得到该段的周期图然后把所有段的周期图平均。因为不同段的噪声是相对独立的平均之后方差大约除以段数。你牺牲了一点频率分辨率换来了谱估计方差的显著下降这个交易在绝大多数工程场景里都是值得的。MATLAB里用现成函数就能实现fs 1000; % 采样率 t 0:1/fs:10; % 10秒数据 x sin(2*pi*50*t) 0.5*randn(size(t)); N 1024; % 每段长度 overlap 50; % 重叠百分比 w hann(N); % 汉宁窗 [pxx, f] pwelch(x, w, overlap, N, fs); plot(f, 10*log10(pxx))5.3 参数怎么定频率分辨率和平均次数的权衡Welch法里有几个关键参数需要根据自己的场景来定它们互相制约不能拍脑袋随便选。首先是每段长度 (N)它直接决定频率分辨率。分辨率 (\Delta f f_s / N)段越长频域越精细。如果你要分辨两个相差5Hz的频率分量采样率是1000Hz那么 (N) 至少要200点实际使用最好取256或512以上。但段越长整段数据能分出的段数越少平均次数下降方差又会变大。然后是重叠率。相邻段完全不重叠当然可以但数据利用率低。50%重叠是常用折中75%重叠在窗函数衰减较快时能进一步降低方差但计算量也会上升。我一般是先用50%看谱线还是毛糙就提到75%。最后是窗函数。矩形窗频率分辨率最好但旁瓣泄漏最大汉宁窗是通用选择旁瓣衰减快频率分辨率损失可以接受如果动态范围要求高可以用布莱克曼窗。参数作用典型取值注意事项段长 (N)频率分辨率按频谱细节需求定太短会把靠得近的谱峰糊在一起重叠率方差抑制50%-75%重叠太高计算量大收益递减窗函数抑制谱泄漏Hann默认矩形窗泄漏大不要用于强干扰场景平均段数方差抑制尽量多与段长互相制约实操心得如果你的目标是找周期信号先用矩形窗争取最大分辨率确认谱峰位置后再加窗做精细分析如果目标是比较平稳的随机噪声谱形直接用汉宁窗加高重叠率把谱线磨光滑最重要。两种场景的参数取向完全相反别一套参数走天下。6. 常见误区与排查技巧6.1 别把单次FFT的周期图当成功率谱这是新人最容易犯的错。直接对一段数据做FFT取模平方画出来发现谱线又密又毛糙密密麻麻全是“刺”这是周期图的典型表现。问题的根源在于周期图估计的方差很大它的标准差和信号本身的功率谱值在同一个量级而且不管你数据多长方差都降不下来。解决方式就是平均。把一个长序列切成很多段分别算周期图再取平均。你可能会问平均之后频率分辨率不是下降了吗确实下降了但这是用分辨率换稳定性的合理代价。工程上绝大多数情况下平滑稳定的谱比精密但抖动的谱有用得多。6.2 频率分辨率与谱泄漏两个容易缠在一起的问题谱泄漏是因为截断带来的。有限长数据等于把无限长信号乘了一个矩形窗矩形窗的频谱有很高的旁瓣于是强信号的能量会泄漏到邻近频率上去。两个幅度差别很大的信号弱信号很可能被强信号的泄漏旁瓣淹没。解决谱泄漏的办法是加窗。加Hann窗后主瓣稍微变宽但旁瓣衰减快得多强信号对弱信号的压制会小很多。我的经验是做工程谱分析时默认加窗只有在校准或需要严格频点读数时才考虑矩形窗。6.3 单双边谱、单位换算能出错的地方太多了最后说说单位的问题。功率谱密度PSD的单位是 (V^2/\text{Hz}) 或者 (W/\text{Hz})而用FFT直接算出来的幅值谱单位是 (V)取模平方后如果不做归一化数值上差了采样率和数据长度的倍数关系。用MATLAB的pwelch时默认输出已经是单边PSD单位是 (V^2/\text{Hz})可以直接用。但如果你自己写FFT做周期图就得手动归一化N length(x); X fft(x .* hann(N)); Pxx 2 * abs(X).^2 / (fs * sum(hann(N).^2)); f (0:N/2-1) * fs / N;注意分母里的fs * sum(w.^2)这是把窗函数能量和采样率一起归一化缺一个都会导致数值偏差。直流成分对应的点不要乘2做单边化时只对正频率部分乘2。我每次写完自制的功率谱估计代码都会先用一个已知幅度和频率的正弦信号做验证检查峰值处 (A^2/4) 是否正确时域正弦信号在单边功率谱上的峰值功率理论上是幅度平方的四分之一。这个验证只花两分钟但能避免你觉得程序写对了、其实单位全错的情况。回到功率谱这个概念本身我个人的体会是它本质上就是随机过程在频域的“统计指纹”描述的是平均功率随频率的分布而不是某一次样本的具体频谱。理解这一点之后再看维纳-辛钦定理、再看Welch估计就会觉得顺理成章。刚开始学的时候我一直在纠结为什么不能用傅里叶变换直接处理随机过程后来明白“统计平均”才是随机过程频域分析的灵魂。最后再分享一个小技巧遇到一个陌生的随机过程先画出它的自相关函数和功率谱从两个域对比着看很多性质一下子就清楚了——相关函数的宽度对应谱的带宽衰减速度对应谱的平滑程度这两个视角互相印证学起来会快很多。