python - python中离散功率谱密度的正确归一化实际问题

Question

我正在努力对功率谱密度(及其倒数)进行正确的归一化。

我遇到了一个真正的问题，假设加速度计的读数以功率谱密度 (psd) 的形式表示，单位为 Amplitude^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]: Eq. 2.82. Random Vibrations in Spacecraft Structures Design Theory and Applications, Authors: Wijker, J. Jaap, 2009

Answer 1

scipy 库使用 Welch 的方法来估计 PSD。这种方法比仅仅取离散傅立叶变换的平方模更复杂。简而言之，它的过程如下:

设 x 是包含 N 个样本的输入离散信号。

将 x 拆分为 M 个重叠段，使得每个段 sm 包含 nperseg 样本，并且每两个连续段在 noverlap 样本中重叠，因此 nperseg = K * (nperseg - noverlap)，其中 K 是一个整数(通常 K = 2)。另请注意:
N = nperseg + (M - 1) * (nperseg - noverlap) = (M + K - 1) * nperseg/K

从每个段 sm 中减去它的平均值(这会去除 DC 分量):
tm = sm - 总和(sm)/nperseg

将获得的零均值段 tm 的元素乘以合适的(非对称)窗函数 h(例如 Hann 窗)的元素:
嗯 = tm * h

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

计算所有变换向量的平方模:
am = abs(vm) ** 2,
或更有效地:
am = vm.real ** 2 + vm.imag ** 2

将向量归一化如下:
bm = am/sum(h * h)
bm[1:-1] *= 2(这考虑了负频率)，
其中 h 是包含窗口系数的维度 nperseg 的实向量。在 Hann 窗的情况下，我们可以证明
sum(h * h) = 3/8 * len(h) = 3/8 * nperseg

将 PSD 估计为所有向量 bm 的平均值:
psd = sum(bm)/M
结果是一个维度为 len(psd) = 2p - 1 + 1 的向量。如果我们希望所有 psd 系数的总和与加窗输入数据的均方幅度相匹配(而不是平方幅度之和)，那么向量 psd 也必须除以 nperseg。但是，scipy 例程省略了这一步。无论如何，我们通常在分贝刻度上呈现 psd，这样最终的结果是:
psd_dB = 10 * log10(psd)。

更详细的描述请阅读 original Welch's paper .另见 Wikipedia's page和 Numerical Recipes in C 的第 13.4 章