经典谱估计（功率谱密度）方法及对应python代码

1. 自相关函数估计功率谱（1958年Blackman和Tukey，BT法）

一个离散随机信号在各时间点上的值是不能先验确定的，它的每个实现（样本）往往无法如确定信号那样可以用数学表达式精确地表示，而只能用它的各种统计量来表征。其中，离散随机信号的自相关序列作为时移函数能较完整地表征它的特定统计量。由维纳-欣钦定理可知，一个随机信号的功率谱密度（函数）是自相关序列的傅氏变换。对于一个随机信号来讲，它本身的傅氏变换是不存在的，只能用功率谱密度来表征它的统计平均谱特性。因此功率谱密度是随机信号的一种重要表征形式。

2. 周期图法-直接法

利用对自相关序列的估计进行傅氏变换求得功率谱估计的方法即为经典的谱估计方法。这种方法又可分成二种。一种是间接方法，它先对自相关序列进行估计（一般都需要引入窗函数将自相关序列值加权，以减小自相关序列截段的影响），然后再对其作傅里叶变换得功率谱估计。另一种是直接法，它根据观察到的有限样本数据 ，利用离散傅里叶变换来直接进行功率谱估计而不通过对自相关序列的估计，这种方法称为周期图法。

3. 改进周期图

3.1 平均周期图法（Bartlett 法）

将周期图进行平均处理是先将原始观测数据分成若干段，再求各段周期图的平均值，将其作为功率谱的估计，这就是平均周期图法。

3.2 平均修正周期图法（Welch 法）

平均修正周期图法通过两个步骤实现对平均周期图法的改进。第一步，先将有限长的观测数据分割成重叠的分段数据。这样做的目的是增加平均周期图法中的段数，以减少功率谱估计的方差；同时，还可以减轻上一节提到的谱估计分辨率和谱估计方差对于分段数目有矛盾的依赖。

3.3 平滑周期图法（Blackman-Tukey 法）

平滑周期图法通过引入窗函数对信号的自相关序列有偏估计加权，然后对其进行傅里叶变换，得到功率谱估计，可以改善估计质量。

4. 简单比较

4.1 自相关法：

优势：

• 理论简单： 自相关法的理论基础相对简单，易于理解和实现。
• 适用于非平稳信号： 对于非平稳信号，自相关法有一定的优势，因为它考虑了信号在不同时间点的相关性。

劣势：

• 计算复杂度高： 自相关法的计算复杂度通常较高，特别是在信号较长的情况下。
• 对噪声敏感： 自相关法对噪声较为敏感，可能导致估计结果的不稳定性。

4.2 周期图法：

优势：

• 适用于周期信号： 周期图法在处理周期性信号时表现良好，能够有效地反映信号的周期性特征。
• 频谱解析较好： 在信号具有明显频谱特征时，周期图法能够提供相对准确的频谱解析。

劣势：

• 对窗函数敏感： 周期图法对窗函数的选择敏感，选择不当可能导致频谱泄漏或分辨率不足。
• 不适用于非平稳信号： 对于非平稳信号，周期图法可能无法提供准确的频谱信息。

4.3 Bartlett法：

优势：

• 减小噪声影响： Bartlett法通过对信号分段并取平均，能够有效减小噪声的影响，提高估计的稳定性。

劣势：

• 频率分辨率相对较低： 由于取平均的过程，Bartlett法的频率分辨率相对较低，可能无法捕捉信号的细节特征。

4.4 Welch法：

优势：

• 频谱平滑： 引入窗函数后，Welch法能够在一定程度上平滑功率谱，减小频谱泄漏的问题。
• 适用于非平稳信号： 相对于简单的周期图法，Welch法对非平稳信号的处理更为有效。

劣势：

• 窗口选择敏感： 对于不同信号特性，选择合适的窗口函数可能需要一些经验，选择不当会影响估计的准确性。
• 分辨率降低： 窗口的引入可能导致频率分辨率的降低。

4.5 Blackman-Tukey法：

优势：

• 适用于大样本量： Blackman-Tukey法在样本量较大时可以提供相对准确的功率谱估计。
• 对非平稳信号有一定优势： 通过自相关函数的方式，Blackman-Tukey法相对适用于非平稳信号。

劣势：

• 对信号长度敏感： 信号长度对于Blackman-Tukey法的性能有一定影响，较短的信号可能导致估计的不准确性。
• 计算复杂度较高： 对于大样本量，Blackman-Tukey法的计算复杂度可能相对较高。

5. 代码实现

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 生成示例信号
fs = 1000  # 采样频率
t = np.linspace(0, 1, fs, endpoint=False)  # 时间数组
signal1 = np.sin(2 * np.pi * 50 * t)  # 正弦波
noise = 0.5 * np.random.randn(fs)  # 噪声
signal_with_noise = signal1 + noise  # 合成信号

# 创建纵向多个子图
fig, axes = plt.subplots(nrows=5, ncols=1, figsize=(10, 15), sharex=True)

# 方法1: 周期图法
fft_result = np.fft.fft(signal_with_noise)
frequencies_fft = np.fft.fftfreq(len(fft_result), 1/fs)
psd_fft = (np.abs(fft_result)**2 / len(fft_result))**2
axes[0].semilogy(frequencies_fft[:len(frequencies_fft)//2], psd_fft[:len(psd_fft)//2], label='FFT', color='blue')
axes[0].set_title('FFT Power Spectral Density')
axes[0].set_xlim(0, 500)  # 设置 x 轴范围
axes[0].set_ylim(1e-6, 1e3)  # 设置 y 轴范围

# 方法2: 自相关函数
lags = np.arange(-len(signal_with_noise)//2, len(signal_with_noise)//2)
autocorr = np.correlate(signal_with_noise, signal_with_noise, mode='full')
autocorr /= np.max(autocorr)  # 归一化自相关函数
frequencies_acf = np.fft.fftfreq(len(autocorr), 1/fs)
psd_acf = np.abs(np.fft.fft(autocorr))
axes[1].semilogy(frequencies_acf[:len(frequencies_acf)//2], psd_acf[:len(psd_acf)//2], label='Autocorrelation', color='green')
axes[1].set_title('Autocorrelation Power Spectral Density')
axes[1].set_xlim(0, 500)  # 设置 x 轴范围
axes[1].set_ylim(1e-6, 1e3)  # 设置 y 轴范围

# 方法3: Welch方法
frequencies_welch, psd_welch = signal.welch(signal_with_noise, fs, nperseg=1024, nfft=1024)
axes[2].semilogy(frequencies_welch, psd_welch, label='Welch', color='red')
axes[2].set_title('Welch Power Spectral Density')
axes[2].set_xlim(0, 500)  # 设置 x 轴范围
axes[2].set_ylim(1e-6, 1)  # 设置 y 轴范围

# 方法4: Bartlett方法
frequencies_bartlett, psd_bartlett = signal.welch(signal_with_noise, fs, nperseg=1024, nfft=1024, window='bartlett')
axes[3].semilogy(frequencies_bartlett, psd_bartlett, label='Bartlett', color='purple')
axes[3].set_title('Bartlett Method for Power Spectral Density')
axes[3].set_xlabel('Frequency (Hz)')
axes[3].set_ylabel('Power/Frequency (dB/Hz)')
axes[3].set_xlim(0, 500)  # Set x-axis range
axes[3].set_ylim(1e-6, 1)  # Set y-axis range

# 方法5: Blackman-Tukey方法
frequencies, spectrum = signal.periodogram(signal_with_noise, fs, window='blackman', scaling='spectrum', nfft=1024)
axes[4].semilogy(frequencies, spectrum, label='Blackman-Tukey', color='orange')
axes[4].set_title('Blackman-Tukey Method for Power Spectral Density')
axes[4].set_xlabel('Frequency (Hz)')
axes[4].set_ylabel('Power/Frequency (dB/Hz)')
axes[4].set_xlim(0, 500)  # Set x-axis range
axes[4].set_ylim(1e-6, 1)  # Set y-axis range

# 调整子图布局
plt.tight_layout()
plt.show()

6. 现代谱估计方法（展望）

6.1 基于参数建模的功率谱估计

AR模型-自回归模型

MA模型-移动平均模型

ARMR模型-自回归-移动平均模型

6.2 基于非参数建模的功率谱估计

参考

https://blog.51cto.com/u_16213662/8896248

https://blog.51cto.com/u_16213691/7506006

python计算多窗谱估计工具包

https://zhuanlan.zhihu.com/p/132635425

https://blog.nowcoder.net/n/5927a14685ab4fc6bbf108775ceb3ad4?from=nowcoder_improve