互功率谱密度(CPSD)从定义到Python实战：幅值、相位与相干分析

搞振动测试、声学测量或设备故障诊断的人对自功率谱密度PSD一定不陌生。但很多时候问题不是“这个信号有多强”而是“两个信号之间有什么关系”。比如我测了轴承座和基础的两个加速度响应想判断振动是从哪个方向传过来的再比如做声源识别时布置了一排传声器想知道不同位置的声压在哪些频率上高度相关。这种场景下自谱帮不上忙你需要的是互功率谱密度Cross Power Spectral DensityCPSD。互谱本质上是两个信号互相关函数的傅里叶变换它用一个复数同时给出幅值和相位幅值说明两个信号在某个频率上共同变化的强弱相位说明一个信号相对另一个信号的滞后程度。这篇文章不绕弯子直接从定义、推导、物理意义讲到计算方法最后用一段可运行的 Python 代码演示完整流程。做信号处理的新手、用着商业软件却不太清楚背后原理的工程师都可以照着过一遍。1. 从互相关函数到互功率谱先把定义和推导理清楚1.1 从互相关出发而不是从“功率”出发互功率谱密度不是凭空定义的它的根在时域的互相关函数。对于两个零均值的随机信号 (x(t)) 和 (y(t))互相关函数定义为[ R_{xy}(\tau) E\left[x(t\tau) y^*(t)\right] ]这里的 (*) 表示共轭对实信号可以忽略。拿人话讲就是把 (y) 固定住、把 (x) 在时间轴上左右平移 (\tau)然后看两个信号的平均相似程度。如果某个平移量下两者形态最接近互相关就会出现峰值这个峰值对应的 (\tau) 就是两个信号之间的时间延迟。这个思路在时延估计里很直观但实际问题里信号往往是多个频率成分叠加的时域互相关看起来像一团乱麻峰值不好找也不容易说清楚“哪个频率成分造成了这个延迟”。所以要把互相关函数变到频域。对 (R_{xy}(\tau)) 做傅里叶变换得到的就是互功率谱密度[ S_{xy}(f) \int_{-\infty}^{\infty} R_{xy}(\tau) e^{-j 2\pi f \tau} d\tau ]这是定义式。但要理解它为什么能写成频域信号乘积的期望需要把傅里叶变换代进去推一步。假设 (X(f)) 是 (x(t)) 的傅里叶变换(Y(f)) 是 (y(t)) 的傅里叶变换在平稳随机过程框架下可以证明[ S_{xy}(f) E\left[X(f) Y^*(f)\right] ]注意我这里写的是 (X(f) Y^(f))也就是对 (Y) 取共轭。这样当 (yx) 时退化成自功率谱密度 (S_{xx}(f)E[|X(f)|^2])是一个非负实数。如果换成 (X^(f)Y(f))得到的互谱相位方向会相反后面分析时很容易搞反这点我会在后文反复强调。[ S_{xy}(f) |S_{xy}(f)| e^{j \angle S_{xy}(f)} ]这就是 CPSD 的核心幅值 (|S_{xy}(f)|) 和相位 (\angle S_{xy}(f))两者都是频率 (f) 的函数。互相关是实数但互谱是复数因为它保留了相对相位信息。1.2 有限长信号与谱估计的第一道坎理论上互谱需要无限长时间信号实际工程中永远只能拿到一段有限长数据所以 (S_{xy}(f)) 只能估计、不能精确计算。用有限长数据做傅里叶变换时相当于对无限长信号乘了一个矩形窗频域会和窗函数的频谱做卷积造成频谱泄漏。对于互谱来说泄漏的后果不只是单个频率处幅值不准还会让相邻频率之间的互谱值互相污染。这也是很多人第一次用软件算互谱时容易困惑的地方明明两个信号在 50 Hz 有个很强的相关成分为什么 48 Hz、52 Hz 也出现了不小的互谱幅值大概率就是窗长度不够、泄漏严重、再加上没有做平均。所以真正的工程计算几乎不会用一整段数据直接做一次 FFT 然后相乘而是用分段平均的方法这在第 3 节会详细展开。2. 互谱到底在说什么幅值、相位与相干函数2.1 幅值谱两个信号在哪个频率上“一起变化”(|S_{xy}(f)|) 刻画的是两个信号在频率 (f) 附近共同变化的强度。可以想象两个人合奏如果某个频率上两个人动作幅度都很大并且高度同步互谱幅值就大如果一个人在某频率上很活跃另一个人完全不参与互谱幅值就趋近于零。这与自谱不同——自谱永远只反映单个通道自身的能量无法告诉你“另一个测点是不是也在同一频率上一起动”。有一个容易踩的误区是把互谱幅值和两个自谱幅值的关系理解成简单的几何平均。如果两个信号完全线性相关并且没有噪声确实会成立 (|S_{xy}(f)| \sqrt{S_{xx}(f) S_{yy}(f)})。一旦中间有噪声、非线性或者第三个干扰源互谱幅值会比这个几何平均值小小多少正好反映了“共同部分”所占比例。所以很多人看互谱幅值图时觉得“数值好小”大概率就是两个信号的相关性没那么强而不是计算错了。另外对于实信号互谱具有共轭对称性负频率部分的幅值等于正频率部分但相位互为相反数。工程上通常只分析正频率使用单边谱时还要把正频率部分的幅值乘 2直流分量除外。关于单边谱还是双边谱后面计算部分会再提。2.2 相位谱谁超前谁滞后滞后多少互谱的相位部分是最有价值、也最容易出错的地方。按照我采用的定义 (S_{xy}E[XY^*])它的相位等于 (X(f)) 的相位减去 (Y(f)) 的相位[ \angle S_{xy}(f) \arg X(f) - \arg Y(f) ]如果某个频率处 (\angle S_{xy}(f)) 是正的说明在频率 (f) 上 (x) 超前、(y) 滞后。如果相位是负的则相反。举个具体的例子设 (y(t)x(t-\tau))也就是 (y) 是 (x) 的延迟版本那么理论上 (Y(f)X(f)e^{-j2\pi f\tau})于是[ S_{xy}(f) E\left[X(f) X^*(f) e^{j2\pi f\tau}\right] S_{xx}(f) e^{j2\pi f\tau} ]相位就是 (2\pi f\tau)。只要在某一个频率点测出相位值就能反推时间延迟[ \tau \frac{\angle S_{xy}(f)}{2\pi f} ]这在机械故障诊断、声源定位、超声测速里面很常用。但要注意互谱相位的取值范围是 ([-\pi, \pi])当真实延迟导致的相位移超过 (\pi) 时会发生“相位裹卷”看起来像负的相位直接套上面的公式会算出错误结果。解决办法一是选择较低的频率来换算二是沿频率方向做相位解包裹再用多个频率点做线性拟合。很多新手在 120 Hz 上看到相位变成负值就以为方向反了其实是延迟超过了该频率下 (\pi/2\pi f1/(2f)) 的无模糊范围。2.3 相干函数互谱最重要的“副产品”有了互谱和两个自谱可以顺便计算一个极其重要的指标——幅值相干函数[ \gamma_{xy}^2(f) \frac{|S_{xy}(f)|^2}{S_{xx}(f) S_{yy}(f)} ]它本质上就是频域上的相关系数平方取值范围是 0 到 1。两个信号在某个频率上完全线性相关且没有噪声干扰时相干等于 1完全无关时相干等于 0中间值说明部分相关或者存在非线性、噪声、多处激励叠加等影响。我强烈建议在实际分析中把相干函数和互谱相位画在一起。原因很简单相位谱在没有意义的地方会显示为随机跳动比如两个信号在 2 kHz 处根本不相关互谱相位可能是任意值你去看它会觉得信号关系和实际情况矛盾。但一旦把相干函数放在旁边看到相干在那个频段接近 0就知道这个相位不值得解读。一般经验是只对相干大于 0.5 或 0.6 的频段讨论相位和幅值才有意义。3. 计算方法详解从单次 FFT 到 Welch 平均3.1 为什么不能直接拿一整段 FFT 的结果当互谱有些人图省事直接把两个信号各做一次 FFT然后 (X(k) \cdot Y^(k)) 当作互谱。这在分析确定性正弦信号时还能看但换成随机信号就会出大问题单次 FFT 计算得到的 (XY^) 只是互谱的一个样本它的随机方差和信号本身的能量是同一个量级画出来满屏毛刺完全没法用。这就像只拍一张照片就想判断一个人的平均身高拍到的瞬间姿势、光线、角度都会让结果产生巨大偏差必须有足够的样本做平均。对随机信号的谱估计工程上几乎都用 Welch 平均法。它的核心思想是把长数据切成多段每段做 FFT计算互谱样本再对所有段取平均。这样随机误差会随平均段数的增加而下降互谱的幅值也会更接近真实值。3.2 Welch 平均互谱的完整步骤与归一化具体步骤如下将同步采集到的 (x(t)) 和 (y(t)) 各分成 (K) 段每段长度 (L)相邻段之间允许重叠一定比例。对每一段去掉均值detrend乘上选定的窗函数 (w[n])比如汉宁窗。对每一段分别做 FFT得到该段的频谱 (X_i(f)) 和 (Y_i(f))。计算该段的互谱样本 (P_i(f)X_i(f)Y_i^*(f))。对所有段的互谱样本求平均再按窗能量和采样率归一化得到互功率谱密度估计[ \hat{S}{xy}(f) \frac{1}{K} \sum{i1}^{K} X_i(f) Y_i^*(f) \times \frac{1}{f_s \sum_{n} w^2[n]} ]如果使用单边谱还需要把正频率部分乘 2直流和奈奎斯特频率除外。这里乘的归一化因子是为了保证谱密度单位的正确性也就是让自谱积分后等于信号的总功率。不同软件库可能有细微差别但思路一致。在实际代码里我建议直接用 scipy.signal.csd 函数它会自动处理分段、加窗、重叠、去趋势、归一化。但你必须清楚它默认做的是 (X Y^*)即返回的是 (S_{xy})以第一个输入为参考、单边谱、density 归一化。使用前最好查一下文档或者做一个已知延迟的仿真信号验证方向。3.3 参数选择窗长、重叠率、平均次数怎么搭配Welch 平均的三大参数是窗长nperseg、重叠率noverlap和平均段数。它们的矛盾非常经典窗长 (L) 决定频率分辨率(\Delta f f_s / L)。窗越长谱线越细越能分辨相邻的两个峰。平均段数 (K) 决定随机误差和谱的平滑程度段数越多方差越小谱线越平滑。重叠率可以在不缩短窗长的前提下增加段数但重叠太多相邻段之间的相关性上升额外效益递减。汉宁窗下经验上重叠 50% 是最常见的折中想进一步压低方差可以用 75% 重叠。表里给了几个典型组合假设采样率 (f_s1000) Hz、数据时长 60 秒npersegnoverlap频率分辨率平均段数适用场景10247680.98 Hz约 232 段快速看整体方差小204810240.49 Hz约 57 段均衡推荐204815360.49 Hz约 189 段想要低方差又保持分辨率409620480.24 Hz约 28 段追求频率分辨率段数偏少我的经验是段数尽量不要少于 10最好在 30 段以上。如果为了分辨两个相隔很近的谱峰而被迫使用很长的窗导致平均段数只有个位数那结果基本不可信不如优先保证平均段数、先看清楚整体相干再针对感兴趣频段单独加长窗重新分析。注意Welch 平均后互谱的随机误差不仅来自幅值也来自相位。段数少的时候相位在非相关频段会剧烈跳动这属于正常现象不代表信号真的存在随机相位差更不代表计算错误。4. 手把手实操用 Python 算互谱和相位延迟4.1 构造模拟数据一个延迟相关的双通道信号为了演示完整流程我生成一个 60 秒的模拟信号两个通道共享 50 Hz 和 120 Hz 的正弦分量但第二通道比第一通道整体延迟 5 ms再分别叠加独立噪声。这个场景很像实际测量中同一个振动源传到两个测点的情况。import numpy as np from scipy import signal fs 1000 # 采样率 1 kHz T 60.0 # 60 秒足够做多次平均 t np.arange(int(T * fs)) / fs tau 0.005 # 5 ms 时间延迟 # 源信号50Hz 120Hz 正弦叠加 s 1.0 * np.sin(2 * np.pi * 50 * t) 0.5 * np.sin(2 * np.pi * 120 * t) s_delay np.interp(t - tau, t, s) # 线性插值实现延迟 rng np.random.default_rng(42) x s 0.3 * rng.standard_normal(len(t)) y s_delay 0.35 * rng.standard_normal(len(t))这里我特意把噪声设置成两通道独立因为互谱最大的优点就是与参考信号不相关的噪声在平均后会趋于 0对自谱无能为力的干扰互谱能抑制掉一部分。4.2 用 scipy.signal.csd 计算互谱和相干函数然后调用 scipy.signal.csd 得到互谱同时计算两个自谱再算相干函数f, Sxy signal.csd( x, y, fsfs, windowhann, nperseg2048, noverlap1536, detrendconstant, scalingdensity ) f, Sxx signal.csd(x, x, fsfs, windowhann, nperseg2048, noverlap1536) f, Syy signal.csd(y, y, fsfs, windowhann, nperseg2048, noverlap1536) coh np.abs(Sxy) ** 2 / (Sxx * Syy) coh np.clip(coh, 0.0, 1.0) # 防止数值误差导致轻微越界 mag np.abs(Sxy) phase np.unwrap(np.angle(Sxy))如果你只是想要相干函数也可以直接用 scipy.signal.coherence 函数一步得到。不过自己用互谱和自谱推导一次能帮助你确认对定义的理解有没有偏差。运行完之后幅值谱应该能看到 50 Hz 和 120 Hz 处两个明显的峰相干函数在这两个频率处接近 1其他频率处很低。4.3 从相位谱估计时间延迟并核对结果现在从相位谱里挑一个相干高、幅值大的频率点比如 50 Hz用下面这段代码估计延迟idx50 np.argmin(np.abs(f - 50)) phi50 phase[idx50] delay_est phi50 / (2 * np.pi * f[idx50]) print(f50 Hz 相位: {phi50:.3f} rad) print(f估计延迟: {delay_est * 1e3:.3f} ms真值 5.000 ms)按理论计算5 ms 延迟在 50 Hz 处产生的相位是 (2\pi \times 50 \times 0.005 1.571) rad正好是 90 度处于无模糊范围内所以估计结果应该非常接近 5 ms。如果你在 120 Hz 处也试着算一下会发现 (2\pi \times 120 \times 0.005 3.770) rad已经超过 (\pi)直接把 raw phase 拿来换算会得到错误延迟。这就是前面说的相位裹卷陷阱实际项目中一定要先确认频率点是否满足 (\tau 1/(2f)) 的约束。完整绘图代码我就不贴了但建议至少画出三张图互谱幅值、互谱相位、相干函数三张图共享横轴频率。注意幅值轴用对数坐标否则动态范围太大50 Hz 的峰会把 120 Hz 的小峰压得看不见。5. 实测中那些坑问题现象与排查思路5.1 问题速查表现象 / 可能原因 / 处理思路实测多年我总结过一张互谱分析的问题排查表比较实用现象可能原因处理思路相位谱一团糟毫无规律该频段相干太低信号本质不相关先看相干图只分析相干高的频段相干函数整体偏低噪声太大、存在多源激励或非线性检查测试环境增加平均次数排除干扰源谱峰附近幅值被低估且展宽频谱泄漏增加窗长或换用更强旁瓣抑制的窗相位符号和理论方向相反互谱定义或参数顺序搞反用已知延迟的标定信号验证方向结果数量级差好几倍单边谱/双边谱、density/spectrum 归一化不一致查库文档确认返回值的单位和规范直流或趋势项污染低频未去趋势分段前 detrend必要时高通滤波但注意相位畸变两个通道明明同步采样相位却随频率线性漂移通道间滤波器相移不一致做背对背校准记录并补偿通道相位差相位在某个频率点跳变到负值相位裹卷降低参考频率或用 unwrap 处理5.2 几个容易忽略的细节第一个需要注意的细节是互谱方向定义直接影响相位正负。scipy 的 csd(x, y) 计算的是 (E[XY^*])也就是以 (x) 为基准看 (y) 相对 (x) 的相位。如果你换用 csd(y, x)相位会完全反号。不同软件、不同库的定义可能不同我在处理客户数据时发现过好几例因为没搞清楚方向而把“滞后”分析成“超前”的情况。唯一的处理原则是拿到一个新工具先用一个已知延迟的仿真信号或者给传感器一个已知的机械冲击来标定相位方向。第二个细节是通道校准问题。实际采集系统里两个通道的传感器灵敏度、滤波器的相频特性很难做到完全一致。如果两个通道之间本来就有 5 度的相位失配在做互谱相位分析时就会直接叠加到结果里而且这个误差随频率变化。测量前尽量做一次背对背校准也就是把两个传感器放在同一个振源上同时测量用得到的互谱相位差作为通道修正量再从后续测量的相位里扣除。第三个细节是不要把互谱幅值直接当自谱幅值去解释。互谱幅值的大小取决于两个信号公共成分的强度公共成分占比很低时幅值会远小于任何一个通道的自谱幅值。你有时候看到两个信号各自都有很大的自谱峰但互谱幅值很小这是正常的说明两个峰可能来自彼此独立的振源而不是同一来源。6. 几点实测习惯与经验心得做了这么多年振动信号分析我的习惯是每分析一段数据永远先把相干函数放在最上面看一遍。相干高的频段幅值和相位都值得深入解读相干低的频段再漂亮的相位曲线也只能当噪声处理。这个习惯帮我少走了很多弯路也避免了不少误判。第二点经验是模拟数据是验证计算方法的好帮手。每次在新的软件、新的库、甚至新换一套采集系统时我都会先生成一组包含已知延迟和已知幅值的仿真信号把整个计算流程跑通确认幅值数量级正确、相位方向正确、延迟反推结果在误差范围内。有时候你以为代码没问题实际上一句csd(y, x)写反了就会把自己的分析结论完全带偏。第三点是用互谱做时间延迟估计时尽量选多个频率点做联合拟合而不是只看单个频率。单个频率点受噪声和裹卷影响大多个频率点用线性拟合得到的斜率会稳定得多。如果信号是宽带随机信号还可以考虑用广义互相关方法但互谱相位法依然是理解整个思路的最佳入口。很多时候工程问题卡住的不是数学而是对定义细节和测量链路的理解不够扎实。把 CPSD 真正搞透了再去学传递函数、模态分析、传递路径分析都会顺手很多。