使用 numpy.fft VS 计算信号 PSD 的前置因素。scipy.signal.welch

Question

St信号的功率谱密度u可以计算为信号的 FFT 与其u_fft复共轭 的乘积u_fft_c。在 Python 中，这可以写成：

import numpy as np

u = # Some numpy array containing signal
u_fft = np.fft.rfft(u-np.nanmean(u))

St = np.multiply(u_fft, np.conj(u_fft))

然而，Numpy 中的 FFT 定义需要将结果与因子 相乘1/N，其中N=u.size，以便在 u 与其 FFT 之间实现能量上一致的变换。这导致使用 numpy 的 fft 修正 PSD 的定义：

St = np.multiply(u_fft, np.conj(u_fft))
St = np.divide(St, u.size)

另一方面，Scipy 的函数signal.welch直接根据输入计算 PSD u：

from spicy.signal import welch

freqs_st, St_welch = welch(u-np.nanmean(u), 
     return_onesided=True, nperseg=seg_size, axis=0)

生成的 PSD是通过在大小为 的St_welch数组段中执行多次 FFT 获得的。因此，我的问题是：useg_size

是否应该St_welch乘以 因子1/seg_size才能得到能量一致的 PSD？应该乘以吗1/N？根本不应该相乘吗？

PD：通过对信号执行这两种操作进行比较并不简单，因为韦尔奇方法还引入了信号平滑并改变了频域中的显示。

有关使用时前置因素必要性的信息numpy.fft：

关于此事的期刊文章

Answer 1

参数的定义scale表明scipy.signal.welch适当的缩放由函数 执行：

\n

\n

缩放 : { \xe2\x80\x98密度\xe2\x80\x99, \xe2\x80\x98spectrum\xe2\x80\x99 }, 可选\n在计算功率谱密度 (\xe2\x80\x98密度\xe2\x80 \x99)，其中 Pxx 的单位为 V^2/Hz，并计算功率谱 (\xe2\x80\x98spectrum\xe2\x80\x99)，其中Pxx 的单位为 V^2，如果 x 以 V 为单位测量且 fs 为以赫兹为单位测量。默认为 \xe2\x80\x98密度\xe2\x80\x99

\n

\n

正确的采样频率将作为参数提供，fs以检索正确的频率和准确的功率谱密度。\n要恢复与使用 计算的功率谱类似的功率谱np.multiply(u_fft, np.conj(u_fft))，必须分别提供 fft 帧的长度和应用的窗口作为框架的长度和boxcar（相当于根本没有窗口）。scipy.signal.welch可以通过测试正弦波来检查是否应用了正确的缩放：

\n

import numpy as np\nimport scipy.signal\n\nimport matplotlib.pyplot as plt\n\n\ndef abs2(x):\n    return x.real**2 + x.imag**2\n\nif __name__ == \'__main__\':\n    framelength=1.0\n    N=1000\n    x=np.linspace(0,framelength,N,endpoint=False)\n    y=np.sin(44*2*np.pi*x)\n    #y=y-np.mean(y)\n    ffty=np.fft.fft(y)\n    #power spectrum, after real2complex transfrom (factor )\n    scale=2.0/(len(y)*len(y))\n    power=scale*abs2(ffty)\n    freq=np.fft.fftfreq(len(y) , framelength/len(y) )\n\n    # power spectrum, via scipy welch. \'boxcar\' means no window, nperseg=len(y) so that fft computed on the whole signal.\n    freq2,power2=scipy.signal.welch(y, fs=len(y)/framelength,window=\'boxcar\',nperseg=len(y),scaling=\'spectrum\', axis=-1, average=\'mean\')\n    \n    for i in range(len(freq2)):\n        print i, freq2[i], power2[i], freq[i], power[i]\n    print np.sum(power2)\n    \n    \n    plt.figure()\n    plt.plot(freq[0:len(y)/2+1],power[0:len(y)/2+1],label=\'np.fft.fft()\')\n    plt.plot(freq2,power2,label=\'scipy.signal.welch()\')\n    plt.legend()\n    plt.xlim(0,np.max(freq[0:len(y)/2+1]))\n\n\n    plt.show()\n

\n

对于实数到复数的变换，正确的缩放比例np.multiply(u_fft, np.conj(u_fft))是 2./(u.size*u.size)。事实上， 的缩放比例u_fft是 1./u.size。此外，实数到复数变换仅报告一半的频率，因为 bin Nk 的幅度将是 bin k 的幅度的复共轭。因此，该仓的能量等于仓 k 的能量，并且将其与仓 k 的能量相加。因此，因子为 2。对于幅度为 1 的测试正弦波信号，能量报告为 0.5：它确实是幅度为 1 的平方正弦波的平均值。

\n

如果帧的长度不是信号周期的倍数或者信号不是周期性的，则加窗很有用。如果信号由阻尼波组成，则使用较小的 fft 帧很有用：信号可以被视为在特征时间上具有周期性：选择小于该特征时间但大于波的周期的 fft 帧似乎是明智的。

\n