时间戳

首页 > 解决方案 > python中离散功率谱密度的正确归一化以解决实际问题

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

问题描述

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

我遇到了一个真正的问题，假设加速度计的读数以功率谱密度 (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

标签： pythonfftnormalizationdftspectral-density

解决方案

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(um )，因此 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)。

如需更详细的说明，请阅读Welch 的原始论文。另请参阅维基百科的页面和C 中的数字食谱的第 13.4 章

推荐阅读