# 信号与系统中的相关、互相关、卷积、相干性与谱分析

## 目录

1. [一、概念总览与信号基本操作](#一概念总览与信号基本操作)
2. [二、卷积：定义、计算与 LTI 系统](#二卷积定义计算与-lti-系统)
3. [三、相关与自相关：定义、性质与计算](#三相关与自相关定义性质与计算)
4. [四、卷积、相关与 Fourier/DFT 对偶关系](#四卷积相关与-fourierdft-对偶关系)
5. [五、能量谱、功率谱与互谱](#五能量谱功率谱与互谱)
6. [六、相关性与统计信号分析](#六相关性与统计信号分析)
7. [七、相干性及其扩展](#七相干性及其扩展)
8. [八、谱与相干性的估计方法](#八谱与相干性的估计方法)
9. [九、通信、检测与系统辨识应用](#九通信检测与系统辨识应用)
10. [十、典型例题、常见误区与检查清单](#十典型例题常见误区与检查清单)

—

## 一、概念总览与信号基本操作

在信号处理中，经常会看到类似的乘加表达式：

$$
\sum_n x[n]h[n-k],
\qquad
\sum_n x[n]y[n+k],
\qquad
\sum_n x[n]y^*[n-k].
$$

它们的形式很接近，但含义并不相同。最容易混淆的几个运算是：

| 运算 | 典型表达式 | 主要用途 |
|—|—|—|
| 线性卷积 | $y(t)=\int x(\tau)h(t-\tau)\,d\tau$ | 求 LTI 系统输出、描述滤波 |
| 离散线性卷积 | $y[n]=\sum_k x[k]h[n-k]$ | 数字滤波、系统响应 |
| 自相关 | $R_{xx}(\tau)=\int x(t)x^*(t-\tau)\,dt$ | 测量信号与自身的相似性 |
| 互相关 | $R_{xy}(\tau)=\int x(t)y^*(t-\tau)\,dt$ | 测量两个信号的相似性和相对延迟 |
| 内积 | $\langle x,y\rangle=\int x(t)y^*(t)\,dt$ | 某一个固定对齐位置的相似性 |

可以先记住以下核心区别：

> **卷积包含“反转后平移”，相关包含“共轭、反转后平移”或等价的相关排列。**

对于实值信号，共轭不起作用，所以卷积和相关在图形操作上看起来更加相似；对于复值信号，共轭是相关定义的重要组成部分。

另一个关键区别是：

– 卷积描述一个信号经过系统冲激响应后的累积作用；
– 相关描述两个信号在不同相对位移下的匹配程度。

从 Fourier/DFT 对偶关系看，还应特别区分“相乘”的具体对象：

| 时域目标 | 频域对应关系 | 结果类型 |
|—|—|—|
| 点对点乘法 $z[n]=x[n]y[n]$ | $Z[k]=\frac1N(X\circledast_NY)[k]$ | 频域循环卷积 |
| 线性卷积 $x*h$ | $XH$ | 频域普通逐点乘法 |
| 循环卷积 $x\circledast_Nh$ | $X[k]H[k]$ | $N$ 点频域普通乘法 |
| 互相关 $R_{xy}[m]$ | $X[k]Y^*[k]$ 后 IDFT | 按 lag 展开的相关序列 |
| 零延迟内积 $\sum_nx[n]y^*[n]$ | $\frac1N\sum_kX[k]Y^*[k]$ | 一个标量 |
| 自功率谱 | $XX^*=|X|^2$ | 非负实数 |
| Hermitian 梯度 $\boldsymbol X^H\boldsymbol e$ | $X^*E$ 后 IDFT | 自适应滤波抽头更新方向 |

因此，“时域信号相乘可用频域复共轭相乘替代”不是普遍规律。时域点乘对应频域卷积；频域一方取共轭后逐点相乘，通常对应相关、复内积、互谱或 Hermitian 梯度。看到 $XE^*$ 或 $X^*E$ 时，必须结合相关方向、信号角色和后续是求和还是 IFFT 判断其含义。

—
### 预备知识：信号的平移、反转与共轭

理解卷积和相关之前，必须熟悉三个基本操作。

#### 时间平移

对信号 $x(t)$：

$$
x(t-t_0)
$$

表示向右延迟 $t_0$；

$$
x(t+t_0)
$$

表示向左移动 $t_0$。

离散信号中：

$$
x[n-n_0]
$$

表示延迟 $n_0$ 个采样点。

#### 时间反转

时间反转为

$$
x(-t)
$$

离散形式为

$$
x[-n].
$$

几何上，这相当于把信号关于 $t=0$ 或 $n=0$ 镜像翻转。

#### 复共轭

复信号

$$
x(t)=x_R(t)+jx_I(t)
$$

的共轭为

$$
x^*(t)=x_R(t)-jx_I(t).
$$

复内积通常定义为

$$
\langle x,y\rangle=\int x(t)y^*(t)\,dt
$$

或离散形式

$$
\langle x,y\rangle=\sum_n x[n]y^*[n].
$$

这样定义可以保证

$$
\langle x,x\rangle=\int |x(t)|^2dt\ge 0,
$$

因此它能够表示能量或平方范数。

#### 三个操作的组合

卷积中的被卷积信号通常经历：

1. 反转：$h(\tau)\to h(-\tau)$；
2. 平移：$h(-\tau)\to h(t-\tau)$；
3. 与 $x(\tau)$ 相乘并积分。

相关中的信号通常经历：

1. 共轭：$y(\tau)\to y^*(\tau)$；
2. 反转和平移：$y^*(\tau)\to y^*(\tau-t)$；
3. 与 $x(\tau)$ 相乘并积分。

—
### 能量信号与功率信号

讨论能量谱和功率谱前，必须先区分能量信号与功率信号。这个分类决定了应使用 Fourier 变换还是功率谱密度来描述频率分布。

#### 连续时间信号的能量和平均功率

对信号 $x(t)$，总能量定义为

$$
\boxed{
E_x=\int_{-\infty}^{\infty}|x(t)|^2dt
}.
$$

平均功率定义为长期时间平均：

$$
\boxed{
P_x=\lim_{T\to\infty}\frac{1}{2T}
\int_{-T}^{T}|x(t)|^2dt
}.
$$

#### 离散时间信号的能量和平均功率

离散时间信号的能量为

$$
\boxed{
E_x=\sum_{n=-\infty}^{\infty}|x[n]|^2
}.
$$

平均功率为

$$
\boxed{
P_x=\lim_{N\to\infty}\frac{1}{2N+1}
\sum_{n=-N}^{N}|x[n]|^2
}.
$$

若时间起点不重要，也可以使用 $0$ 到 $N-1$ 的长窗口定义相同的极限。

#### 能量信号

若

$$
00$、$b>0$ 且 $a\ne b$。对于 $t\ge0$：

$$
\begin{aligned}
y(t)
&=\int_0^t e^{-a\tau}e^{-b(t-\tau)}d\tau\\
&=e^{-bt}\int_0^t e^{-(a-b)\tau}d\tau\\
&=\frac{e^{-bt}-e^{-at}}{a-b}u(t).
\end{aligned}
$$

卷积结果的起始时间由两个信号的起始时间共同决定，系统响应的衰减形式由输入和系统极点共同决定。

#### 连续时间卷积的重叠解释

对有限时长信号，互相关

$$
R_{xy}(\tau)=\int x(t)y^*(t-\tau)dt
$$

同样可以理解为两个信号在相对平移后重叠部分的加权面积。两个波形越匹配，重叠乘积的积分越大。

—
### 卷积的重要性质与常用结论

#### 卷积的交换律

$$
x*h=h*x.
$$

这说明卷积中两个信号的角色可以交换。

#### 卷积的结合律

$$
(x*h)*g=x*(h*g).
$$

在串联系统中，可以先合并系统的冲激响应，再与输入卷积。

#### 卷积的分配律

$$
x*(h_1+h_2)=x*h_1+x*h_2.
$$

这对应线性系统的叠加性质。

#### 微分与卷积

在适当边界条件下：

$$
\frac{d}{dt}(x*h)
=\frac{dx}{dt}*h
=x*\frac{dh}{dt}.
$$

这使得可以先微分输入或系统响应，再进行卷积。

#### 因果性与卷积支撑集

若 $x(t)$ 和 $h(t)$ 都是因果信号，即 $t<0$ 时为零，则 $$ (x*h)(t)=0,\qquad t<0. $$ 更一般地，两个信号的非零区间卷积后，其支撑集是原两个支撑集的 Minkowski 和。有限离散序列长度满足 $N+M-1$，就是这一结论的离散体现。 #### 复数信号的卷积：不使用共轭乘法 当信号为复值时，线性卷积的定义不变： $$ \boxed{ y(t)=\int_{-\infty}^{\infty}x(\tau)\,h(t-\tau)\,d\tau, \qquad y[n]=\sum_{k=-\infty}^{\infty}x[k]\,h[n-k] }. $$ 被积函数中的 $x(\tau)\,h(t-\tau)$（或 $x[k]\,h[n-k]$）是普通的复数乘法，不含任何共轭。这一点在信号从实值推广到复值时保持不变。 共轭乘法出现在另一组运算中——互相关、自相关和复内积： | 运算 | 时域乘法形式 | 频域对应 | |---|---|---| | 线性卷积 $(x*h)(t)$ | $x(\tau)\,h(t-\tau)$（普通乘法） | $X(\omega)H(\omega)$ | | 互相关 $R_{xy}(\tau)$ | $x(t)\,y^*(t-\tau)$（共轭乘法） | $X(\omega)Y^*(\omega)$ | | 复内积 $\langle x,y\rangle$ | $x(t)\,y^*(t)$（共轭乘法） | $\frac{1}{2\pi}\int X(\omega)Y^*(\omega)d\omega$ | 频域同样如此：卷积定理给出 $Y=XH$，不含共轭；相关定理给出 $XY^*$，含共轭。 **复指数验证。** 取 $x[n]=e^{j\omega_0 n}$、$h[n]=e^{j\omega_1 n}$（均为复信号）。卷积乘积为 $$ x[k]\,h[n-k]=e^{j\omega_0 k}\cdot e^{j\omega_1(n-k)}=e^{j\omega_1 n}\cdot e^{j(\omega_0-\omega_1)k}. $$ 对 $k$ 求和后，只有 $\omega_0=\omega_1$ 时不抵消，输出频率为 $\omega_0+\omega_1$ 方向的贡献——频率相加。 而互相关乘积为 $$ x[n]\,h^*[n-m]=e^{j\omega_0 n}\cdot e^{-j\omega_1(n-m)}=e^{j(\omega_0-\omega_1)n}\cdot e^{j\omega_1 m}. $$ $\omega_0=\omega_1$ 时乘积不再旋转，相关峰最大——共轭把频率比较变成了减法。 因此，“复数信号需要共轭乘法”是相关的性质，不是卷积因信号变为复数后的自动规则。关于共轭在相关中的作用和必要性，详见第四章“相关是带共轭的卷积变体”和“为什么相关中有共轭”。 --- ### 系统分析中的卷积 #### LTI 系统的冲激响应 线性时不变系统（LTI）的核心性质是：系统完全由冲激响应 $h(t)$ 或 $h[n]$ 描述。 连续时间系统输入为 $x(t)$ 时： $$ \boxed{ y(t)=x(t)*h(t) }. $$ 离散时间系统输入为 $x[n]$ 时： $$ \boxed{ y[n]=x[n]*h[n] }. $$ #### 从冲激分解推导卷积 利用连续时间冲激分解： $$ x(t)=\int_{-\infty}^{\infty}x(\tau)\delta(t-\tau)d\tau. $$ 由于系统对位于 $\tau$ 的冲激响应为 $h(t-\tau)$，根据线性叠加： $$ y(t)=\int_{-\infty}^{\infty}x(\tau)h(t-\tau)d\tau. $$ 离散时间中： $$ x[n]=\sum_kx[k]\delta[n-k], $$ 所以 $$ y[n]=\sum_kx[k]h[n-k]. $$ #### 冲激分解的严格步骤 上述“分解—响应—叠加”的思想可以拆成四个可以逐条检验的步骤，方便理解每一步用到的假设： 1. **采样表示**：把输入写成加权冲激的叠加 $$ x[n]=\sum_{k\in\mathbb{Z}}x[k]\,\delta[n-k]. $$ 这是恒等式，不涉及系统性质。 2. **时不变性**：设系统对 $\delta[n]$ 的响应为 $h[n]$，则对 $\delta[n-k]$ 的响应为 $h[n-k]$。 3. **齐次性**：对系数 $x[k]$ 的输入 $x[k]\delta[n-k]$，系统输出为 $x[k]h[n-k]$。 4. **可加性**：把所有 $k$ 上的响应相加，得到 $$ y[n]=\sum_{k\in\mathbb{Z}}x[k]h[n-k]=(x*h)[n]. $$ 其中步骤 2 用到时不变性，步骤 3、4 合起来对应线性性质。缺失其中任何一条，卷积公式都不成立。连续时间的推导只需把求和换成积分，把 $\delta[n-k]$ 换成 $\delta(t-\tau)$，并保证在合适的函数空间中允许交换积分与线性算子的顺序。 #### 有限序列卷积的 Toeplitz 矩阵形式 对有限长度序列 $x=[x_0,\dots,x_{N-1}]^{\mathsf T}$、$h=[h_0,\dots,h_{M-1}]^{\mathsf T}$，线性卷积 $y=x*h$ 长度为 $N+M-1$，可以写成矩阵–向量乘积： $$ y=H\,x, $$ 其中 $H$ 是由 $h$ 构造的 $(N+M-1)\times N$ **Toeplitz 矩阵**，第 $k$ 列是把 $h$ 向下平移 $k$ 位、其余元素补零得到的向量。以 $N=3$、$M=3$ 为例： $$ H= \begin{bmatrix} h_0 & 0 & 0 \\ h_1 & h_0 & 0 \\ h_2 & h_1 & h_0 \\ 0 & h_2 & h_1 \\ 0 & 0 & h_2 \end{bmatrix}, \qquad y=H\,[x_0,x_1,x_2]^{\mathsf T}. $$ 由交换律，也可以写成 $y=X\,h$，其中 $X$ 是由 $x$ 构造的 $(N+M-1)\times M$ Toeplitz 矩阵。这个观察带来几个直接的用处： - **数值实现**：$H$ 每一条对角线的值相同，无需显式存储整个矩阵；MATLAB 中的 `convmtx`、Python 中的 `scipy.linalg.toeplitz` 就是根据这一结构构造的。 - **反卷积**：估计 $x$ 时可以求解 $\min_x\|Hx-y\|^2$ 或 Tikhonov 正则化问题，比时域递推更稳。 - **循环卷积的类比**：把 Toeplitz 结构“回绕”成 $N\times N$ 的 **循环矩阵**（circulant），就是 $N$ 点循环卷积；循环矩阵总是由 DFT 对角化，这正是第四章 FFT 快速卷积的代数根源。 #### 因果性与卷积支撑集的更精细描述 对连续时间信号，定义支撑集 $\operatorname{supp}(x)=\overline{\{t:x(t)\ne 0\}}$；离散情形类似。则线性卷积满足 $$ \operatorname{supp}(x*h)\subseteq \operatorname{supp}(x)+\operatorname{supp}(h), $$ 其中右侧是 Minkowski 和 $\{a+b:a\in A,b\in B\}$。由此推出： - 若 $x$ 与 $h$ 均因果（支撑集在 $[0,\infty)$ 上），则 $x*h$ 也因果； - 若 $\operatorname{supp}(x)\subseteq[a_1,b_1]$，$\operatorname{supp}(h)\subseteq[a_2,b_2]$，则 $\operatorname{supp}(x*h)\subseteq[a_1+a_2,b_1+b_2]$； - 反因果、双边信号也可以按同样方式判定其卷积输出的支撑区间。 因果 LTI 系统的等价刻画是 $h(t)=0,\ t<0$。对可用输入 $x(t)$ 的响应 $$ y(t)=\int_{0}^{\infty}h(\tau)x(t-\tau)d\tau $$ 只使用过去和现在的输入。工程上判断系统是否因果，最快的办法就是画出 $h(t)$ 或 $h[n]$，看是否含有 $t<0$ 或 $n<0$ 的样本。 #### 频域系统函数 对 LTI 系统做 Fourier 变换： $$ Y(\omega)=X(\omega)H(\omega). $$ 因此 $$ H(\omega)=\frac{Y(\omega)}{X(\omega)} $$ 表示系统对各频率分量的幅度和相位响应。 若输入为复指数 $$ x(t)=e^{j\omega t}, $$ 则输出为 $$ y(t)=H(\omega)e^{j\omega t}. $$ 这说明复指数是 LTI 系统的特征函数，卷积在频域中变成乘法。 #### 级联和并联系统 两个 LTI 系统级联时： $$ h_{\mathrm{eq}}=h_1*h_2, $$ 频域为 $$ H_{\mathrm{eq}}(\omega)=H_1(\omega)H_2(\omega). $$ 两个系统并联时： $$ h_{\mathrm{eq}}=h_1+h_2, $$ 频域为 $$ H_{\mathrm{eq}}(\omega)=H_1(\omega)+H_2(\omega). $$ #### 卷积与系统稳定性 离散时间 BIBO 稳定的 LTI 系统需要满足 $$ \sum_n|h[n]|<\infty. $$ 连续时间对应条件为 $$ \int_{-\infty}^{\infty}|h(t)|dt<\infty. $$ 这是因为对于有界输入 $|x(t)|\le M$： $$ |y(t)| \le M\int|h(\tau)|d\tau. $$ #### BIBO 稳定性的充要性与常见判据 上式给出的是充分性。必要性可以通过构造最坏输入证明：取 $$ x[n]=\operatorname{sgn}(h^*[-n]),\qquad |x[n]|\le 1, $$ 则 $n=0$ 处的输出 $$ y[0]=\sum_k x[k]h[-k]=\sum_k |h[-k]|=\sum_k|h[k]|. $$ 若 $\sum|h[n]|=\infty$，就存在有界输入使输出无界；因此 $h\in\ell^1$ 是离散时间 BIBO 稳定的充要条件。连续时间的证明思路一致，只需保证选取的 $x(t)$ 可测且有界。 判断稳定性时，常用几种角度对照： - **时域**：直接计算或估计 $\sum|h[n]|$、$\int|h(t)|dt$，例如判断 $h[n]=a^nu[n]$ 是否稳定，只需 $|a|<1$； - **传递函数**：连续系统的所有极点位于左半平面（严格 $\operatorname{Re}\{s_i\}<0$）等价于因果 LTI 系统 BIBO 稳定；离散系统对应所有极点位于单位圆内（$|z_i|<1$）； - **频率响应存在性**：$h\in\ell^1$ 保证 DTFT $H(e^{j\hat\omega})$ 作为普通连续函数存在，工程上常反过来用“$H$ 平滑无奇点”作为快速筛查依据； - **零输入响应**：LTI 系统的零输入响应指数衰减到 0，是稳定性的等价体现。 因果性与稳定性是彼此独立的两个性质：可以有稳定但非因果（如理想低通的对称 sinc）、也可以有因果但不稳定（如 $h[n]=2^nu[n]$）的系统；只有两者同时满足，才能实时且可靠地实现。 --- --- ## 三、相关与自相关：定义、性质与计算 ### 相关的基本思想 相关是一种相似性度量。将一个信号相对另一个信号移动，在每个相对位移处计算乘积并累加： - 若两者在当前位移下形状相似，乘积大多同号或相位接近，累加值较大； - 若两者不匹配，乘积会发生正负抵消或相位抵消，累加值较小。 因此，相关函数的峰值通常表示最佳匹配位置。 ### 连续时间互相关 一种常用定义是 $$ \boxed{ R_{xy}(\tau)=\int_{-\infty}^{\infty}x(t)y^*(t-\tau)\,dt }. $$ 这里 $x(t)$ 是参考信号，$y(t)$ 被共轭、反转并平移。 也有文献使用 $$ R_{xy}(\tau)=\int_{-\infty}^{\infty}x^*(t)y(t+\tau)\,dt. $$ 这两个定义在改变相关方向或取共轭后可以相互转换。使用公式前，必须明确所采用的定义。 ### 离散时间互相关 对应的离散时间定义为 $$ \boxed{ R_{xy}[m]=\sum_{n=-\infty}^{\infty}x[n]y^*[n-m] }. $$ 另一种常见写法是 $$ R_{xy}[m]=\sum_n x^*[n]y[n+m]. $$ 只要前后一致，二者都可以使用；但相关峰的正负延迟方向可能不同。 ### 互相关的几何意义 把信号看作向量。对每个延迟 $m$，定义平移后的信号 $$ y_m[n]=y[n-m]. $$ 则 $$ R_{xy}[m]=\langle x,y_m\rangle. $$ 因此，互相关实际上是在扫描不同平移量下的向量内积。内积幅度越大，说明两个信号越相似。 对于复信号，$R_{xy}[m]$ 一般是复数： - 幅度表示匹配程度； - 相位表示相对相位信息； - 峰值位置表示相对时延。 ### 归一化互相关 不同信号的幅度可能不同，仅比较未归一化相关值会受到能量影响。常见的归一化互相关为 $$ \rho_{xy}[m] =\frac{R_{xy}[m]} {\sqrt{R_{xx}[0]R_{yy}[0]}}. $$ 其中 $$ R_{xx}[0]=\sum_n|x[n]|^2, \qquad R_{yy}[0]=\sum_n|y[n]|^2. $$ 由 Cauchy–Schwarz 不等式，有 $$ |\rho_{xy}[m]|\le 1. $$ 归一化后可以更公平地比较不同能量信号的匹配程度。 ### 互相关与内积的关系 内积是零延迟下的相关： $$ \langle x,y\rangle=R_{xy}[0] $$ 或在另一种约定下等于 $R_{yx}[0]$。相关则把一个标量内积扩展为关于延迟的函数。 --- ### 自相关 #### 定义 信号与自身的互相关称为自相关。 连续时间： $$ \boxed{ R_{xx}(\tau)=\int_{-\infty}^{\infty}x(t)x^*(t-\tau)\,dt }. $$ 离散时间： $$ \boxed{ R_{xx}[m]=\sum_n x[n]x^*[n-m] }. $$ #### 零延迟自相关等于能量 在能量有限的情况下： $$ R_{xx}(0)=\int|x(t)|^2dt=E_x $$ 或 $$ R_{xx}[0]=\sum_n|x[n]|^2=E_x. $$ 因此，自相关在零延迟处通常具有最大峰值或至少满足最大幅度不超过零延迟值。 #### 自相关的共轭对称性 对于复信号： $$ \boxed{ R_{xx}(-\tau)=R_{xx}^*(\tau) }. $$ 离散形式为 $$ R_{xx}[-m]=R_{xx}^*[m]. $$ 因此： - 实信号的自相关是偶函数； - 复信号的自相关不一定是实数，但具有共轭对称性； - $R_{xx}(0)$ 必定是非负实数。 #### 周期信号的自相关 对于周期信号，通常使用时间平均形式。例如周期为 $T_0$ 的信号： $$ R_{xx}(\tau) =\frac{1}{T_0}\int_{t_0}^{t_0+T_0}x(t)x^*(t-\tau)\,dt. $$ 周期信号的自相关也是周期函数，周期通常与原信号相同。周期结构越明显，自相关中的周期峰越明显。 #### 自相关的应用 自相关常用于： - 检测周期和基频； - 测量信号持续时间； - 估计噪声白化程度； - 估计随机过程的统计结构； - 语音基音检测； - 雷达和声纳回波分析。 --- ### 相关的常用性质 #### 相关的共轭对称关系 互相关满足 $$ \boxed{ R_{yx}(\tau)=R_{xy}^*(-\tau) }. $$ 离散形式为 $$ R_{yx}[m]=R_{xy}^*[-m]. $$ 对于实信号： $$ R_{yx}(\tau)=R_{xy}(-\tau). $$ #### Cauchy–Schwarz 不等式 相关幅度满足 $$ |R_{xy}(\tau)|^2 \le R_{xx}(0)R_{yy}(0). $$ 因此，相关值的幅度不会超过两个信号能量平方根的乘积。 #### 相关的平移性质 若 $$ y(t)=x(t-t_0), $$ 则互相关会发生相应平移。以 $$ R_{xy}(\tau)=\int x(t)y^*(t-\tau)dt $$ 为定义时： $$ R_{xy}(\tau)=R_{xx}(\tau+t_0). $$ 因为 $y(t)=x(t-t_0)$ 时，$y^*(t-\tau)=x^*(t-\tau-t_0)$，所以 $R_{xy}$ 的峰位为 $\tau=-t_0$。若采用 $R_{yx}$ 或另一种相关方向，峰位会反号。结论是： > 一个信号相对另一个信号延迟，会使相关峰移动；峰移动的正负方向取决于相关定义。

### 离散卷积与相关的计算方法

#### 直接计算离散线性卷积

若 $x[n]$ 长度为 $N$，$h[n]$ 长度为 $M$，则线性卷积长度为

$$
N+M-1.
$$

计算公式为

$$
y[n]=\sum_{k=0}^{N-1}x[k]h[n-k],
$$

其中超出 $h$ 有效下标的项取零。

对于有限序列，通常建立下标范围：

$$
0\le n\le N+M-2.
$$

#### 卷积表格法

计算离散卷积时，可以：

1. 写出 $x[k]$；
2. 写出 $h[n-k]$；
3. 对 $h$ 进行反转并随 $n$ 平移；
4. 找出两个序列的重叠区间；
5. 将重叠样本相乘并求和。

#### 直接计算离散互相关

对于

$$
R_{xy}[m]=\sum_n x[n]y^*[n-m],
$$

计算流程是：

1. 对 $y[n]$ 取共轭；
2. 对共轭后的序列做时间反转；
3. 将其平移；
4. 与 $x[n]$ 重叠相乘并求和。

也可以直接根据公式逐个 $m$ 计算，不必进行图形反转。

#### 使用 FFT 加速卷积

直接计算长度为 $N$ 和 $M$ 序列的卷积，复杂度约为 $O(NM)$。利用 FFT 可以将复杂度降低到近似

$$
O(L\log L),
$$

其中 $L$ 是足够大的 FFT 长度。

线性卷积的 FFT 实现：

1. 取 $L\ge N+M-1$；
2. 将 $x$、$h$ 补零到长度 $L$；
3. 计算 $X=\operatorname{FFT}(x)$ 和 $H=\operatorname{FFT}(h)$；
4. 计算 $Y=XH$；
5. 计算 $y=\operatorname{IFFT}(Y)$；
6. 保留前 $N+M-1$ 个有效样本。

互相关的 FFT 实现：

$$
R_{xy,\mathrm{circ}}
=\operatorname{IFFT}\left\{X[k]Y^*[k]\right\}.
$$

实际使用时必须根据相关定义对结果进行循环移位或重新排列。

#### 相关系数与滑动窗口相关

在实际检测中，常常不是对两个完整序列做相关，而是将模板 $s[n]$ 在观测信号 $r[n]$ 上滑动：

$$
C[m]=\sum_n r[n]s^*[n-m].
$$

为避免模板能量和局部信号能量影响，可以采用局部归一化：

$$
C_{\mathrm{norm}}[m]
=\frac{\sum_n r[n]s^*[n-m]}
{\sqrt{\sum_n|r[n]|^2\sum_n|s[n-m]|^2}}.
$$

这类指标常用于模板匹配和同步检测。

#### 复有限序列互相关的完整例题

设两个长度为 3 的复序列

$$
x[n]=\{1,\ j,\ -1\}_{n=0,1,2},
\qquad
y[n]=\{1,\ 1+j,\ 0\}_{n=0,1,2}.
$$

按本文默认定义 $R_{xy}[m]=\sum_n x[n]y^*[n-m]$，$m$ 的有效范围为 $-(N_y-1)\le m\le N_x-1$，即 $-2\le m\le 2$。为便于计算，先写出 $y^*[n]=\{1,\ 1-j,\ 0\}$，然后对每个 $m$ 找出使 $n$ 与 $n-m$ 同时落在有效区间的项：

– $m=-2$：仅 $n=0$ 满足，
$$
R_{xy}[-2]=x[0]y^*[2]=1\cdot 0=0.
$$
– $m=-1$：$n=0,1$，
$$
R_{xy}[-1]=x[0]y^*[1]+x[1]y^*[2]=(1)(1-j)+(j)(0)=1-j.
$$
– $m=0$：$n=0,1,2$，
$$
R_{xy}[0]=x[0]y^*[0]+x[1]y^*[1]+x[2]y^*[2]=1+(j)(1-j)+0=1+(j+1)=2+j.
$$
– $m=1$：$n=1,2$，
$$
R_{xy}[1]=x[1]y^*[0]+x[2]y^*[1]=j+(-1)(1-j)=j-1+j=-1+2j.
$$
– $m=2$：$n=2$，
$$
R_{xy}[2]=x[2]y^*[0]=-1.
$$

于是

$$
R_{xy}[m]=\{0,\ 1-j,\ 2+j,\ -1+2j,\ -1\}_{m=-2,\dots,2}.
$$

取模平方 $|R_{xy}[m]|^2=\{0,2,5,5,1\}$，最大匹配位置在 $m=0$ 与 $m=1$ 之间并列，反映出 $y$ 中的 $1+j$ 与 $x$ 中 $\{1,j\}$ 的两种对齐方式相似度接近。若同时计算 $R_{yx}[m]$，可以验证 $R_{yx}[m]=R_{xy}^*[-m]$，例如 $R_{yx}[-1]=R_{xy}^*[1]=-1-2j$，与逐项直接算出的结果一致。

这个例子说明：在复信号中，$R_{xy}[m]$ 每个分量的相位携带了对齐后两段信号的相位差信息；只看幅度会丢掉相对相位，因而不能替代复数结果本身。

#### lag 方向与实现约定的对照

不同教材、不同数值库对“lag 的正方向”约定不同，这是初学者最容易出错的地方。对本文默认定义 $R_{xy}[m]=\sum_n x[n]y^*[n-m]$：

| 约定 | 峰位含义 | 数值示例 |
|—|—|—|
| $R_{xy}[m]=\sum_n x[n]y^*[n-m]$（本文） | 若 $y[n]=x[n-m_0]$（$y$ 右移、相对 $x$ 延迟 $m_0>0$），峰在 $m=-m_0$；正 $m$ 表示按本文公式扫描到相反方向 | 用单位脉冲或已知移位序列验证，不凭“正 lag”字面判断 |
| $R_{xy}[m]=\sum_n x^*[n]y[n+m]$ | 与上式等价地描述另一方向，峰位符号随索引定义改变 | 先固定公式，再解释正负 |
| NumPy `np.correlate(x, y, mode=’full’)` | 返回 $\sum_n x[n+m]\,\overline{y[n]}$，正 lag 对应 $x$ 相对 $y$ 超前 | 结果索引从 $-(N_y-1)$ 到 $N_x-1$ 需自行换算 |
| MATLAB `xcorr(x, y)` | 返回 $\sum_n x[n+m]\,y^*[n]$，正 lag 对应 $x$ 超前 $y$ | 索引以中心为零 lag |

工程上避免出错的三步做法：

1. 在报告或代码中写清楚“相关的定义式与 lag 的物理含义”；
2. 用一个已知延迟的合成信号（如 $y=x$ 平移 $k$ 个采样）跑一遍，验证峰位符号；
3. 若不同来源的结果符号相反，直接对整段结果做 $m\to -m$ 或额外取共轭即可换算。

#### 归一化估计与偏/无偏估计

对宽平稳过程 $x[n]$，用有限观测 $\{x[0],\dots,x[N-1]\}$ 估计自相关时有两种常见形式：

– **有偏估计**（NumPy/`scipy.signal.correlate` 默认）
$$
\hat R^{(b)}_{xx}[m]=\frac{1}{N}\sum_{n=0}^{N-1-|m|}x[n+|m|]\,x^*[n],\quad |m|\le N-1.
$$
分母固定为 $N$，$|m|$ 越大参与求和的样本越少，估计值系统性偏低，但整个序列作为正定核更稳定，对应的功率谱估计非负。
– **无偏估计**
$$
\hat R^{(u)}_{xx}[m]=\frac{1}{N-|m|}\sum_{n=0}^{N-1-|m|}x[n+|m|]\,x^*[n].
$$
期望值等于真实自相关，但大 $|m|$ 时方差急剧增大，用它做 Fourier 变换得到的功率谱可能出现负值。

工程实践中：

– 若只关心峰位或波形匹配，两种估计差别不大，用有偏估计更安全；
– 若要进一步做 PSD、Wiener–Khinchin 或参数化谱估计，务必使用有偏估计或加窗后的 Blackman–Tukey 方法，以保证正半定性；
– 当 $N$ 有限、且要看到较大 $|m|$ 的自相关时，可以只信任 $|m|\ll N$ 的部分（例如 $|m|\le N/4$），并在报告中标注估计带宽和有效样本数。

对互相关也有完全类似的定义；此时 $N$ 换成两段信号可用的重叠样本数 $N-|m|$，同样存在偏、无偏两种归一化方式。

—

—

## 四、卷积、相关与 Fourier/DFT 对偶关系

### 相关是带共轭的卷积变体

定义

$$
\tilde y(t)=y^*(-t).
$$

则

$$
\begin{aligned}
(x*\tilde y)(\tau)
&=\int x(t)\tilde y(\tau-t)\,dt\\
&=\int x(t)y^*(t-\tau)\,dt\\
&=R_{xy}(\tau).
\end{aligned}
$$

因此

$$
\boxed{
R_{xy}(\tau)=x(\tau)*y^*(-\tau)
}.
$$

离散形式为

$$
\boxed{
R_{xy}[m]=x[m]*y^*[-m]
}.
$$

这里等式中的卷积变量和结果索引需要按照具体定义对齐，但本质关系是：

> **互相关等于一个信号与另一个信号的“共轭时间反转版本”做卷积。**

### 为什么相关中有共轭

对于复指数信号：

$$
x[n]=e^{j\omega_0n},
\qquad
y[n]=e^{j\omega_1n},
$$

逐点共轭相乘得到

$$
x^*[n]y[n]=e^{j(\omega_1-\omega_0)n}.
$$

当 $\omega_1=\omega_0$ 时，结果为常数，累加不会相互抵消；当频率不同时，结果继续旋转，累加会部分抵消。因此，共轭相乘测量的是相对相位和频率差。

若使用普通乘法：

$$
x[n]y[n]=e^{j(\omega_1+\omega_0)n},
$$

得到的是频率相加，而不是比较两者的相对频率。

### 什么时候卷积和相关形式相同

当信号为实值且满足偶对称时：

$$
y(t)=y(-t),
$$

共轭和反转都不会改变信号，此时相关与卷积可能具有相同的表达形式。

但一般情况下，尤其是信号不对称或为复信号时，不能把相关直接当成卷积。

—
### 卷积和相关的 Fourier 变换关系

采用连续时间 Fourier 变换约定：

$$
X(\omega)=\int_{-\infty}^{\infty}x(t)e^{-j\omega t}\,dt,
$$

$$
x(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}X(\omega)e^{j\omega t}\,d\omega.
$$

#### 卷积定理

若

$$
y(t)=x(t)*h(t),
$$

则

$$
\boxed{
Y(\omega)=X(\omega)H(\omega)
}.
$$

反过来：

$$
\boxed{
\mathcal{F}\{x(t)h(t)\}
=\frac{1}{2\pi}X(\omega)*H(\omega)
}.
$$

因此：

– 时域卷积对应频域逐点相乘；
– 时域逐点相乘对应频域卷积。

#### 相关定理

对于定义

$$
R_{xy}(\tau)=\int x(t)y^*(t-\tau)\,dt,
$$

有

$$
\boxed{
\mathcal{F}\{R_{xy}(\tau)\}=X(\omega)Y^*(\omega)
}.
$$

根据相关方向，也可能写为

$$
\mathcal{F}\{R_{yx}(\tau)\}=Y(\omega)X^*(\omega).
$$

这两个表达式不是简单的“一个正确、一个错误”，而是对应不同的相关定义和方向。

#### 相关定理的两种等价推导

**推导一：直接展开定义。** 记 $Z(\omega)=\mathcal{F}\{R_{xy}(\tau)\}$，则

$$
\begin{aligned}
Z(\omega)
&=\int_{-\infty}^{\infty}\left[\int_{-\infty}^{\infty}x(t)y^*(t-\tau)dt\right]e^{-j\omega\tau}d\tau\\
&=\int_{-\infty}^{\infty}x(t)\int_{-\infty}^{\infty}y^*(t-\tau)e^{-j\omega\tau}d\tau\,dt.
\end{aligned}
$$

在内层积分中令 $u=t-\tau$，$\tau=t-u$，$d\tau=-du$，改变积分限：

$$
\int_{-\infty}^{\infty}y^*(u)e^{-j\omega(t-u)}du
=e^{-j\omega t}\int y^*(u)e^{j\omega u}du
=e^{-j\omega t}\,Y^*(\omega).
$$

代回外层积分：

$$
Z(\omega)=Y^*(\omega)\int x(t)e^{-j\omega t}dt=X(\omega)Y^*(\omega).
$$

**推导二：借助“共轭时间反转”与卷积定理。** 记 $\tilde y(t)=y^*(-t)$，则由基本变换公式

$$
\mathcal{F}\{y^*(-t)\}=Y^*(\omega).
$$

再利用第三章得到的 $R_{xy}(\tau)=(x*\tilde y)(\tau)$ 与卷积定理

$$
\mathcal{F}\{x*\tilde y\}=X(\omega)\cdot Y^*(\omega),
$$

同样得到相关定理。第二条推导凸显了“相关就是与共轭反转版本做卷积”这一结构，把新定理挂到已有结论上，只需一步。

对离散 DTFT 有完全对应的形式：

$$
\mathcal{F}_{\mathrm{DTFT}}\{R_{xy}[m]\}=X(e^{j\hat\omega})\,Y^*(e^{j\hat\omega}).
$$

DFT 版本则涉及循环相关，见下文。

#### Parseval 定理与零延迟相关

在上述 Fourier 约定下：

$$
R_{xy}(0)=\int x(t)y^*(t)\,dt
$$

满足

$$
\boxed{
R_{xy}(0)=\frac{1}{2\pi}
\int X(\omega)Y^*(\omega)\,d\omega
}.
$$

离散 DFT 下，若

$$
X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi kn/N},
$$

则

$$
\sum_{n=0}^{N-1}x[n]y^*[n]
=\frac1N\sum_{k=0}^{N-1}X[k]Y^*[k].
$$

### 时域点乘与频域循环卷积

频域普通逐点相乘对应时域卷积，而不是时域点对点乘法。对长度为 $N$ 的 DFT，若

$$
z[n]=x[n]y[n],
$$

将两个逆 DFT 代入并按频率索引收集，可得

$$
Z[k]=\frac{1}{N}\sum_{r=0}^{N-1}X[r]Y[(k-r)\bmod N].
$$

也就是说，时域点乘对应频域循环卷积：

$$
\boxed{\mathcal F_N\{x[n]y[n]\}=\frac1N(X\circledast_NY)[k]}.
$$

这与

$$
\operatorname{IDFT}\{X[k]Y[k]\}
$$

对应的循环卷积是两种不同的对偶关系。一个快速判断方法是看结果维度：点乘仍是与输入等长的逐样本向量，零延迟内积是标量，相关则是每个 lag 都有一个值的序列。

### 共轭方向、互谱与频域梯度

本文默认的相关定义为

$$
R_{xy}[m]=\sum_nx[n]y^*[n-m].
$$

因此

$$
R_{xy}[m]=\operatorname{IDFT}\{X[k]Y^*[k]\}[m].
$$

如果交换信号角色，或者需要矩阵梯度中的

$$
g[m]=\sum_nx^*[n-m]e[n],
$$

则频域形式为

$$
G[k]=X^*[k]E[k],
\qquad
g[m]=\operatorname{IDFT}\{X^*[k]E[k]\}[m].
$$

常见形式可以对照如下：

| 时域表达式 | 频域乘积 | 常见含义 |
|—|—|—|
| $\sum_nx[n]y^*[n-m]$ | $XY^*$ | $R_{xy}[m]$ 或互功率谱 |
| $\sum_ny[n]x^*[n-m]$ | $YX^*$ | $R_{yx}[m]$ |
| $\sum_nx^*[n-m]e[n]$ | $X^*E$ | $X^He$ 型频域梯度 |
| $\sum_nx[n]y^*[n]$ | $\frac1N\sum_kXY^*$ | 零延迟复内积 |
| $x[n]y[n]$ | $\frac1N X\circledast_NY$ | 时域点乘 |

对于实信号，频谱共轭对称可能使方向差异暂时不明显；对于复信号，$XY^*$ 与 $X^*Y$ 通常互为共轭并伴随 lag 方向变化，不能仅凭记忆互换。

### 相对相位与单频验证

设

$$
x[n]=e^{j\omega_0n},
\qquad
e[n]=Ae^{j(\omega_1n+\phi)}.
$$

共轭乘积为

$$
x^*[n]e[n]=Ae^{j((\omega_1-\omega_0)n+\phi)}.
$$

当 $\omega_1=\omega_0$ 时，乘积相位不再随时间旋转，累加幅度最大；普通乘积 $x[n]e[n]$ 则包含频率和 $\omega_0+\omega_1$。因此共轭的作用是比较相对相位和频率差，而不是把点乘“改写”为频域逐点乘法。
#### 功率谱与 Wiener–Khinchin 关系

自相关的 Fourier 变换为功率谱：

$$
\boxed{
S_{xx}(\omega)=\mathcal{F}\{R_{xx}(\tau)\}
=X(\omega)X^*(\omega)=|X(\omega)|^2
}.
$$

对于随机过程，这一结论称为 Wiener–Khinchin 定理：功率谱密度是自相关函数的 Fourier 变换。

互相关的 Fourier 变换则是互功率谱：

$$
S_{xy}(\omega)=X(\omega)Y^*(\omega).
$$

自功率谱通常为非负实数，而互功率谱一般是复数，其相位包含相对延迟信息。

自功率谱和互功率谱的完整推导、核心性质和数值算例见第五章。

—
### 线性卷积与循环卷积

#### 线性卷积

线性卷积假设信号在有效区间之外为零。对于有限序列 $x[n]$ 和 $h[n]$，输出长度为

$$
N+M-1.
$$

线性卷积没有首尾相接的回绕。

#### 循环卷积

长度为 $N$ 的循环卷积定义为

$$
y_N[n]=\sum_{k=0}^{N-1}x[k]h[(n-k)\bmod N].
$$

它把序列看成长度为 $N$ 的周期序列，因此一个序列末尾移出后会从开头重新出现。

#### DFT 与循环卷积

长度为 $N$ 的 DFT 满足

$$
\operatorname{DFT}\{x\circledast_N h\}
=X[k]H[k].
$$

所以直接对两个长度为 $N$ 的序列做 FFT、频域相乘再 IFFT，得到的是 $N$ 点循环卷积，而不一定是线性卷积。

#### 通过补零得到线性卷积

若 $x$ 长度为 $N$、$h$ 长度为 $M$，选择

$$
L\ge N+M-1
$$

并将二者补零至长度 $L$，则 $L$ 点循环卷积的前 $N+M-1$ 个样本等于线性卷积。

#### 循环相关

同理，直接计算

$$
\operatorname{IFFT}\{X[k]Y^*[k]\}
$$

得到的是循环相关。若需要线性相关，也必须足够补零，并根据相关定义截取有效区间。

#### 循环卷积与线性卷积的定量关系

设 $x[n]$、$h[n]$ 都在 $[0,N-1]$ 之外为零，长度均不大于 $N$。线性卷积 $y_L[n]=(x*h)[n]$ 在 $0\le n\le 2N-2$ 上非零；$N$ 点循环卷积 $y_C[n]=(x\circledast_N h)[n]$ 只覆盖 $0\le n\le N-1$。两者的关系为

$$
\boxed{
y_C[n]=\sum_{r\in\mathbb{Z}}y_L[n+rN],\qquad 0\le n\le N-1
}.
$$

即：**循环卷积等于把线性卷积按周期 $N$ 折叠回主区间。** 若线性卷积长度 $L_{\lin}=N_x+N_h-1\le N$，折叠不会产生重叠，$y_C$ 与 $y_L$ 在 $[0,N-1]$ 上完全一致；若 $L_{\lin}>N$，末端会“绕回”前端，形成时域混叠（time-domain aliasing）。因此“使用 FFT 求线性卷积”的黄金规则是：

$$
N\ge N_x+N_h-1,\qquad 通常取\ N=2^{\lceil\log_2 L_{\lin}\rceil}.
$$

对于相关，同样有

$$
R^{(C)}_{xy}[m]=\sum_r R^{(L)}_{xy}[m+rN],
$$

其中 $R^{(L)}$ 是线性相关。使用 FFT 求线性相关时也要按 $\ge N_x+N_y-1$ 补零，并把 IFFT 结果重排到中心 lag 处：负 lag 位于 $[N-\lfloor(N_y-1)\rfloor,N-1]$，正 lag 位于 $[0,N_x-1]$，可以用 `np.fft.fftshift` 或手工切片。

#### 最小 NumPy 示例：卷积与相关的三种算法互证

下面给出一个可以直接跑通的最小示例，比对“时域直接算”“FFT 加零算”“库函数”三种实现，工程上用作单元测试模板：

“`python
import numpy as np

rng = np.random.default_rng(0)
x = rng.standard_normal(8) + 1j * rng.standard_normal(8)
h = rng.standard_normal(5) + 1j * rng.standard_normal(5)
N_lin = len(x) + len(h) – 1

# 1) 时域直接线性卷积
y_direct = np.convolve(x, h, mode=”full”)

# 2) FFT 补零到 N >= N_lin
N = 1 << (N_lin - 1).bit_length() X = np.fft.fft(x, N) H = np.fft.fft(h, N) y_fft = np.fft.ifft(X * H)[:N_lin] # 3) 相关：R_xy[m] = sum_n x[n] * conj(y[n-m]) y_seq = h.copy() # 直接把 h 当作第二个信号 R_direct = np.correlate(x, y_seq, mode="full") # NumPy 使用 sum x[n+m] * conj(y[n]) R_fft = np.fft.ifft(np.fft.fft(x, N) * np.conj(np.fft.fft(y_seq, N))) R_fft = np.concatenate([R_fft[-(len(y_seq) - 1):], R_fft[: len(x)]]) print(np.allclose(y_direct, y_fft)) # True print(np.allclose(R_direct, R_fft)) # True ``` 关键点： - FFT 长度 $N$ 必须大于等于 $N_x+N_h-1$，否则最后一步的循环卷积会与线性卷积不一致； - 相关的 FFT 版本使用 $X\cdot Y^*$，最后要按定义把 IFFT 结果做 `fftshift` 或手工切分； - 用 `np.allclose` 做数值等价校验，容忍浮点误差，是自查代码正确性的最简单办法。 #### 分块卷积：Overlap-Add 与 Overlap-Save 简表 工程上遇到长信号 $x$（长度 $L\gg 1$）与固定短滤波器 $h$（长度 $M$）时，直接一次做长 FFT 会占用大量内存并引入延迟。分块卷积把 $x$ 切成长度 $L_{\text{blk}}$ 的段，每段做长度 $N=L_{\text{blk}}+M-1$（或 $L_{\text{blk}}\ge M$）的 FFT，将两种主流做法列表对照： | 特性 | Overlap-Add (OLA) | Overlap-Save (OLS) | |---|---|---| | 分块方式 | 把 $x$ 切成不重叠的长度 $L_{\text{blk}}$ 段 | 相邻段有 $M-1$ 个重叠样本 | | FFT 长度 | $N=L_{\text{blk}}+M-1$，各段补零后 FFT | $N=L_{\text{blk}}$，$N\ge M$，各段整体 FFT | | 输出拼接 | 相邻块尾部 $M-1$ 与下一块头部相加 | 丢弃每块前 $M-1$ 样本，剩余直接拼接 | | 时域混叠处理 | 通过补零消除 | 依赖“丢弃头部”把混叠部分裁掉 | | 计算量倾向 | 每块处理量略高，实现直观 | 每块少一次补零操作，常用高性能实现 | | 常见应用 | 一般数字滤波、教学示例 | 频域自适应滤波、GPU/DSP 高吞吐流水线 | 无论 OLA 还是 OLS，其核心都是循环卷积的**折叠公式** $y_C[n]=\sum_r y_L[n+rN]$：OLA 通过补零让折叠不产生重叠，OLS 通过丢弃前 $M-1$ 样本剔除受折叠污染的部分。理解了这个关系，就能在实现中随意切换两种策略而不出错。 --- --- ## 五、能量谱、功率谱与互谱 ### 能量谱的定义 对于能量信号，令 Fourier 变换为 $$ X(f)=\int_{-\infty}^{\infty}x(t)e^{-j2\pi ft}dt. $$ 能量谱密度（energy spectral density，ESD）定义为 $$ \boxed{ \Psi_x(f)=|X(f)|^2=X(f)X^*(f) }. $$ 它表示信号能量在频率轴上的分布。 如果使用角频率 $\omega$： $$ X(\omega)=\int x(t)e^{-j\omega t}dt, \qquad \Psi_x(\omega)=|X(\omega)|^2. $$ 必须注意 $f$ 和 $\omega$ 的变量不同，积分关系中会出现 $d\omega=2\pi df$ 的比例因子。 ### Parseval 定理与总能量 采用 $f$ 为频率的 Fourier 变换约定： $$ \boxed{ E_x=\int_{-\infty}^{\infty}|x(t)|^2dt =\int_{-\infty}^{\infty}|X(f)|^2df =\int_{-\infty}^{\infty}\Psi_x(f)df }. $$ 因此，能量谱密度曲线下面积等于信号总能量。 若采用角频率约定，则 $$ E_x=\frac{1}{2\pi}\int_{-\infty}^{\infty}|X(\omega)|^2d\omega. $$ $1/(2\pi)$ 不是额外的物理修正，而是 Fourier 正逆变换归一化约定带来的因子。 ### 能量谱密度与自相关 对能量信号，自相关定义为 $$ R_{xx}(\tau)=\int x(t)x^*(t-\tau)dt. $$ 相关定理给出 $$ \boxed{ \mathcal F\{R_{xx}(\tau)\}=|X(f)|^2=\Psi_x(f) }. $$ 上述结论可以从相关定理直接推出。自相关是互相关 $R_{xy}$ 在 $y=x$ 时的特殊情形： $$ R_{xx}(\tau)=\int x(t)x^*(t-\tau)\,dt. $$ 代入相关定理 $\mathcal F\{R_{xy}\}=X(f)Y^*(f)$，令 $Y=X$： $$ \mathcal F\{R_{xx}(\tau)\}=X(f)X^*(f)=|X(f)|^2. $$ 自功率谱之所以必然是实非负的，正是因为 $X\cdot X^*=|X|^2\ge 0$。这里的共轭 $X^*$ 来自自相关定义中对第二个信号取共轭 $x^*(t-\tau)$，而不是来自卷积——卷积的频域对应为 $X(f)H(f)$，不含共轭。 因此，能量谱密度是能量信号自相关函数的 Fourier 变换。反变换关系为 $$ R_{xx}(\tau)=\int_{-\infty}^{\infty} \Psi_x(f)e^{j2\pi f\tau}df. $$ 在零延迟处： $$ R_{xx}(0)=\int\Psi_x(f)df=E_x. $$ ### 能量谱的幅度和相位 能量谱密度只保留 Fourier 变换的幅度： $$ \Psi_x(f)=|X(f)|^2. $$ 因此它不包含 $X(f)$ 的绝对相位。不同信号可能具有相同的能量谱，但时域波形不同，这称为同谱或相位不确定性问题。 特别地，时移信号 $$ y(t)=x(t-t_0) $$ 的 Fourier 变换为 $$ Y(f)=X(f)e^{-j2\pi ft_0}, $$ 所以 $$ |Y(f)|^2=|X(f)|^2. $$ 能量谱无法单独确定信号的绝对时间位置。 ### 实信号的双边与单边能量谱 实值信号满足 $$ X(-f)=X^*(f), $$ 所以 $$ \Psi_x(-f)=\Psi_x(f). $$ 双边能量谱包含正、负频率。若只显示 $f\ge0$ 的单边能量谱，则除 DC 和 Nyquist 点外，通常将正频率部分乘以 2，以保持总能量： $$ E_x=\int_0^{\infty}\Psi_{x,\mathrm{one-sided}}(f)df. $$ 不能把双边谱直接截去负频率而不补偿，否则面积会少一半。 ### 离散有限序列的能量谱 对长度为 $N$ 的有限序列 $x[n]$，其 DTFT 为 $$ X(e^{j\omega})=\sum_{n=0}^{N-1}x[n]e^{-j\omega n}. $$ 能量谱密度为 $$ \Psi_x(e^{j\omega})=|X(e^{j\omega})|^2. $$ Parseval 关系为 $$ \boxed{ \sum_{n=0}^{N-1}|x[n]|^2 =\frac{1}{2\pi}\int_{-\pi}^{\pi} |X(e^{j\omega})|^2d\omega }. $$ 如果用 $N$ 点 DFT： $$ X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi kn/N}, $$ 则 $$ \boxed{ \sum_{n=0}^{N-1}|x[n]|^2 =\frac1N\sum_{k=0}^{N-1}|X[k]|^2 }. $$ ### 能量谱与频谱的区别 - $X(f)$ 是复频谱，包含幅度和相位； - $|X(f)|$ 是幅度谱； - $|X(f)|^2$ 是能量谱密度； - 三者的单位和积分意义不同。 能量谱不能恢复一般信号的相位，因此不能简单地把能量谱当作完整频谱。 --- ### 功率谱与功率谱密度 #### 功率谱密度的基本思想 功率信号的总能量通常为无穷大，不能使用 $$ |X(f)|^2 $$ 直接作为有限总能量谱。功率谱密度（power spectral density，PSD）描述平均功率如何分布在频率上。 对连续时间功率信号，可定义截断 Fourier 变换： $$ X_T(f)=\int_{-T/2}^{T/2}x(t)e^{-j2\pi ft}dt. $$ 功率谱密度定义为长时间归一化极限： $$ \boxed{ S_{xx}(f)=\lim_{T\to\infty} \frac{1}{T}\mathbb E\left[|X_T(f)|^2\right] }. $$ 对确定性功率信号，可在存在极限时省略期望；对随机过程，期望用于获得 ensemble 平均。 #### Wiener–Khinchin 定理 对宽平稳随机过程，功率谱密度是自相关函数的 Fourier 变换： $$ \boxed{ S_{xx}(f)=\mathcal F\{R_{xx}(\tau)\} =\int_{-\infty}^{\infty}R_{xx}(\tau)e^{-j2\pi f\tau}d\tau }. $$ 反变换为 $$ R_{xx}(\tau)=\int_{-\infty}^{\infty} S_{xx}(f)e^{j2\pi f\tau}df. $$ 零延迟给出平均功率： $$ \boxed{ P_x=R_{xx}(0)=\int_{-\infty}^{\infty}S_{xx}(f)df }. $$ 因此，PSD 曲线在整个频率轴上的面积等于平均功率，而不是总能量。 #### PSD 的单位 若 $x(t)$ 的单位为伏特，且采用电压平方表示功率，则 $S_{xx}(f)$ 的单位为 $$ \mathrm{V}^2/\mathrm{Hz}. $$ 若包含参考阻抗 $R$，物理功率谱密度可为 $$ S_P(f)=\frac{S_{xx}(f)}{R}, $$ 单位为 $$ \mathrm{W}/\mathrm{Hz}. $$ 对于离散数据，频率单位常为 cycles/sample 或 Hz，必须说明采样率换算。 #### 周期信号的功率谱 周期信号的 Fourier 变换通常由频率冲激组成，因此 PSD 是离散线谱而不是普通连续曲线。 设周期为 $T_0$ 的信号有 Fourier 级数 $$ x(t)=\sum_{k=-\infty}^{\infty}c_ke^{j2\pi kf_0t}, \qquad f_0=\frac1{T_0}. $$ 其平均功率为 $$ P_x=\sum_{k=-\infty}^{\infty}|c_k|^2. $$ 功率谱可以写为 $$ \boxed{ S_{xx}(f)=\sum_{k=-\infty}^{\infty} |c_k|^2\delta(f-kf_0) }. $$ 这表明每条谱线的面积是对应 Fourier 级数分量的平均功率。 #### 正弦信号的功率谱 对 $$ x(t)=A\cos(2\pi f_0t+\phi), $$ 有 $$ x(t)=\frac A2e^{j\phi}e^{j2\pi f_0t} +\frac A2e^{-j\phi}e^{-j2\pi f_0t}. $$ 其双边功率谱为 $$ S_{xx}(f)=\frac{A^2}{4}\delta(f-f_0) +\frac{A^2}{4}\delta(f+f_0). $$ 积分得到 $$ P_x=\frac{A^2}{2}. $$ 相位 $\phi$ 不影响单个正弦的功率谱，但会影响与其他信号的互功率谱相位。 #### 白噪声功率谱 理想连续时间白噪声的 PSD 为常数： $$ S_{xx}(f)=N_0/2 $$ 或根据双边/单边约定写为其他常数形式。白噪声在无限带宽上具有无限总功率，因此物理系统总要受到带宽限制。 经过带宽为 $B$ 的理想滤波器后，噪声功率为 $$ P=\int_{-B}^{B}S_{xx}(f)|H(f)|^2df. $$ 更一般地，LTI 系统输出 PSD 满足 $$ \boxed{ S_{yy}(f)=|H(f)|^2S_{xx}(f) }. $$ #### 离散时间功率谱 离散时间宽平稳过程的自相关为 $$ R_{xx}[m]=\mathbb E[x[n]x^*[n-m]]. $$ DTFT 给出 PSD： $$ \boxed{ S_{xx}(e^{j\omega}) =\sum_{m=-\infty}^{\infty}R_{xx}[m]e^{-j\omega m} }. $$ 平均功率为 $$ P_x=R_{xx}[0] =\frac{1}{2\pi}\int_{-\pi}^{\pi} S_{xx}(e^{j\omega})d\omega. $$ 若频率使用 $f$（Hz），则应将角频率与采样率联系起来： $$ \omega=2\pi f/f_s. $$ #### PSD 的非负性 自功率谱密度满足 $$ S_{xx}(f)\ge0. $$ 这是因为 PSD 是相关函数的 Fourier 变换，并且对于任意滤波器 $g$： $$ \mathbb E\left|\int g^*(t)x(t)dt\right|^2\ge0. $$ 由此可知，自谱不能出现真实的负功率。工程图中出现少量负值通常来自数值误差、去噪处理或绘图使用了不适合的估计量。 #### 自谱与能量谱的根本区别 | 项目 | 能量谱密度 ESD | 功率谱密度 PSD | |---|---|---| | 适用对象 | 能量信号 | 功率信号、随机过程、周期信号 | | 典型定义 | $|X(f)|^2$ | $\lim_{T\to\infty}|X_T(f)|^2/T$ 或 $\mathcal F\{R_{xx}\}$ | | 频率积分 | 总能量 $E_x$ | 平均功率 $P_x$ | | 单位 | 信号单位平方乘时间 | 信号单位平方/Hz，或功率/Hz | | 相位 | 不包含相位 | 自谱也不包含相位 | | 典型谱形 | 有限时长脉冲的连续谱 | 正弦的线谱、噪声的连续谱 | 同一有限数据段既可画 periodogram，也可按有限记录能量解释；但把其纵轴称为 PSD 还是 ESD，取决于归一化和问题对象，不能只看图形形状。 --- ### 互能量谱与互功率谱 #### 互能量谱 对于两个能量信号 $x(t)$、$y(t)$，定义互能量谱密度（cross energy spectral density）为 $$ \boxed{ \Psi_{xy}(f)=X(f)Y^*(f) }. $$ 它是互相关的 Fourier 表示： $$ \Psi_{xy}(f)=\mathcal F\{R_{xy}(\tau)\}. $$ 与自能量谱 $|X(f)|^2$ 不同，互能量谱一般是复数： $$ \Psi_{xy}(f)=|X(f)||Y(f)|e^{j(\angle X(f)-\angle Y(f))}. $$ 其模表示共同能量幅度，复数相位表示两个信号的相对相位，具体正负取决于互相关和共轭的定义。 #### 互能量谱的 Parseval 关系 零延迟互相关为 $$ R_{xy}(0)=\int x(t)y^*(t)dt. $$ 由 Parseval 定理： $$ \boxed{ R_{xy}(0)=\int\Psi_{xy}(f)df }. $$ 因此，互能量谱在频率上的复积分给出两个能量信号的内积。若存在时延，则对互能量谱做逆 Fourier 变换可以恢复不同延迟下的互相关。 #### 互功率谱 对联合宽平稳随机过程 $x(t)$、$y(t)$，互功率谱密度定义为互相关函数的 Fourier 变换： $$ \boxed{ S_{xy}(f)=\mathcal F\{R_{xy}(\tau)\} }. $$ 在截断 Fourier 估计中，常见形式为 $$ S_{xy}(f) =\lim_{T\to\infty}\frac{1}{T} \mathbb E\left[X_T(f)Y_T^*(f)\right]. $$ 互功率谱一般为复数，不能解释为“功率必须为正”的自谱。只有自功率谱 $S_{xx}$、$S_{yy}$ 必须非负；互功率谱可以有正负实部和虚部。 #### 互功率谱的推导与核心性质 **从截断 Fourier 变换推导互功率谱。** 对联合宽平稳过程 $x(t)$、$y(t)$，取观测区间 $[-T/2,T/2]$ 的截断 Fourier 变换 $$ X_T(f)=\int_{-T/2}^{T/2}x(t)e^{-j2\pi ft}dt, \qquad Y_T(f)=\int_{-T/2}^{T/2}y(t)e^{-j2\pi ft}dt. $$ 构造 $$ \frac{1}{T}\mathbb E\bigl[X_T(f)Y_T^*(f)\bigr] =\frac{1}{T}\int_{-T/2}^{T/2}\int_{-T/2}^{T/2} \mathbb E[x(t_1)y^*(t_2)]\,e^{-j2\pi f(t_1-t_2)}dt_1\,dt_2. $$ 由宽平稳性，$\mathbb E[x(t_1)y^*(t_2)]=R_{xy}(t_1-t_2)$。令 $\tau=t_1-t_2$，交换积分顺序并取 $T\to\infty$ 极限： $$ S_{xy}(f)=\lim_{T\to\infty}\frac{1}{T}\mathbb E[X_T(f)Y_T^*(f)] =\int_{-\infty}^{\infty}R_{xy}(\tau)e^{-j2\pi f\tau}d\tau =\mathcal F\{R_{xy}(\tau)\}. $$ 中间步骤的关键是 $Y_T^*$ 中的共轭继承自互相关定义 $R_{xy}(\tau)=\mathbb E[x(t)y^*(t-\tau)]$ 中的 $y^*$，与卷积无关。 **核心性质汇总：** 1. $S_{xy}(f)$ 一般为复数，其模 $|S_{xy}|$ 反映共同功率幅度，相位 $\angle S_{xy}$ 反映频率相关的相对相位。 2. $S_{yx}(f)=S_{xy}^*(f)$（交换信号角色等价于取共轭）。 3. $|S_{xy}(f)|^2\le S_{xx}(f)S_{yy}(f)$（Cauchy–Schwarz 不等式的频域形式），这正是相干性满足 $0\le\gamma_{xy}^2\le1$ 的根本原因。 4. 若 $y(t)=x(t-t_0)$（纯延迟），则 $$ R_{xy}(\tau)=R_{xx}(\tau-t_0), \qquad S_{xy}(f)=S_{xx}(f)\,e^{-j2\pi ft_0}. $$ 互谱相位为线性斜率 $-2\pi ft_0$，其斜率直接给出延迟 $t_0$。自谱 $S_{xx}$ 不包含这一信息。 5. 若 $x$ 与 $y$ 不相关（$R_{xy}(\tau)=0$ 对所有 $\tau$），则 $S_{xy}(f)=0$。不相关信号没有互谱贡献。 还要区分互谱与自适应滤波梯度。按本文约定，$S_{xe}=XE^*$ 表示输入与误差的互谱；Hermitian 最小二乘梯度通常出现 $X^*E$，它对应时域的 $\boldsymbol X^H\boldsymbol e$。二者都含共轭和相对相位，但信号顺序、相位符号和后续处理不同：互谱通常进行分段平均以描述统计关系，梯度则经过 IFFT、抽头截取和步长缩放来更新参数。 归一化频域梯度常写成 $$ \Delta W[k]=\mu\frac{X^*[k]E[k]}{\Phi_x[k]+\varepsilon}. $$ 当 $\Phi_x[k]$ 很小时，除法会放大噪声和舍入误差；因此 $\varepsilon$ 必须与功率量纲一致，并应结合功率平滑、更新限幅或频带门限使用。不能把 $XX^*$ 的幅度直接当作梯度方向，梯度方向还取决于误差的复相位。 #### 互谱的实部和虚部 写成 $$ S_{xy}(f)=S_{xy,R}(f)+jS_{xy,I}(f). $$ - 实部通常表示同相或偶对称相关成分； - 虚部通常表示正交或奇对称相关成分； - 模 $|S_{xy}|$ 表示互谱幅度； - 相位 $\arg S_{xy}$ 表示频率相关的相对相位。 实部和虚部的具体物理解释依赖信号定义、传感器方向、参考方向和互谱约定，不应脱离系统模型机械解释。 #### 互谱的共轭对称关系 互相关满足 $$ R_{yx}(\tau)=R_{xy}^*(-\tau). $$ Fourier 变换后： $$ \boxed{ S_{yx}(f)=S_{xy}^*(f) }. $$ 因此谱矩阵 $$ \boldsymbol S(f)= \begin{bmatrix} S_{xx}(f)&S_{xy}(f)\\ S_{yx}(f)&S_{yy}(f) \end{bmatrix} $$ 是 Hermitian 矩阵。它还必须是半正定的，这给出 $$ |S_{xy}(f)|^2\le S_{xx}(f)S_{yy}(f). $$ 这正是相干性满足 $0\le\gamma_{xy}^2(f)\le1$ 的根本原因。 #### LTI 系统中的互功率谱 若 $$ y(t)=h(t)*x(t)+n(t), $$ 且 $x$ 与 $n$ 不相关，则在采用 $$ S_{xy}=\mathbb E[X_TY_T^*]/T $$ 的约定下： $$ Y(f)=H(f)X(f)+N(f), $$ $$ \boxed{ S_{xy}(f)=H^*(f)S_{xx}(f) }. $$ 若采用相反的互谱定义 $S_{yx}=\mathbb E[Y_TX_T^*]/T$，则 $$ S_{yx}(f)=H(f)S_{xx}(f). $$ 因此传递函数估计可以写成两种对应形式： $$ \widehat H(f)=\frac{S_{yx}(f)}{S_{xx}(f)} $$ 或 $$ \widehat H^*(f)=\frac{S_{xy}(f)}{S_{xx}(f)}. $$ 实际使用时必须确认软件的 cross-spectrum 定义以及通道顺序。 #### 互谱与相干性的关系 互功率谱本身不能直接表示“相关比例”。需要使用 $$ \gamma_{xy}^2(f) =\frac{|S_{xy}(f)|^2}{S_{xx}(f)S_{yy}(f)}. $$ 互谱的绝对值可能很大，只是因为两个信号功率都大；相干性经过自谱归一化后才可以跨频率、跨实验比较线性关联强度。 #### 互谱相位与群时延 在高相干频带，互谱相位可用于估计时延。若系统频响相位为 $\phi(f)$，群时延定义为 $$ \tau_g(f)=-\frac{1}{2\pi}\frac{d\phi(f)}{df}. $$ 对于近似纯延迟，群时延接近常数；对于色散系统、滤波器或结构振动系统，群时延可能随频率变化。 相位缠绕、低相干性和频率响应零点都会使相位斜率不稳定，因此应联合查看相干性和自谱功率。 #### 互能量谱与互功率谱的对照 | 项目 | 互能量谱 $\Psi_{xy}$ | 互功率谱 $S_{xy}$ | |---|---|---| | 适用对象 | 能量信号、有限记录 | 功率信号、随机过程、长时间记录 | | 典型形式 | $X(f)Y^*(f)$ | $\lim_{T\to\infty}X_TY_T^*/T$ 的统计极限 | | 对应时域量 | 互相关的能量版本 | 互相关函数的功率版本 | | 频率积分 | 零延迟内积 | 零延迟互功率或协方差相关量 | | 是否可能为复数 | 是 | 是 | | 典型用途 | 脉冲、模板、有限波形比较 | 系统辨识、相干性、随机振动和噪声分析 | #### 互功率谱的数值算例 取一个长度为 $N=4$ 的复序列及其延迟版本： $$ x[n]=[1,\; j,\; -1,\; -j], \qquad y[n]=x[n-1]=[{-j},\; 1,\; j,\; -1]. $$ 其中 $y$ 是 $x$ 的循环右移 1 个样本。$x[n]$ 恰好是频率 $\omega_0=2\pi/4=\pi/2$ 的单频复指数 $e^{j\pi n/2}$。 做 4 点 DFT。由于 $x[n]=e^{j\pi n/2}$ 是单频信号，其 DFT 集中在一个 bin 上： $$ X[k]=\begin{cases}4,&k=1,\\0,&k\ne1.\end{cases} $$ 延迟 1 样本的 DFT 满足 $Y[k]=X[k]e^{-j2\pi k/4}$，即 $$ Y[k]=\begin{cases}4e^{-j\pi/2}=-4j,&k=1,\\0,&k\ne1.\end{cases} $$ 自功率谱： $$ S_{xx}[k]=X[k]X^*[k]=|X[k]|^2 =\begin{cases}16,&k=1,\\0,&k\ne1.\end{cases} $$ 结果为实非负，不含相位信息。 互功率谱： $$ S_{xy}[k]=X[k]Y^*[k] =\begin{cases}4\cdot(4j)=16j,&k=1,\\0,&k\ne1.\end{cases} $$ 结果为复数。在 $k=1$ 处，$\angle S_{xy}[1]=\pi/2=2\pi\cdot 1/4$，恰好是延迟 1 样本对应的相位斜率 $2\pi k/N$。 若错误地使用卷积形式 $X[k]Y[k]$（不取共轭）： $$ X[k]Y[k] =\begin{cases}4\cdot(-4j)=-16j,&k=1,\\0,&k\ne1.\end{cases} $$ 其相位为 $-\pi/2$，等于 $\angle X[1]+\angle Y[1]=0+(-\pi/2)$，反映的是两个信号相位之和，而不是相位差。这不能用于时延估计。 这个例子清楚地表明： - 自功率谱 $|X|^2$ 丢弃所有相位，只保留能量分布； - 互功率谱 $XY^*$ 通过共轭把相位变成差值 $\angle X-\angle Y$，编码两个信号的相对延迟； - 共轭来自互相关定义中的 $y^*$，不是卷积的组成部分。 --- ### ESD、PSD 与互谱的统一比较 前几节分别定义了能量谱密度 (ESD)、功率谱密度 (PSD)、互能量谱和互功率谱。它们在时域—频域关系上有一个统一骨架： | 时域量 | 对应频域量 | 频率积分给出 | |---|---|---| | $R_{xx}(\tau)=\int x(t)x^*(t-\tau)dt$ | $\Psi_x(f)=|X(f)|^2$ | 总能量 $E_x$ | | $R_{xx}(\tau)=\mathbb E[x(t)x^*(t-\tau)]$ | $S_{xx}(f)$ | 平均功率 $P_x$ | | $R_{xy}(\tau)=\int x(t)y^*(t-\tau)dt$ | $\Psi_{xy}(f)=X(f)Y^*(f)$ | $R_{xy}(0)$ | | $R_{xy}(\tau)=\mathbb E[x(t)y^*(t-\tau)]$ | $S_{xy}(f)$ | $R_{xy}(0)$ | 因此这些谱本质上都是相关函数的 Fourier 变换。差别只在于对象是能量信号还是功率信号／随机过程、是否需要期望，以及归一化中是否除以时长 $T$。对同一段有限记录，若把归一化因子由 $1$ 换成 $1/T$，就可以在“记录能量”与“记录平均功率”两种观点间转换： $$ \widehat S_{xx}(f)\approx\frac{1}{T}\widehat\Psi_x(f), \qquad \widehat S_{xy}(f)\approx\frac{1}{T}\widehat\Psi_{xy}(f). $$ 一旦记录中含有窗函数或有效带宽被修改，还需要额外的窗能量归一化，这在第八章详细讨论。 --- ### 采样与归一化：一个数值算例 考虑离散正弦 $x[n]=A\cos(2\pi f_0 n/f_s)$，取 $N$ 点使得 $f_0 N/f_s$ 恰为整数（无泄漏），做 $N$ 点 DFT： $$ X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi kn/N}. $$ 在 $k=\pm k_0$（$k_0=f_0N/f_s$）两个 bin 上，$|X[k_0]|=|X[-k_0]|=NA/2$，其他 bin 为零。三种常见归一化的意义如下： | 表达式 | 数值 | 意义 | |---|---|---| | $\sum_k|X[k]|^2$ | $N^2A^2/2$ | 未归一化 DFT 能量 | | $\dfrac1N\sum_k|X[k]|^2$ | $NA^2/2$ | 满足 Parseval：等于 $\sum_n|x[n]|^2$ | | $\dfrac{1}{N^2}\sum_k|X[k]|^2$ | $A^2/2$ | 等于该正弦的平均功率 $P_x$ | 若把 $|X[k]|^2/N$ 视为 periodogram 频点数据，则两条谱线各自的“功率”均为 $A^2/4$；加起来才等于双边总功率 $A^2/2$。若使用单边谱，需要把非 DC、非 Nyquist 处的功率翻倍，得到 $k_0$ 处功率为 $A^2/2$，与单边线谱一致。 上述结果与是否使用角频率、是否含 $2\pi$、以及是否除以 $f_s$ 都相关。粗略的自检方法是：把估计谱在整个正频段（单边约定）积分或求和乘频率间隔，看是否恢复出理论功率 $A^2/2$。不能仅凭曲线形状判断纵轴是幅度谱、能量谱还是功率谱。 --- ### 谱矩阵的 Hermitian 与半正定性质 对 $M$ 个联合宽平稳过程 $x_1(t),\dots,x_M(t)$，将其在频率 $f$ 处的所有自谱和互谱排为 $M\times M$ 谱矩阵 $$ \boldsymbol S(f)= \begin{bmatrix} S_{11}(f)&\cdots&S_{1M}(f)\\ \vdots&\ddots&\vdots\\ S_{M1}(f)&\cdots&S_{MM}(f) \end{bmatrix}, \qquad S_{ij}(f)=\mathcal F\{R_{ij}(\tau)\}. $$ 由 $R_{ji}(\tau)=R_{ij}^*(-\tau)$ 可得 $S_{ji}(f)=S_{ij}^*(f)$，即 $$ \boldsymbol S(f)=\boldsymbol S(f)^{\mathrm H}. $$ 对任意复向量 $\boldsymbol a\in\mathbb C^M$，令 $z(t)=\sum_i a_i^*x_i(t)$，其功率谱 $$ S_{zz}(f)=\boldsymbol a^{\mathrm H}\boldsymbol S(f)\boldsymbol a\ge0. $$ 因此谱矩阵在每个频率上都是 Hermitian 半正定的。它的特征分解 $$ \boldsymbol S(f)=\sum_{k=1}^{M}\lambda_k(f)\boldsymbol u_k(f)\boldsymbol u_k^{\mathrm H}(f) $$ 具有直接的物理意义：$\lambda_k(f)$ 是该频率上第 $k$ 个正交“空间模式”所占的功率，$\boldsymbol u_k(f)$ 是对应的复合成方向。这是多通道谱主成分分析、盲源分离和阵列信号处理的基础。 两两情形回到 $$ \det\boldsymbol S(f) =S_{xx}(f)S_{yy}(f)-|S_{xy}(f)|^2\ge0, $$ 正是第七章相干性上界 $\gamma_{xy}^2(f)\le1$ 的谱矩阵版本。 --- ### LTI 系统频率响应估计与方向性 在噪声环境下，可以用互谱和自谱构造多种频率响应函数（FRF）估计。对模型 $y(t)=h(t)*x(t)+n(t)$（噪声只加在输出，$x\perp n$），最常用的三种估计器为 $$ \widehat H_1(f)=\frac{S_{yx}(f)}{S_{xx}(f)}, \qquad \widehat H_2(f)=\frac{S_{yy}(f)}{S_{xy}(f)}, \qquad \widehat H_v(f)=\sqrt{\widehat H_1(f)\,\widehat H_2(f)}. $$ 这里采用全文约定 $S_{xy}=\mathbb E[X Y^*]$，因此对 $y=H*x+n$ 且输入与噪声不相关，有 $S_{yx}=H S_{xx}$；若软件把互谱通道顺序定义为 $S_{xy}=\mathbb E[Y X^*]$，则相应公式中的下标需要整体交换。 它们的偏差性质与噪声位置有关： - $\widehat H_1$：在输入无噪、输出有噪时无偏；输入端有噪时会低估幅值； - $\widehat H_2$：在输出无噪、输入有噪时无偏；输出端有噪时会高估幅值； - $\widehat H_v$：几何平均，介于两者之间，常在两端都有噪声时给出更稳健结果。 三者与相干性满足严格关系 $$ \boxed{ \frac{\widehat H_1(f)}{\widehat H_2(f)} =\widehat\gamma_{xy}^2(f) }. $$ 因此 $H_1$ 与 $H_2$ 的比值差异本身就是相干性缺陷的度量。选择 FRF 估计器时应结合噪声主导端和相干性形状，不能只按“看起来更平滑”决定。 方向性方面，$\widehat H_1$ 依赖“输入→输出”方向定义。反过来把角色互换，会得到 $1/\widehat H_1$ 或 $\widehat H_2$ 的复共轭之类的量，其偏差方向也随之翻转。在多输入、多输出场合，需要用 $\boldsymbol H(f)=\boldsymbol S_{yx}(f)\boldsymbol S_{xx}^{-1}(f)$ 等矩阵形式来推广，且要求 $\boldsymbol S_{xx}$ 在感兴趣频率上良态可逆。 --- --- ## 六、相关性与统计信号分析 “相关”常指相关函数这一具体运算，而“相关性”通常是更宽泛的概念：它描述两个量是否一起变化，以及这种共同变化有多强。不同学科中“相关性”可能指内积、互相关函数、协方差、Pearson 相关系数或其他依赖性指标，因此必须先明确对象和定义。 ### 随机变量的期望、方差与协方差 设随机变量 $X$ 和 $Y$ 的均值分别为 $$ \mu_X=\mathbb E[X], \qquad \mu_Y=\mathbb E[Y]. $$ 方差为 $$ \sigma_X^2=\mathbb E\left[|X-\mu_X|^2\right], \qquad \sigma_Y^2=\mathbb E\left[|Y-\mu_Y|^2\right]. $$ 复随机变量常采用以下互协方差定义： $$ \boxed{ C_{XY}=\mathbb E\left[(X-\mu_X)(Y-\mu_Y)^*\right] }. $$ 它去除了均值，测量两个变量围绕各自均值的共同线性变化。若不去均值，则得到互相关矩： $$ R_{XY}=\mathbb E[XY^*]. $$ 二者满足 $$ C_{XY}=R_{XY}-\mu_X\mu_Y^*. $$ 因此，“相关为零”和“协方差为零”只在至少一个变量为零均值时等价。 ### Pearson 相关系数 对实随机变量，Pearson 相关系数定义为 $$ \boxed{ \rho_{XY} =\frac{\operatorname{Cov}(X,Y)}{\sigma_X\sigma_Y} }. $$ 只要 $\sigma_X>0$、$\sigma_Y>0$，就有

$$
-1\le\rho_{XY}\le1.
$$

其典型解释为：

– $\rho_{XY}=1$：完全正线性关系；
– $\rho_{XY}=-1$：完全负线性关系；
– $\rho_{XY}=0$：无 Pearson 线性相关，但仍可能存在非线性依赖；
– $|\rho_{XY}|$ 越接近 1，线性关联越强。

对复随机变量，可以定义复相关系数

$$
\rho_{XY}
=\frac{C_{XY}}{\sqrt{C_{XX}C_{YY}}},
\qquad
|\rho_{XY}|\le1.
$$

此时 $|\rho_{XY}|$ 表示线性关联强度，$\angle\rho_{XY}$ 描述平均相对相位。

### 样本相关系数

给定 $N$ 对实值观测 $(x_i,y_i)$，样本均值为

$$
\bar x=\frac1N\sum_{i=1}^{N}x_i,
\qquad
\bar y=\frac1N\sum_{i=1}^{N}y_i.
$$

样本 Pearson 相关系数为

$$
\boxed{
r_{xy}
=\frac{\sum_{i=1}^{N}(x_i-\bar x)(y_i-\bar y)}
{\sqrt{\sum_{i=1}^{N}(x_i-\bar x)^2}
\sqrt{\sum_{i=1}^{N}(y_i-\bar y)^2}}
}.
$$

它是总体相关系数的估计。有限样本中的 $r_{xy}$ 会随机波动，因此不能只凭一个数值判断总体关系；通常还要考虑样本量、置信区间、异常值和数据分布。

### 相关系数的性质

若 $a,c>0$，则

$$
\operatorname{corr}(aX+b,cY+d)=\operatorname{corr}(X,Y).
$$

平移和正比例缩放不改变 Pearson 相关系数；负比例缩放会改变符号但不改变绝对值。

相关系数没有物理单位，这是其便于比较的原因；协方差则保留两个变量单位的乘积。

Pearson 相关系数只测量线性关系。例如令 $X$ 关于零对称，$Y=X^2$，则 $Y$ 完全由 $X$ 决定，但可能有

$$
\operatorname{Cov}(X,Y)=\mathbb E[X^3]=0,
$$

从而 Pearson 相关系数为零。

### 独立、不相关与正交

三个概念需要严格区分：

1. **统计独立**：联合分布可以分解为边缘分布的乘积；
2. **不相关**：$C_{XY}=0$；
3. **正交**：在给定内积下 $\langle x,y\rangle=0$。

对于二阶矩存在的随机变量：

$$
\text{独立}\Longrightarrow\text{不相关},
$$

但一般没有反向推论。若 $X$、$Y$ 联合 Gaussian，则不相关可推出独立。

在随机信号理论中，如果采用

$$
\langle X,Y\rangle=\mathbb E[XY^*],
$$

那么零均值随机变量的“不相关”与这种均方意义下的“正交”一致。

### 随机过程的相关性

对随机过程 $x(t)$ 和 $y(t)$，互相关定义为

$$
R_{xy}(t_1,t_2)
=\mathbb E\left[x(t_1)y^*(t_2)\right].
$$

若过程是联合宽平稳的，互相关只依赖时间差：

$$
R_{xy}(\tau)
=\mathbb E\left[x(t)y^*(t-\tau)\right].
$$

对应的互协方差为

$$
C_{xy}(\tau)
=R_{xy}(\tau)-\mu_x\mu_y^*.
$$

宽平稳过程的均值不随时间变化，自相关只依赖时差，并满足

$$
R_{xx}(-\tau)=R_{xx}^*(\tau).
$$

### 时间平均与集合平均

理论定义中的 $\mathbb E[\cdot]$ 是集合平均：需要对同一随机实验的许多实现求平均。但工程中往往只有一条有限记录，于是用时间平均估计：

$$
\widehat R_{xy}(\tau)
=\frac1T\int_0^T x(t)y^*(t-\tau)\,dt.
$$

只有在适当的遍历性条件下，长时间平均才收敛到集合平均。平稳性并不自动保证遍历性，因此用单条记录估计统计量时，应明确这一假设。

### Pearson、Spearman 与非线性依赖

常见相关性指标包括：

| 指标 | 测量对象 | 主要特点 |
|—|—|—|
| Pearson 相关 | 线性关系 | 对异常值敏感；要求二阶矩存在 |
| Spearman 秩相关 | 单调关系 | 基于秩，对非线性单调关系和异常值更稳健 |
| Kendall $\tau$ | 次序一致性 | 解释为样本对一致与不一致的差异 |
| 互信息 | 一般统计依赖 | 能发现非线性依赖，但估计更困难 |
| 距离相关 | 一般依赖 | 在一定条件下为零当且仅当独立 |

信号与系统中最常用的是互相关函数和谱相干性，但当关系明显非线性时，Pearson 相关或普通相干性可能不足。

### 相关不等于因果

即使两个信号的相关性很高，也不能仅凭相关推出因果。常见原因包括：

– 两者由同一个隐藏输入驱动；
– 一个共同趋势造成伪相关；
– 两个测量通道存在串扰；
– 反馈使影响方向双向存在；
– 时间序列自相关降低了有效独立样本数。

因果判断通常还需要实验设计、物理模型、时间方向、控制混杂变量，以及专门的因果推断方法。

—
### 平稳性与遍历性的层次

在用有限样本估计总体统计量之前，必须明确随机过程的假设层次。常用概念按由弱到强排列为：

1. **一阶平稳**：$\mathbb E[x(t)]$ 与 $t$ 无关；
2. **宽平稳 (WSS)**：一阶平稳，且 $R_{xx}(t,t-\tau)$ 只依赖于 $\tau$；
3. **联合宽平稳**：多个过程都 WSS，且互相关 $R_{xy}(t,t-\tau)$ 只依赖 $\tau$；
4. **严格平稳**：所有阶联合分布对时间平移不变；
5. **均值/相关遍历**：时间平均以概率或均方收敛到集合平均；
6. **严格遍历**：所有可积函数的时间平均等于集合平均。

WSS 只是二阶矩层面的条件，并不自动保证遍历性：一个宽平稳过程可能在不同实现之间取值分布不同（例如“抽签固定”的随机常数过程 $x(t)=A$）。工程上使用

$$
\widehat R_{xx}(\tau)
=\frac1T\int_0^T x(t)x^*(t-\tau)dt
$$

去估计 $R_{xx}(\tau)$ 时，实际上暗含 **均值遍历、协方差遍历** 之类的假设：需要一段足够长的记录内，样本充分“探索”了随机过程的可能状态。若被测系统包含慢时变、模式切换或非各态历经的稳态，长记录时间平均并不逼近集合平均，谱和相干性估计就会出现系统性偏差。

一个简单判断做法是把长记录切成若干互不重叠的子段，分别估计均值、方差或谱，比较它们的差异是否在期望的统计涨落范围内。差异远超统计涨落时，遍历性假设应被认为不成立，应改按工况分段处理，而不是继续增加平均次数。

—
### 样本方差与协方差的分母：$N$ 还是 $N-1$

样本 Pearson 相关系数

$$
r_{xy}
=\frac{\sum_i(x_i-\bar x)(y_i-\bar y)}
{\sqrt{\sum_i(x_i-\bar x)^2}\sqrt{\sum_i(y_i-\bar y)^2}}
$$

对分子分母共用同一个 $N$ 或 $N-1$ 因子时结果不变，因此常见教科书写法看似有歧义，但比值本身与分母选择无关。真正需要小心的是**方差和协方差本身**的估计：

$$
\widehat\sigma_X^2=\frac1N\sum_{i=1}^N(x_i-\bar x)^2
\quad\text{（最大似然、偏估计）},
$$

$$
s_X^2=\frac{1}{N-1}\sum_{i=1}^N(x_i-\bar x)^2
\quad\text{（无偏估计，Bessel 修正）}.
$$

对协方差同理，$N-1$ 分母抵消“均值也是估计出的”这一自由度损失。若使用 $N$ 分母，会系统性低估方差和协方差；如果之后再基于这些量做卡方检验、F 检验或置信区间，就会得到过度自信的结论。

在信号处理软件中：

– Python `numpy.var/cov` 默认 `ddof=0`（分母 $N$），需要显式 `ddof=1` 才对应 $N-1$；
– Python `numpy.corrcoef` 内部会归一化，因此不受该差异影响；
– MATLAB `var/cov` 默认使用 $N-1$；`var(x,1)` 才是 $N$。

在报告相关性、协方差矩阵、白化因子或 Mahalanobis 距离时，必须写清所用分母，以免下游把偏估计当作无偏估计使用。

—
### $Y=X^2$：非线性依赖使相关系数消失

第六章正文提到，$Y=X^2$ 时 Pearson 相关系数可能为零。这里给出一个可直接验证的算例：设 $X\sim\mathcal U[-1,1]$，则

$$
\mu_X=\mathbb E[X]=0,
\qquad
\mu_Y=\mathbb E[X^2]=\frac13,
$$

$$
\operatorname{Cov}(X,Y)
=\mathbb E[XY]-\mu_X\mu_Y
=\mathbb E[X^3]=0,
$$

因此 $\rho_{XY}=0$。但 $Y$ 完全由 $X$ 决定，两者显然存在函数依赖。

同样的构造在信号中经常出现。例如输入正弦 $x(t)=\cos(2\pi f_0 t)$ 经过平方器 $y(t)=x^2(t)$ 得到

$$
y(t)=\tfrac12+\tfrac12\cos(2\pi\cdot 2f_0\cdot t).
$$

在同频 $f_0$ 上 $y$ 无功率，因此 MSC $\gamma_{xy}^2(f_0)$ 接近 0；能量转移到直流和二倍频。在 $2f_0$ 上 $y$ 有强功率，但 $x$ 没有，二阶谱相干性同样趋于 0。这类跨频关系需要用双相干性 (bicoherence)、非线性系统辨识或互信息才能揭示。

结论：Pearson 相关系数和二阶谱相干性只对**同频线性关系**敏感。当怀疑存在整流、平方、绝对值、包络、幅相耦合等非线性时，应改用适当的非线性指标或引入解析信号／包络后再做二阶分析。

—
### 共同驱动源与伪相关

即使 $x$ 与 $y$ 之间没有直接联系，若两者被同一个信号 $z$ 驱动，它们之间的相关系数或相干性都可能显著非零。举一个显式例子：设 $z$ 为零均值单位方差的宽平稳过程，$n_1,n_2$ 是与 $z$ 及彼此都不相关的白噪声，令

$$
x(t)=z(t)+n_1(t),
\qquad
y(t)=z(t)+n_2(t).
$$

则

$$
R_{xy}(\tau)=R_{zz}(\tau),
\qquad
S_{xy}(f)=S_{zz}(f).
$$

而

$$
S_{xx}(f)=S_{zz}(f)+S_{n_1n_1}(f),
\qquad
S_{yy}(f)=S_{zz}(f)+S_{n_2n_2}(f).
$$

因此普通相干性

$$
\gamma_{xy}^2(f)
=\frac{S_{zz}^2(f)}
{(S_{zz}(f)+S_{n_1n_1}(f))(S_{zz}(f)+S_{n_2n_2}(f))}
$$

会在 $z$ 主导的频带显著大于 0，但 $x$、$y$ 之间并不存在直接因果链。这就是“共同驱动源伪相关”的一种最简单形式。

要在数据分析中区分“直接联系”与“共同驱动”，通常需要：

– 引入并测量候选共同源 $z$，使用第七章的**偏相干性**将其影响回归掉；
– 引入外部扰动或干预实验，人为改变一端而不改变共同源；
– 结合物理模型或时间滞后结构，例如 Granger 因果、结构 VAR。

单凭高相关或高相干得出因果结论几乎总是错误的；把 $z$ 测出来并纳入分析，是最基本也最有效的排除步骤。

—

—

## 七、相干性及其扩展

“相干性”在不同语境下有不同含义。本节重点介绍随机信号和谱分析中最常用的 **magnitude-squared coherence，MSC**。它描述两个信号在每一个频率上的线性关联强度，可以看成“频率分辨的平方相关系数”。

### 互功率谱和自功率谱

对联合宽平稳过程，定义互功率谱密度

$$
S_{xy}(f)=\mathcal F\{R_{xy}(\tau)\}.
$$

自功率谱密度为

$$
S_{xx}(f)=\mathcal F\{R_{xx}(\tau)\},
\qquad
S_{yy}(f)=\mathcal F\{R_{yy}(\tau)\}.
$$

互功率谱通常是复数：

$$
S_{xy}(f)=|S_{xy}(f)|e^{j\phi_{xy}(f)}.
$$

其中：

– $|S_{xy}(f)|$ 表示该频率上共同变化的幅度；
– $\phi_{xy}(f)$ 表示该频率上的平均相位差；
– $S_{xx}(f)$、$S_{yy}(f)$ 表示各信号在该频率附近的功率密度。

### 复相干度与平方相干性

归一化复互谱定义为复相干度：

$$
\boxed{
C_{xy}(f)
=\frac{S_{xy}(f)}{\sqrt{S_{xx}(f)S_{yy}(f)}}
}.
$$

其模平方称为 magnitude-squared coherence：

$$
\boxed{
\gamma_{xy}^2(f)
=\frac{|S_{xy}(f)|^2}{S_{xx}(f)S_{yy}(f)}
}.
$$

由谱矩阵的半正定性或 Cauchy–Schwarz 不等式：

$$
0\le\gamma_{xy}^2(f)\le1.
$$

解释如下：

– $\gamma_{xy}^2(f)\approx1$：在该频率上存在稳定、可重复的强线性关系；
– $\gamma_{xy}^2(f)\approx0$：在该频率上没有稳定的线性关联，或者关联被不相关噪声、非线性、时变关系或估计误差削弱；
– 中间值：只有部分输出功率能由与另一个信号线性相关的成分解释。

复相干度还保留相位：

$$
\angle C_{xy}(f)=\angle S_{xy}(f).
$$

而 MSC 只有 $0$ 到 $1$ 的实数强度，不保留相位符号。

### 为什么相干性需要归一化

互功率谱 $S_{xy}(f)$ 的大小同时受到两个信号幅度影响，不便直接比较不同频率或实验。除以

$$
\sqrt{S_{xx}(f)S_{yy}(f)}
$$

后，相干性变成无量纲指标，并限制在 $[0,1]$ 内。

这与时域相关系数

$$
\rho_{XY}=\frac{C_{XY}}{\sigma_X\sigma_Y}
$$

的归一化思想相同：相关系数按总方差归一化，相干性则逐频率按谱功率归一化。

### 相干性与频率响应估计

考虑线性系统模型

$$
y(t)=h(t)*x(t)+n(t),
$$

其中噪声 $n(t)$ 与输入 $x(t)$ 不相关。频域中

$$
Y(f)=H(f)X(f)+N(f).
$$

于是

$$
S_{xy}(f)=H^*(f)S_{xx}(f)
$$

或在相反互谱定义下为 $H(f)S_{xx}(f)$。同时

$$
S_{yy}(f)=|H(f)|^2S_{xx}(f)+S_{nn}(f).
$$

代入 MSC：

$$
\gamma_{xy}^2(f)
=\frac{|H(f)|^2S_{xx}(f)}
{|H(f)|^2S_{xx}(f)+S_{nn}(f)}.
$$

因此，在上述理想模型下，相干性可以解释为该频率输出功率中由输入的线性响应贡献的比例。若输出噪声为零，则相干性为 1；噪声越强，相干性越低。

### 相干性与信噪比

定义输出端该频率的信噪比为

$$
\operatorname{SNR}(f)
=\frac{|H(f)|^2S_{xx}(f)}{S_{nn}(f)}.
$$

则

$$
\boxed{
\gamma_{xy}^2(f)
=\frac{\operatorname{SNR}(f)}{1+\operatorname{SNR}(f)}
}.
$$

反过来：

$$
\boxed{
\operatorname{SNR}(f)
=\frac{\gamma_{xy}^2(f)}{1-\gamma_{xy}^2(f)}
}.
$$

这只在单输入线性系统、噪声与输入不相关等假设成立时才适用，不能把它当成所有场景下的通用 SNR 公式。

### 相位谱、相位锁定与时延

互谱相位为

$$
\phi_{xy}(f)=\arg S_{xy}(f).
$$

若两个信号仅相差固定时延 $\tau_0$：

$$
y(t)=x(t-\tau_0),
$$

则根据互谱定义，互谱相位通常呈线性变化：

$$
\phi_{xy}(f)=\pm2\pi f\tau_0.
$$

正负号取决于互谱定义。由相位斜率可估计时延：

$$
\tau_0=\pm\frac{1}{2\pi}\frac{d\phi_{xy}}{df}.
$$

使用这一方法前应先进行 phase unwrapping，并只在相干性足够高的频带解释相位；低相干频点上的相位通常不稳定、没有可靠物理意义。

### 相干性为 1 的含义和限制

理论上的 $\gamma_{xy}^2(f)=1$ 表示，在二阶统计意义下，该频率处一个信号可由另一个信号经过确定的线性频率响应完全解释。但它不意味着：

– 两个时域波形完全相同；
– 频率响应的幅值必须为 1；
– 相位差必须为 0；
– 一个信号必然因果地导致另一个信号；
– 系统在所有频率上都线性。

例如，理想无噪声 LTI 系统可以任意改变幅值和相位，只要输出完全由输入线性产生，对应频率上的相干性仍可为 1。

### 相干性降低的原因

低相干性可能由以下因素造成：

1. 输出中存在与输入不相关的噪声；
2. 输入测量本身有噪声；
3. 还有未测量的其他输入影响输出；
4. 输入与输出之间存在非线性关系；
5. 系统参数随时间变化；
6. 输入与输出没有正确对时；
7. 数据段过短，谱估计方差过大；
8. 泄漏、混叠或传感器饱和；
9. 频率分辨率不合适，把不同动力学混在一个 bin；
10. 不同数据段中的相位关系不稳定，平均后互谱相互抵消。

因此，低相干性不是“完全没有关系”的充分证据，而是“没有检测到稳定线性二阶关联”或估计条件不足。

### 相干函数与脉冲相干性的区别

“coherence”还可能指其他概念：

– **谱相干性**：本节定义的 $\gamma_{xy}^2(f)$；
– **复相干度**：保留相位的 $C_{xy}(f)$；
– **相位锁定值（PLV）**：只分析跨试次的相位差稳定性，不直接考虑幅度；
– **光学相干性**：描述电磁场在时间或空间上的相位关联，有一阶、二阶相干函数；
– **小波相干性**：在时间—频率平面上分析局部共同变化；
– **偏相干性、多重相干性**：控制其他通道或处理多个输入。

这些量虽然都涉及“稳定关系”，但定义、取值和统计性质不同，不能混用。

### 偏相干性

若 $x$ 和 $y$ 都受第三个信号 $z$ 驱动，普通相干性可能很高，但这不代表 $x$ 与 $y$ 存在直接联系。偏相干性用于在控制 $z$ 后分析剩余线性关联。

以复相干度表示，一种形式为

$$
C_{xy\cdot z}(f)
=\frac{C_{xy}(f)-C_{xz}(f)C_{zy}(f)}
{\sqrt{\left(1-|C_{xz}(f)|^2\right)
\left(1-|C_{yz}(f)|^2\right)}}.
$$

其模平方给出偏平方相干性。该公式适用于相应谱矩阵可逆且统计估计可靠的场景；多变量情况下通常直接使用谱矩阵及其逆矩阵计算。

### 多重相干性

多重相干性衡量多个输入共同线性解释某个输出的能力。设输入向量为 $\boldsymbol x$，输出为 $y$，其谱矩阵为 $\boldsymbol S_{xx}$，输入与输出的互谱向量为 $\boldsymbol S_{xy}$，则多重相干性可写为

$$
\gamma_{y:\boldsymbol x}^2(f)
=\frac{\boldsymbol S_{yx}(f)
\boldsymbol S_{xx}^{-1}(f)
\boldsymbol S_{xy}(f)}{S_{yy}(f)}.
$$

它常用于多输入单输出系统辨识、振动分析和多传感器数据融合。

### 小波相干性

普通 MSC 假设统计关系在分析时段内近似稳定，只给出频率维度。如果信号明显非平稳，可以使用小波相干性，在时间—尺度或时间—频率平面中观察局部相干。

小波相干性的一般形式为

$$
R_{xy}^2(a,b)
=\frac{|S\{a^{-1}W_{xy}(a,b)\}|^2}
{S\{a^{-1}|W_x(a,b)|^2\}
S\{a^{-1}|W_y(a,b)|^2\}},
$$

其中 $W_x$、$W_y$ 为连续小波变换，$W_{xy}=W_xW_y^*$，$S\{\cdot\}$ 表示适当的时间和尺度平滑。没有平滑时同样可能产生退化结果。

—
### 相关性与相干性的区别和联系

#### 一张表看清区别

| 比较项 | 时域相关/相关系数 | 谱相干性 |
|—|—|—|
| 自变量 | 延迟 $\tau$，或单个总体统计量 | 频率 $f$ |
| 典型定义 | $R_{xy}(\tau)$、$\rho_{XY}$ | $\gamma_{xy}^2(f)$ |
| 是否归一化 | 互相关通常不归一化；相关系数归一化 | 按频率归一化 |
| 取值 | 互相关可有单位、可为复数；实 Pearson 在 $[-1,1]$ | MSC 在 $[0,1]$ |
| 主要回答 | 是否共同变化、最佳时延在哪里 | 哪些频率具有稳定线性关系 |
| 相位信息 | 复互相关可包含整体延迟/相位信息 | 复相干度或互谱相位逐频率给出相位 |
| 对非平稳性的处理 | 滑动相关可做局部分析 | 普通 MSC 需近似平稳；可改用短时或小波相干性 |
| 对因果性的结论 | 不能单独证明因果 | 同样不能单独证明因果 |

#### 相干性是频率分辨的归一化关联

从数学结构看：

$$
\rho_{XY}
=\frac{\mathbb E[(X-\mu_X)(Y-\mu_Y)^*]}
{\sigma_X\sigma_Y}
$$

按总方差归一化，而

$$
C_{xy}(f)
=\frac{S_{xy}(f)}{\sqrt{S_{xx}(f)S_{yy}(f)}}
$$

按每个频率上的功率归一化。相干性可以理解为将总的线性相关拆分到各个频率后分别考察。

但是，相干性并不等于互相关的 Fourier 变换本身。互相关的 Fourier 变换是互功率谱 $S_{xy}$；只有再除以两个自功率谱的几何平均，才得到复相干度或 MSC。

#### 高相关但局部相干结构不同

两个宽带信号可能具有较高的零延迟相关系数，但这种关系可能主要由少数强能量频带贡献。相干性可以显示究竟哪些频率相关、哪些频率被噪声淹没。

反过来，一个窄频带上的相干性可能接近 1，但该频带占总能量很小，因此整体时域相关系数仍可能不高。

#### 固定时延对两者的影响

固定时延会使零延迟相关降低，但互相关函数在正确延迟处仍可出现高峰。频域中，理想固定时延不会降低理论相干性，只会引入线性相位：

$$
|C_{xy}(f)|=1,
\qquad
\angle C_{xy}(f)=\pm2\pi f\tau_0.
$$

因此，如果两个信号形状相同但存在时延：

– 只看零延迟 Pearson 相关可能误判为关系较弱；
– 看互相关峰可以估计时延；
– 看相干性和互谱相位可以识别相关频带和频率相关的相位关系。

#### 非线性关系对两者的影响

普通 Pearson 相关和 MSC 都主要刻画线性二阶关系。若

$$
y(t)=x^2(t),
$$

即使 $y$ 完全由 $x$ 决定，$x$ 与 $y$ 的普通相关和同频 MSC 也可能很低。非线性会把能量转移到谐波和互调频率，这种跨频关系可用双相干性（bicoherence）、高阶谱或互信息等方法研究。

#### 选择指标的建议

– 要找相对时延：使用互相关、GCC 或互谱相位；
– 要比较整体线性共同变化：使用 Pearson 相关系数；
– 要判断哪些频率存在稳定线性联系：使用 MSC；
– 要保留幅相关系：同时查看互功率谱、复相干度和传递函数；
– 要分析时间变化的频率关系：使用短时相干性或小波相干性；
– 要排除共同第三变量：使用偏相关或偏相干性；
– 要分析非线性依赖：考虑互信息、距离相关、高阶谱或非线性系统辨识。

—
### 相干性上界 $\gamma_{xy}^2(f)\le1$ 的推导

MSC 满足 $0\le\gamma_{xy}^2(f)\le1$ 有两种等价的推导路径，值得记住。

**路径 A（Cauchy–Schwarz）**：对任意固定频率 $f$，在带通滤波极窄化的意义下把 $x$、$y$ 分别视为该频率的复解析成分 $\tilde X(f)$、$\tilde Y(f)$。互谱与自谱可解释为

$$
S_{xy}(f)=\mathbb E[\tilde X(f)\tilde Y^*(f)],
\qquad
S_{xx}(f)=\mathbb E[|\tilde X(f)|^2],
\qquad
S_{yy}(f)=\mathbb E[|\tilde Y(f)|^2].
$$

对二阶矩内积应用 Cauchy–Schwarz：

$$
|\mathbb E[\tilde X\tilde Y^*]|^2
\le\mathbb E[|\tilde X|^2]\mathbb E[|\tilde Y|^2],
$$

两边除以右端得 $\gamma_{xy}^2(f)\le1$，等号仅在 $\tilde Y(f)=c(f)\tilde X(f)$ 几乎必然成立时取到，即两个复窄带信号相差一个确定复常数（幅度和相位）。

**路径 B（谱矩阵半正定）**：把 $2\times2$ 谱矩阵

$$
\boldsymbol S(f)=
\begin{bmatrix}
S_{xx}(f)&S_{xy}(f)\\S_{xy}^*(f)&S_{yy}(f)
\end{bmatrix}
$$

代入 $\boldsymbol S(f)\succeq 0$，即所有主子式非负。$1\times1$ 主子式非负给出 $S_{xx}(f)\ge 0$、$S_{yy}(f)\ge 0$；$2\times2$ 行列式非负给出

$$
\det\boldsymbol S(f)=S_{xx}(f)S_{yy}(f)-|S_{xy}(f)|^2\ge0,
$$

再次得到 $\gamma_{xy}^2(f)\le1$。这一路径可自然推广到多通道的**多重相干性**和**偏相干性**：都对应把谱矩阵按 Schur 补分块后仍保持半正定的结论。

两种路径也告诉我们：$\gamma_{xy}^2(f)=1$ 意味着在该频率上 $\tilde Y(f)=H(f)\tilde X(f)$ 几乎必然成立，$H(f)$ 是**复系数**，可以任意改变幅度和相位。因此高相干并不要求两个信号在时域一致，也不要求存在因果关系。

—
### LTI 噪声模型下的相干性—信噪比公式

第七章正文给出了

$$
\gamma_{xy}^2(f)=\frac{\operatorname{SNR}(f)}{1+\operatorname{SNR}(f)}
$$

这一结果，本节把它拆解成不同噪声位置下的通用形式，方便实际测量。仍设输入 $x$、输出 $y$、噪声 $n$，$x\perp n$，且系统为 LTI。

**情形 1：噪声只加在输出。** $y=h*x+n_{\mathrm o}$。定义输出信噪比

$$
\operatorname{SNR}_{\mathrm o}(f)
=\frac{|H(f)|^2 S_{xx}(f)}{S_{n_\mathrm o n_\mathrm o}(f)},
$$

则

$$
\gamma_{xy}^2(f)
=\frac{\operatorname{SNR}_\mathrm o(f)}
{1+\operatorname{SNR}_\mathrm o(f)}.
$$

这是最常用的解释，且当输出全部由输入线性响应产生时相干性为 1。

**情形 2：噪声只加在输入观测。** 真实输入 $x$，观测到 $\tilde x=x+n_\mathrm i$；输出 $y=h*x$。此时

$$
S_{\tilde x y}(f)=H^*(f)S_{xx}(f),
\qquad
S_{\tilde x\tilde x}(f)=S_{xx}(f)+S_{n_\mathrm i n_\mathrm i}(f).
$$

对应相干性

$$
\gamma_{\tilde xy}^2(f)
=\frac{S_{xx}(f)}
{S_{xx}(f)+S_{n_\mathrm i n_\mathrm i}(f)}
=\frac{\operatorname{SNR}_\mathrm i(f)}
{1+\operatorname{SNR}_\mathrm i(f)}.
$$

**情形 3：输入、输出观测都有噪声。**

$$
\gamma_{\tilde x\tilde y}^2(f)
=\frac{\operatorname{SNR}_\mathrm i\operatorname{SNR}_\mathrm o}
{(1+\operatorname{SNR}_\mathrm i)(1+\operatorname{SNR}_\mathrm o)}.
$$

因此双端噪声会把可达相干性上限压低，即使真实系统是完全线性的。

从这几个公式可以读出两点工程含义：

1. 观察到的相干性缺陷可以来自**任一端**的噪声或未建模非线性，不能单独归咎于系统；
2. 若要用相干性推 SNR，应先确认噪声主要在哪一端，再选择 $H_1/H_2/H_v$ 三种 FRF 估计器之一，避免相干性—SNR 关系被系统性地误解释。

—
### 固定时延下的相位斜率与群时延限制

若 $y(t)=x(t-\tau_0)$（噪声为零的理想固定时延），则 $S_{xy}(f)=S_{xx}(f)e^{-j2\pi f\tau_0}$（在 $S_{xy}=\mathbb E[X_TY_T^*]/T$ 约定下），互谱相位

$$
\phi_{xy}(f)=-2\pi f\tau_0
$$

关于频率呈严格线性。据此可由相位斜率反推时延：

$$
\tau_0=-\frac{1}{2\pi}\frac{d\phi_{xy}}{df}.
$$

**限制一：主值裹绕。** 相位一般以 $(-\pi,\pi]$ 或 $[0,2\pi)$ 主值给出。当 $|2\pi f\tau_0|$ 超过 $\pi$ 时，相位会跨过 $\pm\pi$ 的边界发生跳变。有效做法是先做 phase unwrapping，再取斜率；否则会得到接近 0 的错误“群时延”。

**限制二：奈奎斯特相位模糊。** 数字系统只能分辨 $|\phi|\le\pi$ 的相位，因此单频观察到的时延存在整周期模糊：

$$
\tau_0\in\left\{-\frac{\phi+2\pi k}{2\pi f}\right\}_{k\in\mathbb Z}.
$$

要消除模糊，需要利用多个频率上的相位联合估计，或在时域用互相关峰做粗定位。

**限制三：宽带极限。** 数字宽带系统能可靠估计的最大时延受奈奎斯特速率约束。若相位在相邻 FFT 频点间跳变超过 $2\pi$（即 $2\pi\Delta f\,\tau_0>2\pi$，$\tau_0>1/\Delta f$），unwrapping 也无法恢复原始斜率。此时应先在时域使用 GCC 定位一个粗糙延迟并对齐，再用高相干频带的相位做精细估计。

**限制四：低相干频段不可信。** 群时延

$$
\tau_g(f)=-\frac{1}{2\pi}\frac{d\phi_{xy}(f)}{df}
$$

对相位的数值微分十分敏感。在 $\gamma_{xy}^2(f)$ 较低的频段，相位估计方差按

$$
\operatorname{var}\{\widehat\phi_{xy}(f)\}
\approx\frac{1-\gamma_{xy}^2(f)}{2K\gamma_{xy}^2(f)}
$$

放大（$K$ 为等效平均次数），导致群时延曲线出现巨大噪刺。工程报告应在群时延图上叠加相干性阈值，只解释 $\gamma_{xy}^2\ge\gamma^2_{\mathrm{crit}}$ 频段内的群时延。

—

—

## 八、谱与相干性的估计方法

理论谱密度是无限数据或统计期望下的量，工程中必须从有限数据估计。相干性是谱估计的比值，对分段、窗函数、平均次数、同步和预处理十分敏感。

### Welch 方法的基本步骤

给定同步采样的 $x[n]$ 和 $y[n]$：

1. 去除无效样本，保证时间轴严格对齐；
2. 根据需要去均值、去趋势；
3. 将记录划分为 $K$ 个长度为 $L$ 的数据段；
4. 相邻数据段可以重叠，例如 50%；
5. 每段乘相同窗函数 $w[n]$；
6. 分别计算每段 FFT：$X_r[k]$、$Y_r[k]$；
7. 形成每段自谱和互谱；
8. 对 $K$ 段求平均；
9. 用平均后的谱计算相干性。

采用本文互谱约定：

$$
\widehat S_{xy}[k]
=\frac{1}{K}
\sum_{r=1}^{K}X_r[k]Y_r^*[k],
$$

$$
\widehat S_{xx}[k]
=\frac{1}{K}
\sum_{r=1}^{K}|X_r[k]|^2,
$$

$$
\widehat S_{yy}[k]
=\frac{1}{K}
\sum_{r=1}^{K}|Y_r[k]|^2.
$$

最后：

$$
\boxed{
\widehat\gamma_{xy}^2[k]
=\frac{|\widehat S_{xy}[k]|^2}
{\widehat S_{xx}[k]\widehat S_{yy}[k]}
}.
$$

必须先对谱平均，再做比值；不能先计算每段相干性再简单平均。

### 为什么不能只用一个数据段

若只有一个未经平滑的 FFT 段：

$$
\frac{|X[k]Y^*[k]|^2}
{|X[k]|^2|Y[k]|^2}=1
$$

只要分母非零，所有频点都会得到 1。这只是代数恒等式，不代表真实信号完全相干。

相干性估计需要谱平均或频率平滑，使跨段不稳定的相位关系在互谱中抵消，而稳定关系得到保留。

### 数据段长度与平均次数的权衡

段长 $L$ 决定近似频率分辨率：

$$
\Delta f=\frac{f_s}{L},
$$

其中 $f_s$ 是采样率。

– 段越长：频率分辨率越高，但可平均段数越少，估计方差越大；
– 段越短：平均次数更多，估计更平滑，但相邻频率成分可能无法分开。

应根据系统模态带宽、信号持续时间和所需统计置信度选择段长，而不是固定套用某个值。

### 窗函数与谱泄漏

有限截断会导致谱泄漏。Hann 窗常用于一般谱分析，因为它在主瓣宽度和旁瓣抑制之间较均衡。输入和输出必须采用同一窗和一致的分段方式。

窗函数会使相邻频点相关，并改变有效自由度；重叠段也不是完全独立的。因此，不能把实际段数 $K$ 直接视为严格独立样本数。

### 显著性阈值和置信区间

在“真实相干性为零”且有 $K$ 个独立等权谱平均的理想近似下，显著性水平 $\alpha$ 对应的 MSC 阈值常写为

$$
\boxed{
\gamma_{\mathrm{crit}}^2
=1-\alpha^{1/(K-1)}
}.
$$

例如选择 $\alpha=0.05$，估计值高于该阈值才被视为显著偏离零相干。

但实际 Welch 重叠、窗函数、频率平滑会改变有效自由度，公式中的 $K$ 应由有效独立平均次数或等效自由度替代。更严谨时可使用软件给出的置信区间、jackknife、bootstrap 或基于 surrogate data 的显著性检验。

### 常见预处理

计算相关性或相干性前，应检查：

– 两通道是否使用相同采样率和时间戳；
– 是否存在恒定或变化的采样延迟；
– 是否去除均值、线性趋势和明显漂移；
– 是否存在丢样、插值和时钟漂移；
– 抗混叠滤波是否一致；
– 传感器是否饱和或量化不足；
– 是否需要对工频、机械转频等已知干扰单独处理；
– 数据是否近似平稳，还是应分工况分析。

不要为了获得“更漂亮”的相干曲线而任意滤波。任何滤波都应同时考虑对幅度、相位、自由度和统计解释的影响。

### 相干性图的正确阅读方式

一张完整的相干性分析图通常应联合展示：

1. $S_{xx}(f)$ 和 $S_{yy}(f)$：确认该频率是否有足够能量；
2. $|S_{xy}(f)|$：观察共同谱成分；
3. $\gamma_{xy}^2(f)$：观察归一化线性关联强度；
4. $\angle S_{xy}(f)$：只在高相干区解释相位；
5. 显著性阈值或置信区间；
6. 采样率、段长、窗、重叠率和平均次数。

在两个自谱都接近噪声底时，即使相干性偶然偏高，也应谨慎解释。相干性是归一化比值，不能替代绝对功率分析。

### 相干性分析的报告模板

工程报告至少应写清：

– 数据时长与采样率；
– 输入和输出通道的单位及校准信息；
– 去趋势、滤波和异常值处理；
– 窗函数、段长、重叠率、FFT 长度；
– 单边谱还是双边谱；
– 谱平均方法和有效自由度；
– 相干性定义及互谱共轭约定；
– 显著性阈值或置信区间；
– 哪些频带可解释，以及哪些频带因低功率或低相干不可解释。

—
### 谱估计方法与工程实践

#### Periodogram

有限记录 $x[n]$ 的 periodogram 可写成

$$
\widehat S_{xx}^{\mathrm{per}}(f)
=\frac{1}{Nf_s}\left|
\sum_{n=0}^{N-1}x[n]e^{-j2\pi fn/f_s}
\right|^2.
$$

这里采用的归一化使谱对频率积分近似等于信号的平均平方值。若使用窗函数 $w[n]$：

$$
\widehat S_{xx}(f)
=\frac{1}{f_s U}
\left|\sum_{n=0}^{N-1}w[n]x[n]e^{-j2\pi fn/f_s}\right|^2,
$$

其中

$$
U=\frac1N\sum_{n=0}^{N-1}|w[n]|^2.
$$

不同软件可能把 $N$、$f_s$、窗功率和 FFT 长度放在不同位置，因此应以“积分是否恢复平均功率”为归一化检查。

#### Bartlett 与 Welch 方法

Bartlett 方法将数据分成不重叠的若干段，分别计算 periodogram 后平均。Welch 方法进一步允许：

– 段间重叠；
– 每段使用窗函数；
– 对窗功率做归一化。

平均可以降低谱估计方差，但会牺牲部分频率分辨率。Welch 方法同样适用于互功率谱：

$$
\widehat S_{xy}(f)
=\frac{1}{K}\sum_{r=1}^{K}
X_r(f)Y_r^*(f).
$$

#### 自回归谱估计

参数化方法假设数据可由 AR 模型描述：

$$
x[n]+\sum_{k=1}^{p}a_kx[n-k]=u[n],
$$

其中 $u[n]$ 为白噪声。若其方差为 $\sigma_u^2$，则 PSD 形式为

$$
S_{xx}(e^{j\omega})
=\frac{\sigma_u^2}
{|1+\sum_{k=1}^{p}a_ke^{-j\omega k}|^2}.
$$

AR 谱可以在短数据下提供较高的频率分辨率，但模型阶数不合适时会产生伪峰或过度平滑。它不应与非参数 Welch 谱无条件混用。

#### ESD、PSD 与 periodogram 的归一化检查

对有限记录，需要先明确目标：

– 若要估计这段记录携带的总能量，使用 ESD，并使频率积分得到能量；
– 若要估计长期平均功率或随机过程 PSD，使用按记录时长、采样率和窗能量归一化的 PSD；
– 若只是比较相对峰值，可以使用未绝对校准的 periodogram，但不能把纵轴直接解释为物理功率密度。

离散频率间隔为

$$
\Delta f=\frac{f_s}{N_{\mathrm{FFT}}}.
$$

若 PSD 的单位为 $\mathrm{V}^2/\mathrm{Hz}$，频带功率应近似为

$$
P_{[f_1,f_2]}
\approx\sum_{f_k\in[f_1,f_2]}
\widehat S_{xx}(f_k)\Delta f.
$$

少乘一个 $\Delta f$ 会把“谱密度”误当成“频点功率”。

#### 单边谱与双边谱

对于实信号，双边谱在正负频率上对称。单边 PSD 只保留 $0\le f\le f_s/2$，并将非 DC、非 Nyquist 频点的功率乘 2：

$$
\widehat S_{\mathrm{one-sided}}(f)
=2\widehat S_{\mathrm{two-sided}}(f).
$$

这一操作保证单边积分与双边积分相同。对复信号或解析信号，不能默认乘 2，因为其负频率不一定是正频率的重复信息。

#### 频率分辨率、窗主瓣与泄漏

FFT 频点间隔

$$
\Delta f=\frac{f_s}{N_{\mathrm{FFT}}}
$$

不等于真正的可分辨频率间隔。真正分辨能力还取决于有效观测时长和窗函数主瓣宽度。

– 增大零填充：让曲线采样更密，但不提升真实分辨率；
– 增大有效观测时长：通常提升频率分辨率；
– 降低窗旁瓣：减少泄漏，但通常增宽主瓣；
– 两个接近频率的峰能否分开，取决于窗和数据长度，不仅取决于 FFT 点数。

#### 互谱估计的对齐要求

计算互功率谱前应保证：

1. 两个通道来自同一时间区间；
2. 采样率和采样时刻一致；
3. 分段边界、窗函数和 FFT 长度一致；
4. 通道延迟没有被未知缓存或硬件滤波改变；
5. 未在一个通道上使用不同的去趋势方式；
6. 频率轴和单边/双边约定相同。

幅度看起来合理但通道错位，可能导致互谱相位错误、相干性降低和传递函数偏差。

#### 频带积分和带内功率

对于 PSD，频带 $[f_1,f_2]$ 内的平均功率为

$$
P_{[f_1,f_2]}
=\int_{f_1}^{f_2}S_{xx}(f)df.
$$

离散估计中使用频率间隔加权：

$$
\widehat P_{[f_1,f_2]}
=\sum_{k:f_k\in[f_1,f_2]}
\widehat S_{xx}(f_k)\Delta f.
$$

幅度谱的峰值不能直接当作带内功率；带内功率来自平方量的积分。

#### dB 表示

功率谱常使用

$$
S_{\mathrm{dB}}(f)=10\log_{10}\frac{S_{xx}(f)}{S_{\mathrm{ref}}(f)}.
$$

幅度谱常使用

$$
A_{\mathrm{dB}}(f)=20\log_{10}\frac{|X(f)|}{A_{\mathrm{ref}}}.
$$

因为功率与幅度平方成正比，所以 $10\log_{10}$ 功率和 $20\log_{10}$ 幅度在参考量匹配时等价。对 PSD 使用 dB/Hz 时，要明确参考 PSD；对带内积分，不能先简单平均 dB 值再当作线性功率。

#### 常见谱估计故障

– 把 $|X|$ 当作 PSD，漏掉平方和归一化；
– 将有限记录 ESD 与长期 PSD 混称；
– 忘记乘频率间隔，导致积分功率错误；
– 单边谱没有正确折叠负频率；
– 误以为零填充增加了能量或分辨率；
– 使用单段 FFT 估计相干性，得到虚假的 1；
– 互谱通道顺序和理论公式相反，导致传递函数相位取反；
– 忽略窗函数的功率归一化；
– 在低自谱功率频带过度解释互谱相位和相干性；
– 对非平稳信号使用整段单一 PSD，掩盖工况变化。

—
### Periodogram 的偏差与方差性质

单段 periodogram 是最基础的 PSD 估计，但其统计性质并不理想。设 $x[n]$ 为零均值宽平稳过程、真实谱为 $S_{xx}(f)$，加窗 periodogram

$$
\widehat S_{xx}^{\mathrm{per}}(f)
=\frac{1}{f_sU N}\left|\sum_{n=0}^{N-1}w[n]x[n]e^{-j2\pi fn/f_s}\right|^2
$$

有如下常用近似性质：

– **偏差**：$\mathbb E[\widehat S_{xx}^{\mathrm{per}}(f)]=(S_{xx}*W)(f)$，其中 $W$ 是窗的功率谱。窗主瓣宽度决定了尖峰的展宽，旁瓣决定了远离峰值处的泄漏偏差。$N\to\infty$ 且窗满足一定条件时估计渐近无偏。
– **方差**：对连续谱区段的高斯过程，

$$
\operatorname{var}\{\widehat S_{xx}^{\mathrm{per}}(f)\}\approx S_{xx}^2(f),
$$

与 $N$ **无关**。因此增加数据长度可以改善分辨率，但**不能**单靠加长记录降低方差。

– **相邻频率相关性**：单段 periodogram 的 $\widehat S_{xx}(f_k)$ 在不同 FFT 频点上近似不相关（矩形窗、白噪声）。这是 Bartlett/Welch 之所以能通过频段平滑降低方差的基础。

后果是：单段 periodogram 在真谱上下浮动约 $\pm100\%$，看上去“毛刺很密”，容易被误读为许多小峰。降低方差的唯一手段是**平均**（分段时间平均或频域邻域平均），这正是下一节 Welch 方法的动机。

一致性方面，Welch 方法在 $K\to\infty$、每段长 $L$ 固定的极限下方差按 $\sim S^2/K$ 衰减，但主瓣宽度由 $L$ 决定不再变化；只有再让 $L\to\infty$（同时保证 $K/L\to0$ 一类条件），Welch 估计才是既无偏又一致的。工程上无法同时把 $L$ 和 $K$ 推到无穷，必须在**分辨率**与**方差**之间权衡。

—
### Welch 参数选择的实用准则

Welch 方法有四个主要旋钮：段长 $L$、重叠率 $\alpha$、窗函数 $w[n]$、FFT 长度 $N_{\mathrm{FFT}}$。选择时可按以下准则：

**段长 $L$**：由目标频率分辨率反推。若要求 3 dB 分辨率 $\Delta f_{3\mathrm{dB}}$、选择的窗主瓣系数为 $\beta_w$（Hann 约 1.5、Hamming 约 1.36、Blackman–Harris 约 2.0）：

$$
L\ge\frac{\beta_w f_s}{\Delta f_{3\mathrm{dB}}}.
$$

窄带精细谱应取大 $L$，宽带 PSD 或系统辨识可取更小的 $L$ 换取更多平均。

**重叠率 $\alpha$**：Hann 窗常用 50%，Hamming 常用 50%，Blackman–Harris 常用 75%。重叠可提高对边缘样本的利用率并减小方差，但不能视为增加独立样本数。等效独立段数近似为

$$
K_{\mathrm{eff}}
\approx\frac{K}{1+2\sum_{r=1}^{R-1}\rho_w^2(r)},
$$

其中 $\rho_w(r)$ 是窗在偏移 $r$ 段长的归一化重叠系数；Hann 50% 重叠时 $K_{\mathrm{eff}}\approx K/1.34$。

**窗函数 $w[n]$**：Hann 主瓣较宽、旁瓣衰减快，常用于一般谱分析；Hamming 主瓣略窄但旁瓣不衰减更快，适合谐波密集场合；Blackman–Harris/Flattop 适合幅度校准；矩形窗几乎只在正弦频率对齐 bin 时使用。**输入、输出通道必须使用相同的窗**，否则互谱幅相都会被引入系统误差。

**FFT 长度 $N_{\mathrm{FFT}}$**：$\ge L$，通常取 $L$ 的 2 的幂零填充倍数。零填充只让曲线视觉上更密，不提升可分辨率。

**总记录时长**：由所需自由度反推。若要求等效自由度 $\nu$（Welch 每段自由度约为 2）：

$$
N_{\mathrm{total}}
\ge\frac{\nu L(1-\alpha)^{-1}}{2},
$$

再对照实际能获取的记录时长确认是否可行。

上述准则应在具体系统上做小规模试算：先用短记录扫一遍 Welch，检查主要峰值位置和宽度是否合理，再定最终参数。

—
### 先平均谱、再求相干：一个反例

Welch 方法必须**先对多段互谱和自谱求平均，再做比值**：

$$
\widehat\gamma_{xy}^2[k]
=\frac{|\overline{X_rY_r^*}|^2}
{\overline{|X_r|^2}\cdot\overline{|Y_r|^2}}.
$$

若换成“先算每段的相干性再平均”，即

$$
\overline{\gamma_r^2}[k]
=\frac{1}{K}\sum_{r=1}^{K}
\frac{|X_r[k]Y_r^*[k]|^2}{|X_r[k]|^2|Y_r[k]|^2},
$$

分子分母同段严格相消（第八章开头已指出），每段的比值恒为 1，其算术平均仍为 1。任何两个非零信号无论是否真实相关都会得到 $\overline{\gamma_r^2}\equiv1$。

反例：取 $x_r[n]$ 与 $y_r[n]$ 都是独立高斯白噪声，段间彼此独立。真实 $\gamma_{xy}^2(f)\equiv0$。

– **正确 Welch**：$K$ 段互谱在不同段之间随机相位平均后趋于 0，$|\overline{X_rY_r^*}|^2$ 按 $1/K$ 缩小，$\widehat\gamma_{xy}^2\to0$；
– **错误顺序**：每段代数上得 1，平均后仍是 1，与真值完全相反。

因此“先谱平均，后作比值”不仅是数值细节，还是相干性估计**能否收敛到真值**的结构性要求。若使用了预平均等价物（例如频域邻域平滑或多锥形谱），也要在**同一个平滑核**下对分子和分母分别平滑，再做比值，不能在比值层面平滑。

—
### MSC 估计的分布与有效自由度

对独立平均 $K$ 段、真实相干为 $\gamma^2$ 的近似高斯谱假设下，Welch MSC 估计的均值近似满足

$$
\mathbb E[\widehat\gamma_{xy}^2]
\approx\gamma^2+\frac{1-\gamma^2}{K}
$$

（存在一个 $1/K$ 阶的**正偏差**，尤其在真值接近 0 时明显），方差近似为

$$
\operatorname{var}\{\widehat\gamma_{xy}^2\}
\approx\frac{2\gamma^2(1-\gamma^2)^2}{K}.
$$

这些近似依赖以下**有效自由度条件**：

1. 各段近似独立：Hann 50% 重叠时 $K_{\mathrm{eff}}\approx K/1.34$ 而非 $K$；
2. 每段窗内数据近似平稳；
3. 谱在窗主瓣内近似平坦（否则窗泄漏会导致低相干区被高相干区“污染”）；
4. 各段互谱和自谱使用同一窗、同一分段方式；
5. 高斯近似成立，即感兴趣频率上的谱系数近似 Gauss，且不受强线谱主导。

在“真实相干为零”假设下，$K$ 独立平均的显著性阈值

$$
\gamma^2_{\mathrm{crit}}=1-\alpha^{1/(K-1)}
$$

也应把 $K$ 替换为 $K_{\mathrm{eff}}$，否则会得到过于宽松的阈值，把噪声起伏当作显著。

置信区间可通过 Fisher-Z 型变换或 jackknife/bootstrap 取得。对方差非常敏感的场合（例如生物医学、结构健康监测），建议在结果图上把有效自由度和阈值线一并画出，让读者知道判读的置信程度。

—
### scipy 参考代码示例

以下 Python 例子演示如何用 `scipy.signal` 一致地计算 PSD、CSD 和 MSC。它可以作为工程实现的起点，但生产环境应显式写清窗、重叠、单双边约定和归一化。

“`python
import numpy as np
from scipy import signal

# 假设 x, y 为等长的一维实信号，采样率 fs
fs = 1000.0 # Hz
nperseg = 1024 # 段长 L
noverlap = nperseg // 2 # 50% 重叠
window = ‘hann’ # Hann 窗
nfft = nperseg # 与段长相同（不额外零填充）

# 自谱、互谱、相干性用同一套参数
f, Pxx = signal.welch(x, fs=fs, window=window,
nperseg=nperseg, noverlap=noverlap,
nfft=nfft, detrend=’constant’,
return_onesided=True, scaling=’density’)
_, Pyy = signal.welch(y, fs=fs, window=window,
nperseg=nperseg, noverlap=noverlap,
nfft=nfft, detrend=’constant’,
return_onesided=True, scaling=’density’)
_, Pxy = signal.csd(x, y, fs=fs, window=window,
nperseg=nperseg, noverlap=noverlap,
nfft=nfft, detrend=’constant’,
return_onesided=True, scaling=’density’)

# MSC：注意与 signal.coherence 的定义一致
Cxy = np.abs(Pxy) ** 2 / (Pxx * Pyy)

# 也可直接调用高层 API 进行交叉检查
_, Cxy_ref = signal.coherence(x, y, fs=fs, window=window,
nperseg=nperseg, noverlap=noverlap,
nfft=nfft, detrend=’constant’)

# 相干性显著性阈值（假设真实相干为 0）
K = (len(x) – noverlap) // (nperseg – noverlap) # 段数
alpha = 0.05
gamma_crit = 1 – alpha ** (1 / (K – 1))
“`

要点提示：

– `scaling=’density’` 给出 PSD（单位 $\mathrm{V}^2/\mathrm{Hz}$）；换成 `’spectrum’` 得到分段功率谱（单位 $\mathrm{V}^2$）；两者积分/求和的物理意义不同。
– `return_onesided=True` 只对实信号有效；复信号必须传 `False`。
– `signal.csd` 与 `signal.welch` 使用**相同的分段和窗**才能保证 `Pxy/sqrt(Pxx*Pyy)` 与 `signal.coherence` 结果一致。
– 有效独立段数 $K_{\mathrm{eff}}$ 小于名义 $K$；显著性阈值 `gamma_crit` 上式偏乐观，需要按窗自相关修正。
– 若做因果分析或群时延估计，应从 `Pxy` 复数值取相位并进行 `np.unwrap`，并只在 `Cxy >= gamma_crit` 频段解释相位。

—

—

## 九、通信、检测与系统辨识应用

### 模板匹配

已知模板 $s[n]$，观测信号为 $r[n]$。通过计算

$$
C[m]=\sum_n r[n]s^*[n-m]
$$

可以在 $r[n]$ 中寻找模板出现的位置。最大峰值的位置给出模板的估计延迟。

### 匹配滤波器

匹配滤波器的冲激响应通常与模板的共轭时间反转有关：

$$
h[n]=s^*[N-1-n].
$$

将接收信号与该滤波器卷积，等价于计算接收信号与模板的相关。匹配滤波器在加性白噪声下能够最大化某一采样时刻的输出信噪比。

#### AWGN 下匹配滤波的最优 SNR 推导

设已知模板 $s(t)$（能量 $E_s=\int|s(t)|^2dt$），接收信号为

$$
r(t)=s(t)+w(t),
$$

其中 $w(t)$ 是双边 PSD 为 $N_0/2$ 的加性白 Gaussian 噪声。用任意线性滤波器 $h(t)$ 对 $r$ 滤波，令输出

$$
y(t)=(h*r)(t)=y_s(t)+y_w(t).
$$

固定采样时刻 $t_0$，输出信号分量的瞬时”功率”为 $|y_s(t_0)|^2$，输出噪声方差为

$$
\sigma_w^2=\frac{N_0}{2}\int|h(\tau)|^2d\tau.
$$

采样时刻的输出 SNR 定义为

$$
\operatorname{SNR}_{\mathrm{out}}
=\frac{|y_s(t_0)|^2}{\sigma_w^2}
=\frac{\left|\int h(\tau)s(t_0-\tau)d\tau\right|^2}
{\frac{N_0}{2}\int|h(\tau)|^2d\tau}.
$$

由 Cauchy–Schwarz 不等式：

$$
\left|\int h(\tau)s(t_0-\tau)d\tau\right|^2
\le\left(\int|h(\tau)|^2d\tau\right)
\left(\int|s(t_0-\tau)|^2d\tau\right)
=E_s\int|h(\tau)|^2d\tau,
$$

等号当且仅当

$$
h(\tau)=k\,s^*(t_0-\tau)
$$

时成立（$k$ 为常数）。代入得

$$
\boxed{\operatorname{SNR}_{\mathrm{out},\max}=\frac{2E_s}{N_0}}.
$$

结论：最优线性滤波器就是模板的共轭时间反转（离散域即 $h[n]=s^*[N-1-n]$），最优 SNR 只依赖模板能量与噪声 PSD，与波形形状无关。这也解释了为什么雷达和通信中通过增加脉冲能量（更宽脉冲、脉冲压缩、编码增益）来改善检测性能，而不必改变载频或匹配滤波结构本身。

### 同步与时延估计

发送信号 $s(t)$，接收信号近似为

$$
r(t)=a s(t-\tau_0)+w(t),
$$

其中 $a$ 是幅度或复衰落系数，$w(t)$ 是噪声。计算

$$
R_{rs}(\tau)=\int r(t)s^*(t-\tau)dt
$$

时，在 $\tau\approx\tau_0$ 附近通常出现峰值。

#### 广义互相关与 GCC-PHAT

在多传感器或未知信号场景（如麦克风阵列、被动声纳、TDOA 定位），两路观测常写成

$$
x_1(t)=s(t)+n_1(t),
\qquad
x_2(t)=\alpha\,s(t-\tau_0)+n_2(t).
$$

广义互相关（Generalized Cross-Correlation, GCC）通过在频域引入权函数 $\Psi(f)$ 增强峰值锐度：

$$
R^{\mathrm{GCC}}_{x_1x_2}(\tau)
=\int_{-\infty}^{\infty}\Psi(f)\,S_{x_1x_2}(f)\,e^{j2\pi f\tau}\,df.
$$

不同权函数对应不同方法：$\Psi\equiv1$ 就是普通互相关；Roth 加权 $\Psi=1/S_{x_1x_1}$；SCOT 加权 $\Psi=1/\sqrt{S_{x_1x_1}S_{x_2x_2}}$。工程中最常用的是相位变换（Phase Transform, PHAT）加权：

$$
\boxed{\Psi_{\mathrm{PHAT}}(f)=\frac{1}{|S_{x_1x_2}(f)|}},
$$

代入后

$$
R^{\mathrm{PHAT}}_{x_1x_2}(\tau)
=\int\frac{S_{x_1x_2}(f)}{|S_{x_1x_2}(f)|}e^{j2\pi f\tau}df
=\int e^{j\angle S_{x_1x_2}(f)}e^{j2\pi f\tau}df.
$$

也就是只保留互谱的相位、幅度归一化为 1。理想窄带无噪时 $\angle S_{x_1x_2}(f)=-2\pi f\tau_0$，积分给出

$$
R^{\mathrm{PHAT}}_{x_1x_2}(\tau)\propto\delta(\tau-\tau_0),
$$

峰位与信号谱形状无关。DFT 实现：

$$
\widehat R^{\mathrm{PHAT}}[m]
=\operatorname{IDFT}\!\left\{
\frac{X_1[k]X_2^*[k]}{|X_1[k]X_2^*[k]|+\varepsilon}
\right\}.
$$

工程要点：

– 与普通互相关相比，PHAT 在混响、色噪声、宽带非平稳信号下峰值更锐、旁瓣更低；
– 但因抛弃了幅度信息，低 SNR 或稀疏谱段会因除以小数而被噪声主导，需要加正则化 $\varepsilon$ 或谱阈值；
– 亚采样精度可通过在峰值附近做抛物线插值或频域相位斜率拟合获得；
– 多对麦克风给出多组 $\tau_{ij}$，可联立求解声源坐标（TDOA 定位）。

### 系统辨识流程

第五章给出了单输入线性系统的传递函数估计

$$
\widehat H(f)=\frac{S_{yx}(f)}{S_{xx}(f)},
$$

第七章给出了相干性 $\gamma^2_{xy}(f)$。把二者组合起来，一个可复现的非参数系统辨识流程通常包括：

1. **激励设计**：选取带宽覆盖被测系统关注频段的输入 $x[n]$。常用宽带白噪声、伪随机二进制序列（PRBS）、多正弦（multisine）或调频扫频（chirp）。要求：
– $S_{xx}(f)$ 在关注频带非零；
– 单次实验时间足够长，能提供 $K$ 段独立平均；
– 幅值不激发被测系统的非线性区间。
2. **同步采集**：$x[n]$ 与 $y[n]$ 用同一时钟采样，尽量避免通道间抗混叠滤波器不一致带来的固定相位偏差；若有偏差，应在事后用已知参考做校正。
3. **预处理**：去均值、去线性趋势、可选带通滤波、必要时对异常段进行剔除或加权。
4. **分段谱估计**：使用第八章的 Welch 方法，段长 $L$、重叠 $50\%\sim75\%$、加 Hann/Hamming 窗，得到

$$
\widehat S_{xx}(f),\ \widehat S_{yy}(f),\ \widehat S_{yx}(f).
$$
5. **传递函数估计**：常用三个估计量各有偏差方向

$$
\widehat H_1=\frac{\widehat S_{yx}}{\widehat S_{xx}},\quad
\widehat H_2=\frac{\widehat S_{yy}}{\widehat S_{xy}},\quad
\widehat H_v=\sqrt{\widehat H_1\widehat H_2}.
$$
$\widehat H_1$ 对输出噪声鲁棒，$\widehat H_2$ 对输入噪声鲁棒，$\widehat H_v$ 折衷。
6. **相干性检查**：计算 $\widehat\gamma^2_{xy}(f)$，只在相干性足够高（例如 $\ge 0.8$）的频段内报告 $\widehat H(f)$；相干性显著低于 1 的频段一般对应：输入功率不足、非线性、外部干扰或存在未测量输入。
7. **模型拟合与验证（可选）**：将非参数 $\widehat H(f)$ 拟合成参数模型（如 ARX、状态空间、传递函数），并用未参与辨识的独立数据段做残差与预测检验。

这套流程也可推广到 MIMO 系统：把标量自谱换成谱矩阵，$\widehat{\mathbf H}=\widehat{\mathbf S}_{yx}\widehat{\mathbf S}_{xx}^{-1}$，并使用多重相干性（见 7.3 节）代替单频相干性。

### 扩频通信中的解扩

扩频码与接收信号相关，可以将目标码片序列的能量累积起来，而与其他码或噪声不匹配的分量趋于抵消。

### 雷达和声纳

发射已知波形，接收目标回波后与发射波形做匹配滤波或互相关：

– 峰值位置估计传播时延；
– 时延乘以传播速度可以估计距离；
– 峰值相位或多普勒结构可以用于估计运动信息。

### 自适应滤波中的相关梯度

自适应滤波把相关从“测量两个信号的相似性”变成“寻找使误差下降的参数更新方向”。共轭的位置由输出模型、误差定义和求导变量共同决定。

#### 实数 LMS：从链式法则开始

令输入回归向量为

$$
\boldsymbol{x}_n=
\begin{bmatrix}
x[n]&x[n-1]&\cdots&x[n-L+1]
\end{bmatrix}^{T},
$$

采用

$$
y[n]=\boldsymbol w^T\boldsymbol x_n,
\qquad
e[n]=d[n]-y[n],
\qquad
J[n]=\frac12e^2[n].
$$

对第 $\ell$ 个抽头：

$$
\frac{\partial J[n]}{\partial w[\ell]}
=\frac{\partial J[n]}{\partial e[n]}
\frac{\partial e[n]}{\partial w[\ell]}
=e[n]\cdot[-x[n-\ell]].
$$

因此

$$
\nabla_{\boldsymbol w}J[n]=-\boldsymbol x_ne[n],
$$

梯度下降给出

$$
\boxed{
\boldsymbol w[n+1]=\boldsymbol w[n]+\mu\boldsymbol x_ne[n]
}.
$$

若把误差定义成 $e[n]=y[n]-d[n]$，更新式的符号会相应改变。加号或减号不能脱离误差定义和代价函数单独判断。

#### 复数参数与 Wirtinger 导数

复数代价函数

$$
J[n]=\frac12e[n]e^*[n]
$$

是实值函数，但同时依赖 $w$ 和 $w^*$。令 $w=w_R+jw_I$，Wirtinger 导数定义为

$$
\frac{\partial}{\partial w}
=\frac12\left(\frac{\partial}{\partial w_R}
-j\frac{\partial}{\partial w_I}\right),
\qquad
\frac{\partial}{\partial w^*}
=\frac12\left(\frac{\partial}{\partial w_R}
+j\frac{\partial}{\partial w_I}\right).
$$

推导时把 $w$ 和 $w^*$ 暂时视为形式上独立的变量。对实值损失，常用复梯度约定为

$$
\nabla_{\boldsymbol w}J
=2\frac{\partial J}{\partial\boldsymbol w^*}.
$$

因子 2 也可以并入步长，因此不同资料中的系数可能不同，但更新方向必须一致。

#### Hermitian 模型：$y=\boldsymbol w^H\boldsymbol x$

复数自适应滤波通常采用

$$
y[n]=\boldsymbol w^H\boldsymbol x_n,
\qquad
e[n]=d[n]-\boldsymbol w^H\boldsymbol x_n,
\qquad
J[n]=\frac12|e[n]|^2.
$$

由于 $\boldsymbol w^H=(\boldsymbol w^*)^T$，有

$$
e[n]=d[n]-\boldsymbol w^H\boldsymbol x_n,
\qquad
e^*[n]=d^*[n]-\boldsymbol x_n^H\boldsymbol w.
$$

对 $\boldsymbol w^*$ 求偏导时，$e^*[n]$ 不含 $\boldsymbol w^*$：

$$
\begin{aligned}
\frac{\partial J[n]}{\partial\boldsymbol w^*}
&=\frac12\left(
e^*[n]\frac{\partial e[n]}{\partial\boldsymbol w^*}
+e[n]\frac{\partial e^*[n]}{\partial\boldsymbol w^*}
\right)\\
&=\frac12e^*[n]
\frac{\partial(d[n]-\boldsymbol w^H\boldsymbol x_n)}
{\partial\boldsymbol w^*}\\
&=-\frac12\boldsymbol x_ne^*[n].
\end{aligned}
$$

采用 $2\partial J/\partial\boldsymbol w^*$ 的复梯度约定：

$$
\boxed{
\nabla_{\boldsymbol w}J[n]=-\boldsymbol x_ne^*[n]
}.
$$

所以复数 LMS 更新为

$$
\boxed{
\boldsymbol w[n+1]=\boldsymbol w[n]+\mu\boldsymbol x_ne^*[n]
}.
$$

逐抽头为

$$
w[n+1,\ell]=w[n,\ell]+\mu x[n-\ell]e^*[n].
$$

这里的共轭来自输出使用 $\boldsymbol w^H\boldsymbol x$ 以及对 $\boldsymbol w^*$ 求 Wirtinger 导数，而不是来自某个把任意时域点乘转换成频域逐点乘法的规则。

#### $y=\boldsymbol w^T\boldsymbol x$ 参数化的对照

若系统使用

$$
y[n]=\boldsymbol w^T\boldsymbol x_n,
\qquad e[n]=d[n]-\boldsymbol w^T\boldsymbol x_n,
$$

则 $e[n]$ 依赖 $\boldsymbol w$，$e^*[n]$ 依赖 $\boldsymbol w^*$。对 $\boldsymbol w^*$ 求导得到

$$
\frac{\partial J[n]}{\partial\boldsymbol w^*}
=-\frac12\boldsymbol x_n^*e[n],
$$

从而

$$
\boldsymbol w[n+1]=\boldsymbol w[n]+\mu\boldsymbol x_n^*e[n].
$$

两种常见模型的对应关系是：

| 输出模型 | 梯度方向 | 常见更新 |
|—|—|—|
| $y=\boldsymbol w^H\boldsymbol x$ | $-\boldsymbol xe^*$ | $\boldsymbol w\leftarrow\boldsymbol w+\mu\boldsymbol xe^*$ |
| $y=\boldsymbol w^T\boldsymbol x$ | $-\boldsymbol x^*e$ | $\boldsymbol w\leftarrow\boldsymbol w+\mu\boldsymbol x^*e$ |
| 实数 $y=\boldsymbol w^T\boldsymbol x$ | $-\boldsymbol xe$ | $\boldsymbol w\leftarrow\boldsymbol w+\mu\boldsymbol xe$ |

#### 块矩阵梯度与 $X^*E$

把一个数据块写成

$$
\boldsymbol y=\boldsymbol X\boldsymbol w,
\qquad
\boldsymbol e=\boldsymbol d-\boldsymbol X\boldsymbol w,
\qquad
J=\frac12\boldsymbol e^H\boldsymbol e.
$$

若 $\boldsymbol X\in\mathbb C^{M\times L}$，则

$$
\boldsymbol w\in\mathbb C^L,
\qquad
\boldsymbol e\in\mathbb C^M,
\qquad
\boldsymbol X^H\boldsymbol e\in\mathbb C^L.
$$

对 $\boldsymbol w^*$ 求导并采用复梯度约定：

$$
\frac{\partial J}{\partial\boldsymbol w^*}
=-\frac12\boldsymbol X^H\boldsymbol e,
\qquad
\nabla_{\boldsymbol w}J=-\boldsymbol X^H\boldsymbol e.
$$

第 $\ell$ 个元素为

$$
[\boldsymbol X^H\boldsymbol e]_\ell
=\sum_nx^*[n-\ell]e[n],
$$

它正是输入和误差的相关型乘加。当卷积矩阵由 FFT 对角化时，Fourier 坐标中的逐点形式为

$$
G[k]=X^*[k]E[k],
\qquad
\boldsymbol g=\operatorname{IFFT}\{X^*[k]E[k]\}.
$$

这就是“时域抽头梯度中出现共轭乘法”的严格来源。具体的循环移位、误差补零位置和 IFFT 后有效抽头区间仍由 overlap-save、分区延迟和 MDF 数据排列决定，不能只看 $X^*E$ 三个符号。

#### 块频域 LMS 的归一化和 MDF 注意事项

常见的分区频域更新为

$$
W_p^{(m+1)}[k]
=W_p^{(m)}[k]
+\mu\frac{X_p^*[k]E[k]}
{\Phi_x[k]+\varepsilon}.
$$

其中 $\Phi_x[k]$ 是输入功率估计，$\varepsilon$ 用于避免低功率频点除零和噪声放大。工程中应确认功率估计与 FFT 缩放一致，并观察归一化后更新量是否出现异常峰值。

MDF 常用 $2N$ 点 FFT 和 $N$ 个新样本。若有效误差为 $e[0:N-1]$，一种常见排列是

$$
E_{\mathrm{time}}[0:N-1]=0,
\qquad
E_{\mathrm{time}}[N:2N-1]=e[0:N-1].
$$

放在前半段还是后半段取决于分区输入历史、处理延迟和有效抽头截取。误差块错移一个样本会在频域引入线性相位，导致更新错误的延迟分区。与 OLS/OLA 的基础边界规则一样，MDF 必须通过单位脉冲和已知延迟逐点验证。

#### 最小 NumPy 互证

下面同时验证直接线性相关、补零 FFT 相关和 Parseval 内积。相关定义为
$R_{xy}[m]=\sum_nx[n]y^*[n-m]$：

“`python
import numpy as np

rng = np.random.default_rng(7)
x = rng.normal(size=5) + 1j * rng.normal(size=5)
y = rng.normal(size=3) + 1j * rng.normal(size=3)
lags = np.arange(-(len(y) – 1), len(x))
direct = np.array([
sum(x[n] * np.conj(y[n – lag])
for n in range(len(x))
if 0 <= n - lag < len(y)) for lag in lags ]) L = len(x) + len(y) - 1 R = np.fft.ifft(np.fft.fft(x, L) * np.conj(np.fft.fft(y, L))) fft_linear = np.concatenate((R[-(len(y) - 1):], R[:len(x)])) np.testing.assert_allclose(fft_linear, direct, rtol=1e-12, atol=1e-12) X = np.fft.fft(x) y_padded = np.pad(y, (0, len(x) - len(y))) Y = np.fft.fft(y_padded) np.testing.assert_allclose(np.sum(X * np.conj(Y)) / len(x), np.vdot(y_padded, x), rtol=1e-12, atol=1e-12) print("correlation and Parseval verified") ``` 此处 `ifft` 自动包含 `1/L`；若使用未归一化的逆变换，必须手动补上该因子。 --- --- ## 十、典型例题、常见误区与检查清单 ### 例题一：两个有限序列的卷积 设 $$ x[n]=[1,2], \qquad h[n]=[3,4,5]. $$ 线性卷积长度为 $2+3-1=4$。逐项计算： $$ \begin{aligned} y[0]&=1\cdot3=3,\\ y[1]&=1\cdot4+2\cdot3=10,\\ y[2]&=1\cdot5+2\cdot4=13,\\ y[3]&=2\cdot5=10. \end{aligned} $$ 因此 $$ \boxed{x*h=[3,10,13,10]}. $$ ### 例题二：两个序列的互相关 设 $$ x[n]=[1,2], \qquad y[n]=[3,4]. $$ 取实信号相关定义 $$ R_{xy}[m]=\sum_nx[n]y[n-m]. $$ 零延迟相关为 $$ R_{xy}[0]=1\cdot3+2\cdot4=11. $$ 不同延迟下只保留重叠部分，可得到一个长度为 $2+2-1=3$ 的相关序列。采用不同的延迟索引排列时，结果可能写成 $$ [4,11,6] $$ 或其反向排列 $$ [6,11,4]. $$ 这不是数值矛盾，而是相关定义和延迟索引方向不同造成的。 ### 例题三：单位脉冲定位延迟 设模板为 $$ s[n]=\delta[n], $$ 接收信号为延迟版本 $$ r[n]=\delta[n-n_0]. $$ 使用 $$ R_{rs}[m]=\sum_nr[n]s[n-m] $$ 有 $$ \begin{aligned} R_{rs}[m] &=\sum_n\delta[n-n_0]\delta[n-m]\\ &=\delta[m-n_0]. \end{aligned} $$ 所以相关峰位于 $m=n_0$。若把相关定义写成 $s[n+m]$，峰的符号方向会随之改变。 ### 例题四：矩形脉冲的自相关 设 $$ x(t)=u(t)-u(t-T). $$ 其自相关为两个相同矩形脉冲的相关，结果是三角形： $$ R_{xx}(\tau)= \begin{cases} T-|\tau|, & |\tau|\le T,\\ 0, & |\tau|>T.
\end{cases}
$$

零延迟处

$$
R_{xx}(0)=T
$$

等于矩形脉冲的能量。

### 综合例题五：有限序列 FFT 与时域相关互证

设

$$
x[n]=[1,2,3,0],\qquad y[n]=[0,1,2,3].
$$

其中 $y[n]=x[n-1]$（右移 1 个样本）。用两条路径计算互相关并互相验证。

**路径 A：时域直接计算。** 取定义

$$
R_{xy}[m]=\sum_n x[n]y^*[n-m]=\sum_n x[n]y[n-m]
$$

（实序列）。有效 $m\in\{-3,-2,-1,0,1,2,3\}$，取重叠区求和：

$$
\begin{aligned}
R_{xy}[-3]&=x[0]y[3]=1\cdot3=3,\\
R_{xy}[-2]&=x[0]y[2]+x[1]y[3]=1\cdot2+2\cdot3=8,\\
R_{xy}[-1]&=x[0]y[1]+x[1]y[2]+x[2]y[3]=1+4+9=14,\\
R_{xy}[0]&=\textstyle\sum_n x[n]y[n]=0+2+6+0=8,\\
R_{xy}[1]&=x[1]y[0]+x[2]y[1]+x[3]y[2]=0+2+0=2,\\
R_{xy}[2]&=x[2]y[0]+x[3]y[1]=0+0=0,\\
R_{xy}[3]&=x[3]y[0]=0.
\end{aligned}
$$

峰值在 $m=-1$，对应”$y$ 相对 $x$ 延迟 $-1$”，即 $x$ 落后于 $y$ 反过来即 $y=x[n-1]$（延迟 1）：注意本约定下 $R_{xy}$ 在 $m=-\tau_0$ 处达到峰（复核见第四章“共轭方向、互谱与频域梯度”小节）。

**路径 B：FFT 频域计算。** 线性相关有效长度 $2N-1=7$，取 $L=8$ 点 FFT（$L\ge 2N-1$，方便对齐）。将 $x,y$ 各补零到长度 8：

$$
\tilde x=[1,2,3,0,0,0,0,0],\quad
\tilde y=[0,1,2,3,0,0,0,0].
$$

计算 $X[k]=\operatorname{FFT}(\tilde x)$、$Y[k]=\operatorname{FFT}(\tilde y)$，再由相关定理

$$
\widehat R_{xy}[m]=\operatorname{IFFT}\{X[k]\,Y^*[k]\}.
$$

对本例可以直接验算：$Y[k]=X[k]e^{-j2\pi k/L}$（$y$ 相对 $x$ 循环右移 1），所以

$$
X[k]Y^*[k]=|X[k]|^2 e^{j2\pi k/L}
=\operatorname{FFT}\{R_{xx}\text{ 循环左移 }1\}.
$$

IFFT 后即得 $R_{xx}$ 的循环左移；因补零足够长，循环结果与线性结果一致：峰位于 IFFT 输出的 $k=L-1=7$ 处，重排 $m\in\{-3,\dots,3\}$ 得峰在 $m=-1$，其他位置与路径 A 数值逐一对应。

**结论与互证要点。**

1. FFT 长度必须 $L\ge N_x+N_y-1$，否则前后位置会因循环回绕相互污染（对应第四章”线性 vs. 循环”警示）。
2. 频域取共轭放在参考信号上：与本约定 $R_{xy}[m]=\sum x[n]y^*[n-m]$ 保持一致的是 $X\cdot Y^*$。若写成 $X^*\cdot Y$，得到的是 $R_{yx}[m]=R_{xy}^*[-m]$，峰位方向相反。
3. IFFT 输出的自然索引从 0 开始；要读出 $m<0$ 的相关值，需要按 $m\equiv k\pmod L$ 重排。 ### 综合例题六：固定时延加噪声——互谱、相干性与群时延一体化 设宽平稳零均值实信号 $x(t)$，观测两路 $$ x_1(t)=x(t)+n_1(t),\qquad x_2(t)=x(t-\tau_0)+n_2(t), $$ 其中 $n_1,n_2$ 互不相关，也与 $x$ 不相关，PSD 分别为 $N_1(f),N_2(f)$。 **互谱。** 由 $x_2(t)$ 的频域延迟因子 $$ S_{x_1x_2}(f)=S_{xx}(f)\,e^{-j2\pi f\tau_0}, $$ 自谱为 $$ S_{x_1x_1}(f)=S_{xx}(f)+N_1(f),\qquad S_{x_2x_2}(f)=S_{xx}(f)+N_2(f). $$ **相干性。** $$ \gamma^2_{x_1x_2}(f) =\frac{|S_{xx}(f)|^2} {[S_{xx}(f)+N_1(f)][S_{xx}(f)+N_2(f)]} =\frac{1}{(1+1/\rho_1)(1+1/\rho_2)}, $$ 其中 $\rho_i(f)=S_{xx}(f)/N_i(f)$ 为该通道的谱域 SNR。特例： - 无噪声 $N_1=N_2=0$：$\gamma^2=1$，任意频率都是理想线性关系； - 一路信号被强噪声完全淹没：$\gamma^2\to0$； - 单通道白噪 $N_1=0,N_2=N_0$：$\gamma^2(f)=\rho_2/(1+\rho_2)$，与第七章"输出加不相关噪声"公式一致。 **群时延。** 互谱相位为 $$ \phi(f)=\angle S_{x_1x_2}(f)=-2\pi f\tau_0\pmod{2\pi}, $$ 因此群时延 $$ \tau_g(f)=-\frac{1}{2\pi}\frac{d\phi(f)}{df}=\tau_0 $$ 在所有 $\gamma^2(f)$ 足够大的频段应为常数 $\tau_0$。工程上按第八章“相干性分析的报告模板”执行： 1. 用 Welch 得 $\widehat S_{x_1x_1},\widehat S_{x_2x_2},\widehat S_{x_1x_2}$； 2. 计算 $\widehat\gamma^2$，只对 $\widehat\gamma^2\ge \gamma^2_{\min}$ 的 bin 保留相位； 3. 对保留的相位解缠绕后做线性回归，斜率的 $-1/(2\pi)$ 即时延估计 $\widehat\tau_0$； 4. 与 GCC-PHAT 结果互相比对：GCC-PHAT 通过对幅度归一化后 IFFT 求峰给出 $\widehat\tau_0$，与相位斜率法在高相干带内应一致；差异过大提示存在色散、非线性或多径。 本例把第五章互谱、第七章相干性和第九章 GCC-PHAT 与系统辨识用同一模型串起来：$\gamma^2$ 告诉我们哪里可信，互谱相位告诉我们时延是多少，两者缺一不可。 ### 综合例题七：白噪声经 LTI——PSD、频率响应与功率积分 设离散白噪声输入 $x[n]$，$R_{xx}[m]=\sigma_x^2\delta[m]$，故 $S_{xx}(e^{j\hat\omega})=\sigma_x^2$（双边），采样率 $f_s$。通过一阶 IIR $$ y[n]=x[n]+a\,y[n-1],\qquad |a|<1. $$ **频率响应与输出 PSD。** $$ H(e^{j\hat\omega})=\frac{1}{1-a e^{-j\hat\omega}},\qquad |H|^2=\frac{1}{1-2a\cos\hat\omega+a^2}. $$ 由第二章 LTI 结果 $S_{yy}=|H|^2 S_{xx}$： $$ S_{yy}(e^{j\hat\omega}) =\frac{\sigma_x^2}{1-2a\cos\hat\omega+a^2}. $$ **输出总功率（时域法）。** 由 $R_{yy}[m]=\sigma_x^2 a^{|m|}/(1-a^2)$ 得 $$ P_y=R_{yy}[0]=\frac{\sigma_x^2}{1-a^2}. $$ **输出总功率（频域法）。** 由 Parseval，输出功率是双边 PSD 的积分： $$ P_y=\frac{1}{2\pi}\int_{-\pi}^{\pi}\!S_{yy}(e^{j\hat\omega})\,d\hat\omega =\frac{\sigma_x^2}{2\pi}\int_{-\pi}^{\pi}\!\frac{d\hat\omega}{1-2a\cos\hat\omega+a^2} =\frac{\sigma_x^2}{1-a^2}. $$ 两条路径给出同一数值，交叉验证成功。 **辨识 $H$ 的三种谱估计量对照。** 若同步采集 $x[n],y[n]$，Welch 得 $\widehat S_{xx},\widehat S_{yy},\widehat S_{yx}$： $$ \widehat H_1=\frac{\widehat S_{yx}}{\widehat S_{xx}},\quad \widehat H_2=\frac{\widehat S_{yy}}{\widehat S_{xy}},\quad \widehat H_v=\sqrt{\widehat H_1\widehat H_2}. $$ 在本纯 LTI 无输出噪声场景下三者理论一致；若在输出加互不相关噪声 $v[n]$（PSD $N_v$），则 $\widehat H_1$ 无偏而 $\widehat H_2$ 会被高估： $$ \mathbb E[\widehat H_1]=H,\qquad \mathbb E[\widehat H_2]=H\cdot\left(1+\frac{N_v}{|H|^2 S_{xx}}\right). $$ 相干性 $$ \gamma^2(f)=\frac{|H|^2 S_{xx}}{|H|^2 S_{xx}+N_v} $$ 给出各频段 $\widehat H_1$ 的可信度，与第九章"系统辨识流程"步骤 5、6 完全对应。 **工程要点。** 白噪激励使 $S_{xx}$ 平坦，任何频段都能被辨识；换成窄带激励时，激励能量不足的频带即使数学上仍可计算 $S_{yx}/S_{xx}$，也应根据相干性把这些结果从报告中剔除，避免"数值有值 = 结果可信"的误判。 --- ### 常见误区 #### 误区：时域相乘总是对应频域相乘 错误。正确关系是： $$ \mathcal F\{x(t)h(t)\} =\frac{1}{2\pi}X(\omega)*H(\omega). $$ 时域卷积才对应频域逐点相乘。 #### 误区：相关就是卷积 不完全正确。相关可以表示成带有共轭和时间反转的卷积： $$ R_{xy}=x*y^*(-t) $$ 但相关的物理目标是匹配和相似性，卷积的典型目标是系统响应和滤波。 #### 误区：相关中共轭位置永远固定 常见形式包括： $$ X(\omega)Y^*(\omega), \qquad X^*(\omega)Y(\omega). $$ 共轭位置取决于： 1. 互相关的定义； 2. 谁作为参考信号； 3. 延迟写成 $t-\tau$ 还是 $t+\tau$； 4. Fourier 变换的正负号约定。 #### 误区：相关峰的正负方向可以脱离定义判断 不能。必须把信号模型、相关公式和延迟定义放在一起。例如“接收信号比模板晚到”在不同相关约定下可能表现为正峰或负峰。 #### 误区：FFT 频域相乘直接得到线性卷积 直接使用相同长度 FFT 后频域相乘得到的是循环卷积。若需要线性卷积，FFT 长度至少要满足 $$ L\ge N+M-1. $$ #### 误区：复数信号的卷积需要共轭乘法 错误。线性卷积的定义——无论信号是实的还是复的——始终使用普通复数乘法： $$ y[n]=\sum_k x[k]h[n-k]. $$ 共轭出现在相关和内积的定义中，不是卷积因信号变为复数后自动获得的。详见第二章“复数信号的卷积：不使用共轭乘法”。 #### 误区：自相关在所有情况下都是偶函数 只有实信号的自相关通常是偶函数。复信号满足的是共轭对称： $$ R_{xx}(-\tau)=R_{xx}^*(\tau). $$ #### 误区：相关峰越高，信号一定越相似 未归一化相关会受到信号能量影响。比较不同幅度信号时，应考虑归一化相关系数。 #### 误区：相关值只能是实数 实信号相关通常为实数；复信号互相关一般为复数。其虚部承载相对相位信息，不应随意丢弃。 #### 误区：相关只适用于确定性信号 相关同样是随机过程分析中的核心工具。随机过程的自相关函数、互相关函数和功率谱密度用于描述统计依赖、平稳性和频率能量分布。 #### 误区：相关系数为零就说明两个变量独立 一般不成立。零相关只表示不存在由协方差刻画的线性统计关联；除非变量联合 Gaussian，或者另有更强的分布条件，否则零相关不能推出统计独立。独立且二阶矩存在通常可以推出零协方差，反方向则通常不成立。 #### 误区：相干性高就证明一个信号导致另一个信号 错误。相干性反映两个信号在某个频率上的稳定线性关联，但不能单独证明因果关系。共同输入、串扰、同步时钟、反馈回路或处理链中的共同滤波都可能产生高相干性。 #### 误区：单次 FFT 就能可靠估计相干性 用同一段数据形成未经平均的周期图时，分子和分母可能发生代数抵消，使估计的 magnitude-squared coherence 退化为 1。可靠估计通常需要多个近似独立的数据段、适当加窗、重叠和谱平均，并报告频率分辨率与置信水平。 #### 误区：GCC-PHAT 一定优于普通互相关 PHAT 通过对幅度做归一化提高了峰值锐度，但在低 SNR 频段幅度小，除法会放大噪声相位；对稀疏谱信号或强色散信道也可能引入伪峰。因此需要正则化 $\varepsilon$、频带阈值或与普通互相关、Roth、SCOT 加权做对比。方法选择应由信号-噪声频谱结构决定，而不是"越复杂越好"。 #### 误区：$\widehat H_1$ 和 $\widehat H_2$ 在真实数据上应完全相等 只有噪声与假设完全匹配（$\widehat H_1$ 对应输入无噪，$\widehat H_2$ 对应输出无噪）时二者才一致。真实数据两端都可能有噪声，二者往往不同；差异大小可以反过来估计噪声分布方向，也是相干性显著低于 1 的显式表征。 #### 误区：匹配滤波器的最优 SNR 与波形形状有关 在 AWGN 下，匹配滤波输出 SNR 为 $2E_s/N_0$，只依赖模板总能量和噪声 PSD，与脉冲形状无关。改变形状影响的是时延分辨率、多普勒容差和旁瓣结构，而不是峰值 SNR 本身。 #### 误区：相位斜率法可以在任意频段读时延 只有在相干性足够高、$S_{xy}$ 幅度足够大的频段，互谱相位才稳定；低相干或谱功率不足处相位由噪声主导，斜率毫无意义。工程上应先按 $\gamma^2$ 门限筛频段，再解缠绕做线性回归。 --- ### 相关性与相干性的典型例子 #### 例一：固定时延的同一信号 设 $$ y(t)=a\,x(t-\tau_0), $$ 且没有噪声。若 $x$ 是宽平稳过程，则互相关是自相关的平移和缩放： $$ R_{xy}(\tau)=a^*R_{xx}(\tau-\tau_0) $$ 或按另一约定出现相反平移。互相关峰给出 $\tau_0$。 频域中： $$ Y(f)=aX(f)e^{-j2\pi f\tau_0}. $$ 理论 MSC 为 $$ \gamma_{xy}^2(f)=1 $$ （在 $S_{xx}(f)>0$ 的频率），而互谱相位是一条关于频率的直线。这个例子说明，固定时延会改变零延迟相关和相位，但不会降低理想线性相干性。

#### 例二：输出叠加不相关噪声

设

$$
y(t)=x(t)+n(t),
$$

其中 $n$ 与 $x$ 不相关。则

$$
S_{xy}=S_{xx},
\qquad
S_{yy}=S_{xx}+S_{nn}.
$$

所以

$$
\gamma_{xy}^2(f)
=\frac{S_{xx}(f)}{S_{xx}(f)+S_{nn}(f)}.
$$

在信号主导频带，相干性接近 1；在噪声主导频带，相干性接近 0。这比单个总体相关系数更清楚地显示噪声在哪些频率破坏了线性联系。

#### 例三：共同输入造成伪直接关系

设

$$
x=h_1*u+n_x,
\qquad
y=h_2*u+n_y.
$$

即 $x$、$y$ 都由隐藏源 $u$ 驱动。二者可能具有很高的互相关和相干性，但并不表示 $x$ 直接导致 $y$。若同时测得 $u$，可进一步计算偏相干性，分析控制 $u$ 后是否仍存在剩余关联。

#### 例四：非线性平方关系

设零均值对称过程 $x(t)$，并令

$$
y(t)=x^2(t).
$$

$y$ 完全由 $x$ 决定，但线性相关可能为零，同频 MSC 也可能很低，因为平方运算把能量搬移到 DC、二次谐波和频率和差项。此时应考虑高阶谱、bicoherence 或非线性依赖指标。

#### 例五：同频但相位跨段随机

假设每个短数据段中，$x$ 和 $y$ 都含频率 $f_0$ 的正弦分量，但每段的相位差随机变化。两个自谱在 $f_0$ 都有明显峰值，单段看起来也高度匹配；然而跨段平均时，互谱的复相位相互抵消，因此总体相干性可能很低。

这说明相干性要求的是跨平均样本稳定的相位和线性关系，而不仅是两个信号“都在同一频率有能量”。

#### 例六：部分频带线性相关

设输出由低通后的输入和高频独立噪声组成：

$$
y=h_{\mathrm{LP}}*x+n_{\mathrm{HF}}.
$$

则低频段 $y$ 主要由 $x$ 线性解释，相干性较高；高频段由独立噪声主导，相干性较低。一个全频带 Pearson 相关系数会把两部分压缩成单个数字，而相干性能够定位关联所在频带。

—
### 能量谱、功率谱与互谱的典型例子

#### 有限矩形脉冲：能量谱

设

$$
x(t)=A,\quad 0\le t **相关函数描述时延域中的共同变化，能量谱描述有限能量在频率上的分布，功率谱描述长期平均功率在频率上的分布，互能量谱和互功率谱则进一步保留两个信号之间的相对幅度与相位。**

二者都包含乘法和求和，但卷积的核心是系统响应与反转平移，相关的核心是匹配、内积、共轭和相对延迟。理解这一区别，就能正确解释滤波、匹配滤波、同步、谱分析和自适应算法中的大多数相关公式。

—