预因子计算具有 numpy.fft VS 的信号的 PSD。 scipy.signal.welch

Question

信号 u 的功率谱密度 St 可以计算为信号的 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 直接从输入 u:

计算 PSD

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 是通过在数组 u 的大小为 seg_size 的段中执行多个 FFT 获得的。因此，我的问题是：

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

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

关于使用前因子必要性的信息 numpy.fft :

Journal article on the matter

Answer 1

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

执行

scaling : { ‘density’, ‘spectrum’ }, optional Selects between computing the power spectral density (‘density’) where Pxx has units of V^2/Hz and computing the power spectrum (‘spectrum’) where Pxx has units of V^2, if x is measured in V and fs is measured in Hz. Defaults to ‘density’

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

import numpy as np
import scipy.signal

import matplotlib.pyplot as plt


def abs2(x):
    return x.real**2 + x.imag**2

if __name__ == '__main__':
    framelength=1.0
    N=1000
    x=np.linspace(0,framelength,N,endpoint=False)
    y=np.sin(44*2*np.pi*x)
    #y=y-np.mean(y)
    ffty=np.fft.fft(y)
    #power spectrum, after real2complex transfrom (factor )
    scale=2.0/(len(y)*len(y))
    power=scale*abs2(ffty)
    freq=np.fft.fftfreq(len(y) , framelength/len(y) )

    # power spectrum, via scipy welch. 'boxcar' means no window, nperseg=len(y) so that fft computed on the whole signal.
    freq2,power2=scipy.signal.welch(y, fs=len(y)/framelength,window='boxcar',nperseg=len(y),scaling='spectrum', axis=-1, average='mean')

    for i in range(len(freq2)):
        print i, freq2[i], power2[i], freq[i], power[i]
    print np.sum(power2)


    plt.figure()
    plt.plot(freq[0:len(y)/2+1],power[0:len(y)/2+1],label='np.fft.fft()')
    plt.plot(freq2,power2,label='scipy.signal.welch()')
    plt.legend()
    plt.xlim(0,np.max(freq[0:len(y)/2+1]))


    plt.show()

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

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

预因子计算具有 numpy.fft VS 的信号的 PSD。 scipy.signal.welch

Prefactors computing PSD of a signal with numpy.fft VS. scipy.signal.welch

python

numpy

scipy

spectral-density