【发布时间】：2020-04-11 20:37:21
【问题描述】：
我正在努力正确归一化功率谱密度（及其倒数）。
我遇到了一个真正的问题，假设加速度计的读数以功率谱密度 (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
【问题讨论】：
标签： python fft normalization dft spectral-density