对于实际问题，Python 中离散功率谱密度的正确归一化

Question

我正在努力解决功率谱密度（及其倒数）的正确归一化问题。

我遇到了一个实际问题，假设加速度计的读数以功率谱密度（psd）的形式（以幅度^2/Hz为单位）。我想将其转换回随机时间序列。然而，首先我想了解 PSD 的“前进”方向，即时间序列。

根据[1]，时间序列x(t)的PSD可以通过以下公式计算：

PSD(w) = 1/T * abs(F(w))^2 = df * abs(F(w))^2

其中T是x(t)的采样时间，F(w)是x(t)的傅里叶变换，df=1/T是傅里叶空间中的频率分辨率。然而，我得到的结果并不等于我使用 scipy Welch 方法得到的结果，请参见下面的代码。

第一个代码块取自 scipy.welch 纪录片：

from scipy import signal
import matplotlib.pyplot as plt

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 += np.random.normal(scale=np.sqrt(noise_power), size=time.shape)

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.show()

我注意到的第一件事是绘制的 psd 随变量 fs 变化，这对我来说似乎很奇怪。（也许我需要相应地调整 nperseg 参数？为什么 nperseg 没有自动设置为 fs ？）

我的代码如下：（请注意，我定义了自己的 fft_full 函数，该函数已经处理了正确的傅立叶变换归一化，我通过检查 Parsevals 定理对其进行了验证）。

import scipy.fftpack as fftpack

def fft_full(xt,yt):
    dt = xt[1] - xt[0]
    x_fft=fftpack.fftfreq(xt.size,dt)
    y_fft=fftpack.fft(yt)*dt
    return (x_fft,y_fft)

xf,yf=fft_full(time,x)
df=xf[1] - xf[0]
psd=np.abs(yf)**2 *df
plt.figure()
plt.semilogy(xf, psd)
#plt.ylim([0.5e-3, 1])
plt.xlim(0,)
plt.xlabel('frequency [Hz]')
plt.ylabel('PSD [V**2/Hz]')
plt.show()

不幸的是，我还没有被允许发布图像，但是这两个图看起来不一样！

如果有人能向我解释我错在哪里并一劳永逸地解决这个问题，我将不胜感激:)

[1]：等式。2.82。航天器结构设计理论与应用中的随机振动，作者：Wijker, J. Jaap，2009

Answer 1

scipy 库使用 Welch 方法来估计 PSD。该方法比仅采用离散傅里叶变换的平方模更复杂。简而言之，它的进展如下：

1. 令 x 为包含 N 个样本的输入离散信号。

2. 将 x 分割为 M 个重叠片段，使得每个片段 s _m包含 nperseg 样本，并且每两个连续片段在 noverlap 样本中重叠，使得 nperseg = K * (nperseg - noverlap)，其中 K 是整数（通常 K = 2） 。另请注意：N = nperseg + (M - 1) * (nperseg - noverlap) = (M + K - 1) * nperseg / K

3. 从每个段 s _m中减去其平均值（这会消除 DC 分量）： t _m = s _m - sum(s _m ) / nperseg

4. _{将获得的零均值段 t m}的元素乘以合适的（非对称）窗函数 h（例如 Hann 窗）的元素： u _m = t _m * h

5. _{计算所有向量 u m}的快速傅里叶变换。_{在执行这些变换之前，我们通常首先向每个向量 u m}附加许多零，使其新维度成为 2 的幂（函数 welch 的 nfft 参数用于此目的）。让我们假设 len(u _m ) = 2 ^p。在大多数情况下，我们的输入向量是实值，因此最好对真实数据应用 FFT。其结果是复值向量 v _m = rfft(u _m )，使得 len(v _m ) = 2 ^{p - 1} + 1。

6. 计算所有变换向量的平方模数：a _m = abs(v _m ) ** 2，或更有效：a _m = v _m .real ** 2 + v _m .imag ** 2

7. 按如下方式对向量 a _{m进行归一化： b m} = a _m / sum(h * h) b _m [1:-1] *= 2 （这考虑了负频率），其中 h 是维度的实数向量nperseg 包含窗口系数。对于汉恩窗，我们可以证明 sum(h * h) = 3 / 8 * len(h) = 3 / 8 * nperseg

8. _{将 PSD 估计为所有向量 b m}的平均值： psd = sum(b _m ) / M 结果是维度为 len(psd) = 2 ^{p - 1} + 1 的向量。如果我们希望所有 psd 的总和系数与加窗输入数据的均方振幅（而不是振幅平方和）匹配，则向量 psd 也必须除以 nperseg。然而，scipy 例程省略了这一步。无论如何，我们通常将 psd 表示为分贝刻度，因此最终结果为：psd_dB = 10 * log10(psd)。

更详细的描述，请阅读韦尔奇的论文原文。另请参阅维基百科页面和C 语言数值食谱第 13.4 章