当前位置：首页 > article >正文

谱分析方法

article 2025/3/19 11:10:21

前言

本文隶属于专栏《机器学习数学通关指南》，该专栏为笔者原创，引用请注明来源，不足和错误之处请在评论区帮忙指出，谢谢！

本专栏目录结构和参考文献请见《机器学习数学通关指南》

ima 知识库

知识库广场搜索：

知识库	创建人
机器学习	@Shockang
机器学习数学基础	@Shockang
深度学习	@Shockang

正文

在这里插入图片描述

自相关函数与功率谱密度：时域与频域的桥梁

在这里插入图片描述

📚 谱分析方法是研究信号频率成分的重要手段，是连接时域特性和频域分析的桥梁，在机器学习和数据科学中有着广泛应用。

🔄 一、自相关函数（Autocorrelation Function, ACF）

1.1 数学定义 📐

对于平稳随机过程 $X (t)$ ，自相关函数定义为：

$R_X(\tau) = E[X(t) X(t+\tau)]$

物理意义：描述信号在时间间隔 $\tau$ 上的自相似性。当 $\tau=0$ 时， $R_X(0)$ 表示信号的均方值（信号总功率）。

1.2 关键性质 🧩

性质	数学表达式	含义
对称性	$R_X(-\tau) = R_X(\tau)$	自相关函数是偶函数
非负定性	相关矩阵为非负定矩阵	保证功率谱密度非负
最大值特性	$\|R_X(\tau)\| \leq R_X(0)$	相关性随时间差增大而衰减

1.3 直观理解 🧠

自相关函数度量了信号与其时移版本的相似度：

高自相关：数据点间存在强依赖关系
低自相关：数据点近似独立
周期性自相关：信号中存在周期性成分

1.4 示例 🌊

对随机相位正弦波 $X(t)=a\cos(\omega t + \Theta)$ ，其中 $\Theta \sim U(0,2\pi)$ ：

$R_X(\tau) = \frac{a^2}{2} \cos(\omega \tau)$

其周期性反映信号具有单一频率分量 $\omega$ 。

📊 二、功率谱密度（Power Spectral Density, PSD）

2.1 数学定义 📐

对平稳过程，功率谱密度是自相关函数的傅里叶变换：

$S_X(\omega) = \int_{-\infty}^{\infty} R_X(\tau) e^{-j\omega \tau} d\tau$

物理意义：表征信号功率在各频率分量上的分布密度。积分全频带得总功率：

$\int_{-\infty}^{\infty} S_X(\omega) d\omega = R_X(0)$

2.2 重要性质 📋

非负性： $S_X(\omega) \geq 0$
偶对称性： $S_X(-\omega) = S_X(\omega)$
维纳-辛钦定理：自相关函数与功率谱密度构成傅里叶变换对

2.3 常见信号的谱密度 📈

信号类型	自相关函数	功率谱密度	特点
白噪声	$R_X(\tau)=\sigma^2 \delta(\tau)$	$S_X(\omega)=\sigma^2$	功率均匀分布在所有频率上
正弦波	$R_X(\tau)=\frac{a^2}{2} \cos(\omega_0 \tau)$	$S_X(\omega)=\pi a^2 [\delta(\omega-\omega_0)+\delta(\omega+\omega_0)]$	在特定频率处有脉冲

🔄 三、自相关函数与功率谱密度的关系

3.1 维纳-辛钦定理（Wiener-Khinchin Theorem）🔗

维纳-辛钦定理建立了自相关函数与功率谱密度的互换关系：

$S_X(\omega) = \mathcal{F}\{R_X(\tau)\} = \int_{-\infty}^{\infty} R_X(\tau) e^{-j\omega \tau} d\tau$

$R_X(\tau) = \mathcal{F}^{-1}\{S_X(\omega)\} = \frac{1}{2\pi}\int_{-\infty}^{\infty} S_X(\omega) e^{j\omega \tau} d\omega$

3.2 分析意义 💡

时域难以观察的周期性和噪声特性，在频域可通过谱峰/平坦度直观识别
功率谱密度将复杂信号分解成不同频率成分，便于特征提取
在机器学习中，这种转换可以揭示数据中的隐藏模式和周期性

🔄 四、各态历经性（Ergodicity）

4.1 定义与意义 📝

各态历经性：若平稳过程满足各态历经性，则可用单一样本的时间平均替代统计平均。

$\lim_{T\to\infty} \frac{1}{2T} \int_{-T}^{T} x(t)x(t+\tau) dt = R_X(\tau)$

4.2 实际应用价值 🛠️

允许从有限长度的单条时间序列估计自相关函数
为现实数据分析提供理论支撑
机器学习中，意味着可以从单一数据序列中推断整个随机过程的统计特性

💻 五、Python实现与应用

5.1 使用Python计算自相关函数 🐍

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import acf

# 生成示例数据：AR(1)过程
np.random.seed(42)
n = 1000
phi = 0.8
ar1 = np.zeros(n)
for i in range(1, n):
    ar1[i] = phi * ar1[i-1] + np.random.normal(0, 1)

# 计算自相关函数
lags = 40
acf_values = acf(ar1, nlags=lags)

# 可视化
plt.figure(figsize=(10, 6))
plt.stem(range(lags + 1), acf_values)
plt.axhline(y=0, linestyle='--', color='gray')
plt.axhline(y=1.96/np.sqrt(n), linestyle='--', color='red')
plt.axhline(y=-1.96/np.sqrt(n), linestyle='--', color='red')
plt.title('AR(1)过程的自相关函数', fontsize=14)
plt.xlabel('滞后(Lag)', fontsize=12)
plt.ylabel('自相关系数', fontsize=12)
plt.tight_layout()
plt.show()

5.2 功率谱密度计算 📊

from scipy import signal

# 生成带噪声的正弦信号
fs = 100  # 采样频率
t = np.arange(0, 10, 1/fs)  # 10秒，采样率100Hz
f1, f2 = 5, 15  # 信号频率
x = np.sin(2*np.pi*f1*t) + 0.5*np.sin(2*np.pi*f2*t) + 0.5*np.random.randn(len(t))

# 计算PSD
f, Pxx = signal.welch(x, fs, nperseg=256)

# 可视化
plt.figure(figsize=(10, 6))
plt.semilogy(f, Pxx)
plt.title('功率谱密度', fontsize=14)
plt.xlabel('频率 [Hz]', fontsize=12)
plt.ylabel('PSD [V^2/Hz]', fontsize=12)
plt.grid(True)
plt.tight_layout()
plt.show()

🚀 六、机器学习中的应用场景

6.1 信号检测与特征提取 🔍

在接收信号 $Y(t)=aS(t-\tau)+N(t)$ 中，若有用信号 $S (t)$ 与噪声 $N (t)$ 独立，其互相关函数:

$R_{SY}(\tau) = a R_S(\tau)$

应用：

雷达和声纳系统中的目标检测
生物医学信号处理中的异常检测
IoT传感器数据中的模式识别

6.2 时间序列分析与预测 📈

ARMA模型：自回归滑动平均模型的谱密度具有有理函数形式，便于参数化建模。

例如AR(2)模型： $X_t = 0.4X_{t-1} + 0.4X_{t-2} + \epsilon_t$

# AR(2)模型参数估计与频谱分析
from statsmodels.tsa.arima.model import ARIMA
import pandas as pd

# 生成AR(2)过程
n = 500
ar2 = np.zeros(n)
ar2[0:2] = np.random.randn(2)
for i in range(2, n):
    ar2[i] = 0.4*ar2[i-1] + 0.4*ar2[i-2] + np.random.normal(0, 1)

# 拟合AR模型
model = ARIMA(ar2, order=(2,0,0))
results = model.fit()
print("AR(2)模型参数估计结果:")
print(results.summary())

# 计算理论谱
ar_params = np.r_[1, -results.params[1:]]
w, h = signal.freqz(1, ar_params, worN=1000)
psd = np.abs(h)**2

# 可视化理论谱
plt.figure(figsize=(10, 6))
plt.plot(w/np.pi, psd)
plt.title('AR(2)模型的理论功率谱', fontsize=14)
plt.xlabel('归一化频率 (×π rad/sample)', fontsize=12)
plt.ylabel('功率谱密度', fontsize=12)
plt.grid(True)
plt.tight_layout()
plt.show()

6.3 深度学习中的应用 🧠

特征工程：将时域数据变换到频域作为神经网络输入特征
时频分析：通过短时傅里叶变换(STFT)获取时变频谱图作为2D卷积网络输入
异常检测：监测功率谱密度变化识别异常行为

🔬 七、高级谱分析技术

7.1 多元谱分析 🔄

交叉功率谱：研究两个信号之间的频率关系
$S_{XY}(\omega) = \mathcal{F}\{R_{XY}(\tau)\}$

相干函数：衡量两信号在特定频率上的线性相关性
$\gamma^2_{XY}(\omega) = \frac{|S_{XY}(\omega)|^2}{S_X(\omega)S_Y(\omega)}$

7.2 现代谱估计方法 📊

方法	优点	缺点	适用场景
周期图法	计算简单	方差大	长数据序列
Welch方法	降低方差	频率分辨率降低	一般应用
AR模型法	频率分辨率高	需模型阶数选择	短数据序列
多锥方法	方差小、分辨率高	计算复杂	高精度要求

# 比较不同谱估计方法
methods = ['periodogram', 'welch', 'lombscargle']
plt.figure(figsize=(12, 8))

# 生成含噪声的多频率信号
np.random.seed(0)
t = np.linspace(0, 10, 1000)
x = np.sin(2*np.pi*5*t) + 0.5*np.sin(2*np.pi*10*t) + 0.3*np.random.randn(len(t))

# 周期图法
f1, pxx1 = signal.periodogram(x, fs=100)
plt.subplot(3, 1, 1)
plt.semilogy(f1, pxx1)
plt.title('周期图法', fontsize=12)
plt.ylabel('PSD')

# Welch方法
f2, pxx2 = signal.welch(x, fs=100, nperseg=256)
plt.subplot(3, 1, 2)
plt.semilogy(f2, pxx2)
plt.title('Welch方法', fontsize=12)
plt.ylabel('PSD')

# Lomb-Scargle方法 (适用于非均匀采样)
f3 = np.linspace(0.1, 50, 500)
pxx3 = signal.lombscargle(t, x, f3)
plt.subplot(3, 1, 3)
plt.semilogy(f3, pxx3)
plt.title('Lomb-Scargle方法', fontsize=12)
plt.xlabel('频率 [Hz]')
plt.ylabel('功率')

plt.tight_layout()
plt.show()