企业官网建设流程全解析

做信号处理的兄弟应该都有这种体会：单看两个传感器的自功率谱，各自频谱都挺干净，可一旦要判断它们之间谁先谁后、谁影响谁、在某几个频率点上的相位差到底是多少，光靠自谱就完全抓瞎了。这时候就得请出互功率谱密度——也就是大家常说的 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, fs=1024, nperseg=1024, noverlap=None, window='hann'): """ 计算 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, dtype=complex) 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, d=1.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 实操中的坑与对策

表格里的问题我基本都踩过。其中最典型的是第一次用互谱算频响函数的时候，发现 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 的整套思路基本讲完了。我个人在实际操作中最深的体会是：刚上手的时候总觉得这是一堆公式和代码的堆叠，做多了才发现真正考验人的永远是“这个频段的互谱为什么长这样”这种看似简单却需要综合判断的问题。下次再遇到相位谱抖动或者相干函数不合格，先别急着调代码，回到传感器布置、采样参数和数据质量上找原因，往往比在算法层面死磕更有效。