做信号处理的兄弟应该都有这种体会单看两个传感器的自功率谱各自频谱都挺干净可一旦要判断它们之间谁先谁后、谁影响谁、在某几个频率点上的相位差到底是多少光靠自谱就完全抓瞎了。这时候就得请出互功率谱密度——也就是大家常说的 CPSD。我最早接触 CPSD 是在给一个钢结构平台做模态测试的时候当时想算激励点和响应点之间的频响函数被频响和相干函数来回折腾后来才把 CPSD 的定义、物理意义和数值计算全部理清。这篇文章不绕弯子直接把这套东西从头到脚拆开讲一遍适合刚接触随机信号分析的研究生也适合在实际测试里被频响函数和相干函数搞到头疼的工程师。1. CPSD 的定义与数学推导从互相关函数到互功率谱1.1 为什么随机信号不能直接做傅里叶变换我们平时习惯先把时域信号做傅里叶变换再讨论频谱。这个操作对确定性信号没有太大问题比如正弦信号、脉冲信号它们的能量有限变换结果在数学上是收敛的。但随机信号是完全另一码事。随机信号的一大特点是样本函数不满足能量有限条件。拿一段实测振动加速度信号来说你取无限长时间信号的能量积分是发散的也就是说它不满足傅里叶变换存在的绝对可积条件。这就意味着直接对单个样本做傅里叶变换得到的所谓“频谱”既不稳定也没有明确的物理意义——换一个样本结果就变了。所以处理随机信号的标准思路不是分析单个样本而是分析统计特性。自相关函数和互相关函数就是在时域上描述随机信号统计特性的工具。它们不关心具体的信号波形而是关心信号在不同时刻取值之间的关联程度。这样一来我们就不需要面对能量发散的问题相关函数通常满足傅里叶变换条件于是就可以顺理成章地进入频域。1.2 互相关函数的定义与 CPSD 的理论推导互相关函数描述的是两个信号 x(t) 和 y(t) 在相差延时 τ 的两个时刻上取值的相关性。对于平稳遍历随机过程可以用时间平均来替代集总平均R_xy(τ) lim(T→∞) (1/T) ∫(-T/2, T/2) x(t) · y(tτ) dt这个式子的物理含义很直白把 y(t) 在时间轴上往后挪了 τ然后看两个信号在每个时刻的乘积平均值。如果 x(t) 和 y(t) 之间存在某种固定的先后关系比如 y 总是比 x 晚 0.01 秒到达某个测点那么在 τ 等于这个延迟量附近互相关函数会出现明显的峰值。互功率谱密度 S_xy(f) 就定义为互相关函数 R_xy(τ) 的傅里叶变换S_xy(f) ∫(-∞, ∞) R_xy(τ) · e^(-j2πfτ) dτ这其实是维纳-辛钦定理在两信号场景下的推广。自功率谱密度 S_xx(f) 是自相关函数 R_xx(τ) 的傅里叶变换互功率谱密度 S_xy(f) 则是互相关函数的傅里叶变换。当 x(t) y(t) 时互功率谱密度自动退化为自功率谱密度所以 CPSD 可以看作 PSD 的推广形式。互相关函数在 τ 0 处的值等于两个信号的互功率也就是时间域内的交叉能量。互谱在整个频率轴上的积分正好等于这个值∫(-∞, ∞) S_xy(f) df R_xy(0)这个式子和 Parseval 定理的推广版是一致的。换句话说互谱在频域各频率点上描述了“共同能量”的分布而整个频段内的积分回到时域中等于两信号在当前时刻上的乘积均值。1.3 为什么 CPSD 是复数而不是实数初学的时候很容易在这里卡住自功率谱密度是实函数为什么互功率谱密度就变成复数了原因在于互相关函数 R_xy(τ) 一般不满足偶函数性质。自相关函数满足 R_xx(τ) R_xx(-τ)这是偶函数偶函数做傅里叶变换得到的是实函数。但互相关函数只满足 R_xy(τ) R_yx(-τ)并不等于 R_xy(-τ)。也就是说x 相对于 y 延迟 τ 的结果跟 y 相对于 x 延迟 τ 的结果是不一样的。一个信号领先另一个必然滞后这种“谁在前谁在后”的方向信息必须在频域里体现出来。于是互功率谱密度 S_xy(f) 的实部和虚部分别承载着不同的物理内涵实部代表两信号在同一频率上同相分量的贡献虚部则代表正交分量相位相差 90 度的贡献。由实部和虚部组合可以得到幅值和相位|S_xy(f)| sqrt(Re² Im²)θ_xy(f) arctan(Im / Re)这个相位 θ_xy(f) 就是两个信号在频率 f 处的相位差也是工程上最常用的信息之一。举个简单例子如果两个振动传感器测得的是同一根轴在两端传递过来的振动CPSD 的相位就能告诉你这个振动从 A 端传到 B 端需要多长时间——根据相位差和频率就能算出时间延迟 Δt θ / (2πf)。2. CPSD 的物理意义与工程解读幅值、相位、相干函数与频响估计2.1 CPSD 幅值到底在描述什么很多人第一次看到 CPSD 的计算结果会困惑为什么两个信号每个频率上的幅值都很大互谱的幅值却很小其实 CPSD 幅值不是两个自谱幅值的简单乘积它描述的是两个信号在该频率处“相干成分”的能量。假如两个信号在同一频率上完全独立、互不相关那么经过大量平均之后互谱的幅值会趋近于零。那些不相关的分量在平均过程中会被抵消掉就像噪声一样。只有当两个信号在某个频率上具有确定的相位关系时它们的乘积在平均后才会留下稳定的贡献。我在实际测试中的一个经验是用两个麦克风测同一个扬声器房间里如果有其他不相关的背景噪声源直接看每个麦克风的自谱背景噪声会明显抬高频谱但算两路信号的 CPSD 之后背景噪声的影响会大大降低。这是因为背景噪声在两路麦克风上是不相关的平均后贡献趋近于零而扬声器本身的声音在两路麦克风上具有很强的相关性和固定的相位关系在互谱中保留得很完整。这就是 CPSD 在噪声环境下能“提纯”信号相关成分的核心原因。幅值大小还和两信号幅度都有关系。如果其中一个信号特别小即使相干性极高互谱幅值也不会大。这个特性提醒我们CPSD 的幅值不能替代自谱去判断某个信号本身的能量大小它就是专门用来描述两个信号之间关联的指标。2.2 相干函数判断 CPSD 结果可信度的标尺光有 CPSD 还不够工程上通常还会算相干函数coherence。相干函数的定义是γ²_xy(f) |S_xy(f)|² / (S_xx(f) · S_yy(f))这个值在 0 到 1 之间。它表示在频率 f 处y(t) 的功率中有多大比例是由 x(t) 的线性作用引起的。如果 γ² 接近 1说明两个信号在这个频率上几乎完全线性相关如果接近 0说明它们几乎没有线性关系。工程中的经验判断大致是这样的在模态测试中激励点与响应点之间的相干函数在共振频率附近一般要求大于 0.9不然测出来的频响函数可信度就低。如果相干函数在某个频段普遍低于 0.5就要警惕以下几种情况一是系统存在严重的非线性比如结构间隙、摩擦或者大变形二是测量过程中存在显著的噪声干扰导致输出信号中包含了大量与输入无关的成分三是信号采集过程中出现了泄漏频谱发生了畸变四是激励能量不足导致响应信号的信噪比太低。还有一种情况要注意就是“虚假高相干”。在共振频率附近如果结构响应很大即使有一点泄漏或者噪声相干函数也可能被拉得很高。这时候不能盲目相信相干系数高就万事大吉还要结合相位曲线和模态置信准则一起判断。2.3 利用 CPSD 估计频响函数H1 与 H2 估计的取舍CPSD 最重要的工程应用之一就是频响函数估计。理论上的频响函数 H(f) 定义为输出响应 y(t) 的傅里叶变换与输入激励 x(t) 的傅里叶变换之比但实际测到的信号里总混有噪声直接做频谱相除结果会非常不稳定。更稳健的做法是利用互谱和自谱的比值。常用的有两类H1 估计H1(f) S_xy(f) / S_xx(f)H2 估计H2(f) S_yy(f) / S_yx(f)H1 假设噪声主要存在于输出端H2 假设噪声主要存在于输入端。大多数测试场景比如锤击模态测试和振动台激励试验激励信号的信噪比通常比响应信号高所以更常使用 H1 估计。它能把输出端不相关的噪声在平均过程中消掉从而获得更平滑的频响函数。这里要特别提醒一点用互谱算频响函数时CPSD 的归一化方式不太重要因为分子分母都有同样的因子约掉了。但如果你直接用 CPSD 的幅值去算传递率或者参与计算其他指标就必须确认归一化方式一致不然结果会差一个常数倍。3. CPSD 的计算方法与实操细节从理论和流程到代码实现3.1 基于 FFT 的数值计算流程实际工程中我们拿到的都是离散采样后的数字信号CPSD 的数值计算主要基于 FFT。标准的计算流程一般是这样的第一步对 x(t) 和 y(t) 做去直流处理也就是减去均值。不然直流分量会在频谱零频处产生一个巨大的尖峰还会通过窗函数泄漏到相邻频点干扰低频段的互谱结果。第二步选择分段长度。整段信号会被切成若干段每段长度记为 N。分段长度直接决定了频率分辨率Δf fs / N其中 fs 是采样率。如果段长 N 为 1024采样率 1024 Hz那么频率分辨率就是 1 Hz。想要更精细的频率分辨率就要增加段长但这会让可平均的段数减少方差增加这是一个绕不开的权衡。第三步对每一段施加窗函数。常见的窗有汉宁窗Hann、汉明窗Hamming、平顶窗Flat top等。对于连续随机信号最常用的是汉宁窗因为它能有效抑制频谱泄漏频率分辨率的损失也可以接受。第四步对加窗后的每一段数据做 FFT分别得到 X_i(k) 和 Y_i(k)。第五步计算每一段的互谱估计P_i(k) X_i(k) · Y_i^*(k)注意这里用了 Y 的复共轭。在 Python 的 numpy 中复数数组的共轭可以直接用 np.conj() 实现。第六步将所有段的互谱估计取平均S_xy(k) (1/M) Σ P_i(k)其中 M 是参与平均的段数。平均的目的是减小随机误差段数越多方差越小。第七步做归一化修正。这里最容易出问题因为不同软件、不同文献的归一化方式差异很大。最常见的做法是除以采样率 fs 和窗函数的功率修正因子。我给出一个经过验证的 Python 实现含注释便于直接复用import numpy as np def compute_cpsd(x, y, fs1024, nperseg1024, noverlapNone, windowhann): 计算 x 和 y 的互功率谱密度CPSD 参数 ---------- x, y : 1D array 输入信号二者长度必须一致 fs : float 采样率单位 Hz nperseg : int 每一段的FFT点数同时也决定了频率分辨率 noverlap : int 相邻分段的重叠点数默认为 nperseg // 2 window : str 窗函数类型支持 hann, hamming, boxcar 返回 ---------- freqs : 1D array 频率轴 cpsd : 1D complex array 单边互功率谱密度 n len(x) if noverlap is None: noverlap nperseg // 2 # 去直流 x x - np.mean(x) y y - np.mean(y) # 选择窗函数 if window hann: win np.hanning(nperseg) elif window hamming: win np.hamming(nperseg) elif window boxcar: win np.ones(nperseg) else: raise ValueError(未知窗函数) # 窗函数的功率修正因子 # 使用单边谱时幅值修正系数约为 2但功率谱使用的修正系数是 n * sum(win^2) / sum(win)^2 # 在计算互谱时我们使用能量修正方式保证总的积分功率一致 scale fs / (np.sum(win**2)) # 分段 stride nperseg - noverlap nseg (n - nperseg) // stride 1 cpsd_sum np.zeros(nperseg // 2 1, dtypecomplex) for i in range(nseg): start i * stride end start nperseg x_seg x[start:end] * win y_seg y[start:end] * win X np.fft.rfft(x_seg) Y np.fft.rfft(y_seg) cpsd_sum X * np.conj(Y) cpsd cpsd_sum / nseg * scale freqs np.fft.rfftfreq(nperseg, d1.0/fs) return freqs, cpsd这段代码的思路和我上文的流程完全一致。实际使用时需要注意如果拿这段代码和商业软件对比幅值上可能会差一个常数倍原因就是互谱归一化方式不同。不同软件对 PSD 和 CPSD 的默认归一化并不统一有的是除以 fs有的是除以分辨率 Δf还有的按周期图法处理。在做数据分析的时候最好固定使用同一套代码不要混用两套不同来源的算法去对比绝对值。3.2 Welch 平均法与关键参数的经验选择上面这段代码实际上就是 Welch 平均法。Welch 法的核心思想是通过分段、加窗、重叠、平均这一套组合拳在频率分辨率和估计方差之间寻找平衡。分段数量的计算公式是M floor((N - nperseg) / (nperseg - noverlap)) 1重叠率越高段数越多方差越小但相邻段之间的相关性也变大了所以边际收益递减。对于汉宁窗50% 重叠基本已经能获得足够的平均段数再往上加到 75%改善已经不明显计算量倒是增加了。我自己做测试数据分析时默认选择 50% 重叠。窗函数的选择也有一些讲究汉宁窗最常用适合宽带随机信号和大多数工程场景。汉明窗和汉宁窗类似但旁瓣略低主瓣稍宽差别不大。平顶窗的幅值精度高适合校准类测试但主瓣很宽频率分辨率差。矩形窗boxcar适合瞬态信号或整周期采样的情况用在连续随机信号上泄漏严重一般不建议。还有一个非常容易被忽视的细节窗函数的幅值修正和能量修正不是一回事。用汉宁窗时信号的能量被压缩了如果直接对加窗后的信号做 FFT计算出的功率谱幅值会偏低。如果只关心峰值幅度乘一个幅值修正系数 2 就够了但如果关心某个频段内的总功率就需要用能量修正系数也就是计算窗函数平方和与窗函数和的比值。在互谱计算中我的习惯是用能量修正方式把窗函数的影响尽量均匀分摊到频域各点。3.3 相位计算、相位谱展开与参考通道选择CPSD 的相位谱是另一个重要输出。相位角的计算方法是θ(f) arctan( Im[S_xy(f)] / Re[S_xy(f)] )在实际编程中一定要用四象限反正切 np.angle() 或者 atan2不要用普通的 arctan否则相位会被错误地折叠到 -π/2 到 π/2 之间丢失象限信息。相位谱经常会出现“锯齿状”跳变这是相位值被折叠到 -π 到 π 区间造成的。如果你关心的是传播延迟就需要对相位做解缠unwrap。解缠的本质是在相邻频点之间补偿 2π 的整数倍让相位曲线变成连续函数。在 Python 里直接用 np.unwrap() 即可但要注意解缠只在信噪比高的频段有意义噪声大的频段相位本身就是随机跳动的解缠不会改善结果。参考通道的选择也要留心。CPSD 是有方向性的S_xy(f) 和 S_yx(f) 的相位正好相反。在传递路径分析中选哪个信号作为参考直接决定了相位正负的解释。我的经验是优先选择信噪比高、物理意义明确的通道作为参考。比如做发动机振动传递路径分析时通常以激励源侧信号为参考这样计算出的互谱相位就是响应相对于激励的相位差便于解释。4. 常见问题与排查技巧CPSD 实操中的坑与对策问题现象可能原因排查与解决方式互谱幅值明显偏小窗函数能量损失未修正改用能量修正因子确保乘以 scale fs / sum(win^2)相干函数在高频段大幅波动信号信噪比不足平均段数太少增加重叠率、加长采样时间、提高激励能量相位谱出现严重抖动该频段两信号相干性低只看相干函数大于 0.8 的频段或增加平均次数零频附近出现巨大尖峰信号含直流分量做去均值处理必要时做高通滤波共振峰附近相干虚高但形状怪异泄漏或者双峰结构增加 FFT 点数或者改用多段加窗平均互谱结果随样本变化很大采样时间不够平均段数不足采集时长至少保证 20 个以上平均段表格里的问题我基本都踩过。其中最典型的是第一次用互谱算频响函数的时候发现 H1 估计出的共振峰幅值比预期低了接近一半。排查了很久才发现是窗函数的能量修正没有做对。当时用的汉宁窗但只做了幅值修正没有做能量修正导致频响函数的分子和分母虽然修正因子相同约掉了但拿互谱单独去和其他数据对比时出了问题。另一个常见坑发生在短样本分析中。如果采集数据只有一两秒钟采样率 1024 Hz分段长度设成 1024那就只有 1 到 2 段可以做平均互谱估计的方差非常大相干函数看起来也会很糟。这种情况再怎么优化算法也没用只能去补充采集时长或者接受较低的分辨率换取更多平均段数。还有一次做噪声源识别时两个测点间的 CPSD 相位始终不稳定。后来发现是其中一路传感器的相位响应有偏差两条测量链路的相位特性没有做校准。CPSD 的相位本质上反映的是两个通道之间的相对延迟如果传感器或者信号调理设备本身的相位响应不一致测出的相位差就不是真实的物理量了。所以在高精度测试中对两路传感器做相位一致性校准是非常重要的一步。5. 互谱计算中的几个特别提醒做 CPSD 计算时还有一个容易忽略的点采样率的整倍数频率和奈奎斯特频率处的处理方式。在单边互谱中直流分量和奈奎斯特频率对应的谱线是实数其他频点都是复数。工程中最常用的是单边谱也就是只取 0 到 fs/2 这一半幅值乘以 2直流和奈奎斯特频率除外。Python 的 rfft 自动处理了这种单边输出但 Multi 软件或 MATLAB 的习惯可能不同需要视情况确认。另外CPSD 的计算不要忘了加窗前的信号长度和重叠之间的关系。有人为了省事不设重叠直接分段这样带来的问题是如果使用汉宁窗每段数据在两端都被严重衰减大量有效信息被丢弃有效信号利用率不高方差自然就大了。这也是为什么 Welch 法一定要配合重叠使用的根本原因——不是省时间而是不重叠的话汉宁窗的代价太大了。如果你在计算互谱时使用的是时延估计、波束形成这类应用还要特别注意互谱矩阵的对称性处理。对于多个通道的情况完整的 CPSD 矩阵是共轭对称的也就是说 S_ij(f) S_ji^*(f)。在做阵列信号处理时利用这个对称性可以省一半计算量同时保证矩阵的正定性。最后强调一句互谱永远无法替代对问题的物理理解。我见过很多工具算出来的互谱和相干函数很好看但对照组不做、参考通道选错、传感器相位不一致最后结论完全是错误的。CPSD 是工具不是结论。数据采集阶段的质量控制比如传感器标定、通道匹配、采样同步远比后处理的算法技巧更值得花时间。行文至此CPSD 的整套思路基本讲完了。我个人在实际操作中最深的体会是刚上手的时候总觉得这是一堆公式和代码的堆叠做多了才发现真正考验人的永远是“这个频段的互谱为什么长这样”这种看似简单却需要综合判断的问题。下次再遇到相位谱抖动或者相干函数不合格先别急着调代码回到传感器布置、采样参数和数据质量上找原因往往比在算法层面死磕更有效。
