python数字信号处理之信号功率谱（PSD)计算：welch方法（分段加窗平均周期图）、Bartlett方法（周期图）的原理到代码实现

welch方法原理说明

welch方法[1]通过将数据划分为重叠的段，计算每个段的进行修改(加窗)后的周期图，然后对所有段的周期图求和进行平均，得到最终的功率谱密度。

Python和Matlab中均存在welch函数。welch函数通过配置noverlap为0，可以退化为Bartlett方法。

welch 函数参数说明：

• nperseg: 每一段的点数
• noverlap: 段与段之间重叠的点数，如果没设置就默认设置为noverlap = nperseg // 2 。如果noverlap 为0，则方法转变为Bartlett方法[2]。
• nfft: 如果需要零填充 FFT，则使用的 FFT 长度。如果为“none”，FFT 长度为 “nperseg”。默认为 “none”。
• window: 希望使用的窗。
  • 如果 “window” 是字符串或元组，则是传递给 get_window 以生成窗口值，默认情况下为 DFT 偶数。 可以从get_window函数得到窗和必需的参数列表。
  • 如果窗口是 array_like 它将被使用直接作为窗口及其长度必须为 nperseg。默认值为hanning窗。

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

matplotlib.use('TkAgg')

rng = np.random.default_rng()

#     Generate a test signal, a 2 Vrms sine wave at 1234 Hz, corrupted by
#     0.001 V**2/Hz of white noise sampled at 10 kHz.

fs = 10e3
N = 1e5
amp = 2 * np.sqrt(2)
freq = 1234.0
noise_power = 0.001 * fs / 2
time = np.arange(N) / fs
x = amp * np.sin(2 * np.pi * freq * time)
x += rng.normal(scale=np.sqrt(noise_power), size=time.shape)

# Compute and plot the power spectral density.
f, Pxx_den = signal.welch(x, fs, nperseg=1024)
plt.semilogy(f, Pxx_den)
plt.ylim([0.5e-3, 1])
plt.xlabel('frequency [Hz]')
plt.ylabel('PSD [V**2/Hz]')
plt.grid(True)
plt.show()

运行结果：

---

1. P. Welch, "The use of the fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms", IEEE Trans. Audio Electroacoust. vol. 15, pp. 70-73, 1967... ↩︎

2. M.S. Bartlett, "Periodogram Analysis and Continuous Spectra", Biometrika, vol. 37, pp. 1-16, 1950. ↩︎