43、随机信号的相关性与功率谱分析

随机信号的相关性与功率谱分析

1. 相关函数概述

在研究随机信号时，相关函数是一个重要的概念。常见的相关函数类型在实验时间序列中经常出现，例如：
- 白噪声 ：其自相关函数 (K_{xx}(\Delta) = k\delta(\Delta)) ，这里的 (\delta(\Delta)) 是狄拉克δ函数，意味着只有在 (\Delta = 0) 时存在显著相关性，所以白噪声也被称为 delta 相关噪声。
- 正弦波 ：(K_{xx}(\Delta) = 0.5A^2 \cos(2\pi f_0\Delta)) ，其中 (A) 是振幅，(f_0) 是频率。
- 指数相关噪声 ：(K_{xx}(\Delta) = \exp(-k|\Delta|)) ，其相关函数随 (\Delta) 呈指数衰减。
- 低通滤波白噪声 ：(K_{xx}(\Delta) = k_B\frac{\sin(2\pi B\Delta)}{2\pi B\Delta}) ，其中 (k_B) 和 (B) 为常数。

2. 自相关函数的估计

通常，自相关函数 (C_{xx}(\Delta)) 可以使用快速傅里叶变换算法来计算。不过，从实验数据中测量 (C_{xx}(\Delta)) 时，需要注意以下几点：
- 采样间隔 ：时间序列是在固定时间间隔 (t_1 - t_2) 下采样的，所以 (\Delta = m(t_1 - t_2)) ，其中 (m) 是整数。
- 有限长度 ：收集的时间序列必然是有限长度的，因此估计的自相关函数最好称为样本自相关函数 (c_{xx}(\Delta)) 。那么，样本相关函数 (c_{xx}(\Delta)) 能否很好地估计 (C_{xx}(\Delta)) 呢？
- 端点效应 ：由于 (x(t)) 通常不是周期函数，需要考虑端点效应。可以通过零填充来克服这个问题。为了确定总滞后为 (R) 的相关性，需要在数据集末尾添加长度为 (R) 的零缓冲区。另外，快速傅里叶变换对于长度为 (2^n) （(n) 为正整数）的数据集最为高效，这意味着零填充的长度可以等于原始时间序列的长度。虽然两种方法在理论上都是正确的，但对于非常大的数据集，第一种方法通常能在最短的时间内计算出 (C_{xx}(\Delta)) 。因为需要考虑计算快速傅里叶变换所需的 CPU 时间以及访问和写入内存所需的时间，特别是当数据集的大小接近或超过可用 RAM 时。对于现代计算机，CPU 计算速度较快，而内存存储和访问速度较慢。

下面通过一个 mermaid 流程图展示自相关函数估计的流程：

graph TD
    A[获取实验时间序列数据] --> B[考虑采样间隔]
    B --> C[处理有限长度问题]
    C --> D[处理端点效应（零填充）]
    D --> E[使用快速傅里叶变换计算Cxx(Δ)]

图 1 比较了未去除均值和去除均值后的样本自相关函数。去除均值后的 (K_{xx}(\Delta)) 图呈现出白噪声的特征形状。

3. 相关时间

相关时间 (t_{corr}) 是描述噪声时间序列中测量噪声与白噪声差异程度的第二个指标，其定义为：
[t_{corr} = \frac{1}{C_{xx}(0)}\int_{0}^{\infty} C_{xx}(\Delta)d\Delta]
(C_{xx}(\Delta)) 随着 (\Delta) 趋近于无穷大时趋近于 0 的速率是随机过程记忆性的一种度量。例如，对于指数相关噪声 (C_{xx}(\Delta) = e^{-k\Delta}) ，其相关时间为：
[t_{corr} = \int_{0}^{\infty} e^{-k\Delta}d\Delta = \frac{1}{k}]
这表明相关时间是 (C_{xx}(\Delta)) 下降到其初始值的 (e^{-1}) 所需的时间。

4. 样本自相关函数的偏差

在 Python 中，使用 np.correlate() 计算的样本自相关函数 (c_{xx}(\Delta)) 是原始相关：
[c_{raw_{xx}}(\Delta) = \sum_{n = 0}^{N - |\Delta| - 1} x(n)x(n + \Delta)]
在统计学文献中，(c_{xx}(\Delta)) 通常有两种归一化方式：
- 有偏估计 ：
[c_{biased_{xx}}(\Delta) = \frac{1}{N}\sum_{n = 0}^{N - |\Delta| - 1} x(n)x(n + \Delta)]
- 无偏估计 ：
[c_{unbiased_{xx}}(\Delta) = \frac{1}{N - |\Delta|}\sum_{n = 0}^{N - |\Delta| - 1} x(n)x(n + \Delta)]
有偏估计 (c_{biased_{xx}}(\Delta)) 相对于无偏估计有优势。特别是在计算无偏估计时，相关序列端点处经常会出现随机的大波动，而使用有偏估计可以避免这种情况。当 (\Delta) 趋近于 (N) 时，((N - \Delta)^{-1}) 趋近于无穷大，计算和式时使用的非零项数量减少，因此相关区间端点处会出现随机且可能非常大的波动。而有偏估计 (c_{biased_{xx}}(\Delta)) 只有在 (\Delta = 0) 处存在偏差，所以这种估计方法更可取。

5. 自相关函数的显著性

对于白噪声信号，当时间序列足够长时，(\Delta \neq 0) 时 (C_{xx}(\Delta)) 的波动近似服从均值为零、方差为 (N^{-1}) 的正态（高斯）分布。因此，95% 的置信区间为 (\pm\frac{1.96}{\sqrt{N}}) 。如果对于给定的 (\Delta) ，(C_{xx}(\Delta)) 落在这个置信区间之外，那么原信号是白噪声的零假设将被拒绝。

另一种检验显著相关性的策略是使用“自助法”来估计置信区间。具体步骤如下：
1. 随机打乱时间序列数据（在 numpy 中可以使用 np.random.shuffle() 函数）以去除任何相关性。
2. 对同一时间序列进行多次（例如 100 次）打乱后，计算 (C_{xx}(\Delta)) 的包络，这个包络给出了“偶然”情况下得到的 (C_{xx}(\Delta)) 值。
3. 如果 (C_{xx}(\Delta)) 落在这个经验确定的置信区间之外，那么很可能存在显著的自相关，在这个例子中，置信水平约为 99%。

以下是一个总结表格：
|相关概念|定义或计算方法|特点或用途|
| ---- | ---- | ---- |
|白噪声自相关函数|(K_{xx}(\Delta) = k\delta(\Delta))|只有 (\Delta = 0) 时有显著相关性|
|正弦波自相关函数|(K_{xx}(\Delta) = 0.5A^2 \cos(2\pi f_0\Delta))|与振幅和频率有关|
|指数相关噪声自相关函数|(K_{xx}(\Delta) = \exp(-k|\Delta|))|相关函数随 (\Delta) 指数衰减|
|低通滤波白噪声自相关函数|(K_{xx}(\Delta) = k_B\frac{\sin(2\pi B\Delta)}{2\pi B\Delta})|与 (k_B) 和 (B) 有关|
|自相关函数估计|使用快速傅里叶变换算法|需考虑采样间隔、有限长度和端点效应|
|相关时间|(t_{corr} = \frac{1}{C_{xx}(0)}\int_{0}^{\infty} C_{xx}(\Delta)d\Delta)|衡量噪声与白噪声差异程度|
|样本自相关有偏估计|(c_{biased_{xx}}(\Delta) = \frac{1}{N}\sum_{n = 0}^{N - |\Delta| - 1} x(n)x(n + \Delta))|避免端点大波动|
|样本自相关无偏估计|(c_{unbiased_{xx}}(\Delta) = \frac{1}{N - |\Delta|}\sum_{n = 0}^{N - |\Delta| - 1} x(n)x(n + \Delta))|可能出现端点大波动|
|自相关显著性检验|95% 置信区间 (\pm\frac{1.96}{\sqrt{N}}) ；自助法|判断信号是否为白噪声|

6. 随机信号的功率谱

对于确定性时间序列，使用快速傅里叶变换算法可以方便地估计样本功率谱 (w(f)) 。但对于随机时间序列，这种方法并不适用。因为随机时间序列的样本功率谱 (w(f)) 不会随着样本长度的增加而在统计意义上收敛到一个极限值。这是因为傅里叶分析基于固定振幅、频率和相位的假设，而随机过程的时间序列具有随机变化的振幅、频率和相位。

不过，可以利用自相关函数 (C_{xx}(\Delta)) 和维纳 - 辛钦定理来解决这个问题。维纳 - 辛钦定理表明，样本功率谱密度函数 (w(f)) 和样本自相关函数 (C_{xx}(\Delta)) 构成一个傅里叶变换对：
[w(f) = \int_{-\infty}^{\infty} C_{xx}(\Delta)e^{-j2\pi f\Delta} d\Delta]
取逆傅里叶变换可得：
[C_{xx}(\Delta) = \int_{-\infty}^{\infty} w(f)e^{j2\pi f\Delta} df]
由于 (C_{xx}(\Delta)) 是偶函数，上述公式可以改写为更方便的形式：
[\hat{w}(f) = 4\int_{0}^{\infty} C_{xx}(\Delta)(\cos(2\pi f\Delta))d\Delta]
[C_{xx}(\Delta) = \int_{0}^{\infty} \hat{w}(f)(\cos(2\pi f\Delta))df]
其中 (\hat{w}) 是单边功率谱，即定义在非负频率上的功率谱。自谱密度函数 (S_{xx}(f)) 常用于强调功率谱是使用与定义自相关函数相同的功率度量来定义的，它对应于之前介绍的周期图。

这些傅里叶变换对对于确定性时间序列或随机过程的实现都是有效的。有限时间序列的样本自相关函数 (c_{xx}(\Delta)) 可以很好地估计无限长时间序列的 (C_{xx}(\Delta)) 。可以证明，(c_{xx}(\Delta)) 的均方误差为 (O(\frac{1}{T})) ，随着 (t) 趋近于无穷大，其分布更集中在 (C_{xx}(\Delta)) 附近，这是遍历假设的一种表述。但需要注意的是，遍历性适用于 (c_{xx}(\Delta)) ，而不适用于 (w(f)) 。

下面通过一个 mermaid 流程图展示利用维纳 - 辛钦定理计算功率谱的流程：

graph TD
    A[获取时间序列数据] --> B[计算自相关函数cxx(Δ)]
    B --> C[根据维纳 - 辛钦定理计算w(f)]
    C --> D[得到功率谱]

例如，计算正弦波过程 (x(t)=A\cos(2\pi f_0t)) 的自谱密度函数 (S_{xx}(f)) ：
[S_{xx}(f) = \frac{A^2}{2}\int_{-\infty}^{\infty} e^{-j2\pi f\Delta} \cos(2\pi f_0\Delta)d\Delta]
利用欧拉关系可得：
[S_{xx}(f) = \frac{A^2}{4}\int_{-\infty}^{\infty} (e^{-j2\pi(f + f_0)\Delta} + e^{-j2\pi(f - f_0)\Delta})d\Delta = \frac{A^2}{4} [\delta(f + f_0) + \delta(f - f_0)]]

在实际中，实验室收集的信号通常包含确定性和随机成分。对于这种情况，目前关于估计 (W(f)) 的最佳策略研究较少。如果只是为了估计信号中周期成分的频率，可以使用周期图或功率谱密度。但如果需要准确估计功率，可以先将自相关函数拟合到合适的函数形式，然后使用维纳 - 辛钦定理手动计算 (W(f)) 。

7. 噪声类型

噪声在许多领域都有研究，不同类型的噪声具有不同的特点和应用。以下是几种常见的噪声类型：
- 白噪声 ：自相关函数是 delta 函数，即 (\Delta = 0) 时为 1，(\Delta \neq 0) 时为 0。通常所说的白噪声是均匀分布的，其傅里叶变换是常数。但这种噪声在物理上是不现实的，因为其功率是无限的。不过在小频率范围内，可以近似将噪声视为白噪声。对其进行积分可得到模拟布朗运动的布朗噪声，对其进行微分可得到紫罗兰噪声。另一种常见的白噪声是高斯分布的白噪声，根据中心极限定理，大量统计独立的随机变量之和在一定条件下近似服从高斯分布，因此在系统 - 环境模型中，常将环境表示为高斯分布的白噪声，其傅里叶变换也是高斯的。
- 指数相关噪声 ：相关函数 (C_{xx}(\Delta)) 随 (\Delta) 呈指数衰减，理论上，考虑惯性效应时的布朗运动可以用指数相关噪声来描述，这种情况称为奥恩斯坦 - 乌伦贝克过程。
- 有色噪声 ：相关函数非指数衰减的噪声。特别关注的是幂律噪声，其功率谱形式为 (\frac{1}{f^{\beta}}) ，其中 (\beta) 是常数。在对数 - 对数坐标下，功率与频率的关系是线性的，斜率为 (-\beta) 。颜色术语是基于声音频率谱和光波频率谱的松散类比。例如，白噪声（(\beta = 0) ）的声波模式转换为光波时，观察者会将其视为白光。
- 粉红噪声（(\beta = 1) ） ：功率密度每十倍频程下降 10 dB，每倍频程下降 3 dB。在对数空间中是线性的，在成比例宽度的频率带中具有相等的功率。例如，在 40 - 60 Hz 和 4000 - 6000 Hz 频率带中具有相同的功率。大多数人认为粉红噪声的声音是均匀的，因为它在每个倍频程上具有相等的功率，这表明人类的听觉是在比例空间中进行的。因此，粉红噪声被音频行业广泛用作参考信号，例如在立体声系统的设计中。此外，粉红噪声在低频（DC）和高频端的功率谱中有无限多个对数频带，每个有限能量谱在两端的能量都小于粉红噪声，这是它与其他有色噪声的区别。当 (\beta > 1) 时，功率谱在高频端积分是有限的；当 (\beta < 1) 时，在低频（DC）端积分是有限的。

以下是一个总结表格：
|噪声类型|自相关函数特点|功率谱特点|应用或特点|
| ---- | ---- | ---- | ---- |
|白噪声（均匀分布）|(K_{xx}(\Delta) = k\delta(\Delta))|傅里叶变换是常数|小频率范围近似，积分得布朗噪声，微分得紫罗兰噪声|
|白噪声（高斯分布）| - |傅里叶变换是高斯的|系统 - 环境模型中表示环境|
|指数相关噪声|(C_{xx}(\Delta)=\exp(-k|\Delta|))| - |考虑惯性效应的布朗运动|
|粉红噪声|(S_{xx}(f)\propto\frac{1}{f})|功率密度每十倍频程下降 10 dB，每倍频程下降 3 dB|音频行业参考信号|

8. 互相关函数

互相关函数 (C_{xy}(\Delta)) 的定义与自相关函数 (C_{xx}(\Delta)) 类似，但考虑的是两个函数 (x(t)) 和 (y(t)) ：
[C_{xy}(\Delta) = \int_{-\infty}^{\infty} x(t)y(t - \Delta)dt = \int_{-\infty}^{\infty} x(t + \Delta)y(t)dt]
对于实过程，(C_{xy}(\Delta)) 满足 (C_{xy}(\Delta) = C_{yx}(-\Delta)) ，但一般情况下 (C_{xy}(\Delta) \neq C_{xy}(-\Delta)) 。当 (x) 和 (y) 相关时，(C_{xy}(\Delta)) 通常关于 (\Delta = 0) 不对称，这种不对称的常见原因是存在时间延迟。

互相关函数很有用，它可以利用白噪声来估计脉冲函数和时间延迟。此外，结合自相关函数，可以定义相干函数来确定生物系统的线性程度。

当需要对非零互相关 (C_{xy}(\Delta)) 进行统计检验时，需要从 (x(t)) 和 (y(t)) 中去除所有因果信息，即 (C_{xx}(\Delta)) 和 (C_{yy}(\Delta)) 都必须是 delta 相关的，这个过程称为预白化。常见的预白化策略有两种：
- 差分滤波 ：使用差分滤波器 (d(t) = x(t + 1) - x(t)) （在计算机程序中通常称为 diff() ）。可能需要对时间序列进行多次差分才能得到 delta 相关的 (C_{xx}(\Delta)) 。
- 线性自回归（AR）模型 ：首先将原始时间序列拟合到合适的线性自回归模型 (\hat{X} t = \alpha_0 + \alpha_1X {t - 1} + \alpha_2X_{t - 2} + \cdots + \alpha_nX_{t - n}) ，其中 (\alpha_i) 是常数，(n) 是 (C_{xx}(\Delta)) 中最大的显著滞后。原始数据 (X_t) 与拟合值 (\hat{X}_t) 的残差 (X_t - \hat{X}_t) 构成预白化时间序列。这种预白化形式常用于数据压缩，例如手机不直接传输数字化语音，而是先对语音信号进行白化处理，然后传输白化后的时间序列和从 AR 模型中得到的滤波器系数。当样本量较小时，需要注意这些过程不会引入偏差。

9. 脉冲响应

对于输入为白噪声的线性系统，其输出 (y(t)) 可以表示为：
[y(t) = \int_{0}^{\infty} I(s)x(t - s)ds + \xi(t) = \int_{- \infty}^{t} I(t - v)x(v)dv + \xi(t)]
其中 (x(t)) 是输入，(y(t)) 是输出，(\xi(t)) 是与 (x(t)) 不相关的噪声项。

将上式两边乘以 (x(t - u)) 并取平均可得：
[\langle x(t - \Delta)y(t) \rangle = \int_{0}^{\infty} I(s)\langle x(t - \Delta)x(t - s) \rangle ds + \langle x(t - \Delta)\xi(t) \rangle]
即
[C_{xy}(\Delta) = \int_{0}^{\infty} I(s)C_{xx}(s - \Delta)ds]
由于假设 (x(t)) 和 (\xi(t)) 不相关，最后一项为零。如果输入 (x(t)) 是方差为 (\sigma^2) 的白噪声，则 (C_{xx}(s - \Delta) = \sigma^2\delta(s - \Delta)) ，最终可得：
[I(\Delta) = \frac{C_{xy}(\Delta)}{\sigma^2}]
这意味着可以通过白噪声输入和相关输出的互协方差来获得系统的脉冲函数，这种估计脉冲函数的方法称为互相关分析。

下面通过一个 mermaid 流程图展示脉冲响应估计的流程：

graph TD
    A[输入白噪声x(t)到线性系统] --> B[得到输出y(t)]
    B --> C[计算互相关函数Cxy(Δ)]
    C --> D[根据Cxy(Δ)和σ²计算脉冲函数I(Δ)]

综上所述，随机信号的相关性和功率谱分析在许多领域都有重要应用。通过对自相关函数、互相关函数、功率谱等的研究，可以更好地理解随机信号的特性，为信号处理、系统建模等提供有力支持。同时，不同类型的噪声也具有各自独特的性质和应用场景，在实际应用中需要根据具体情况进行选择和处理。