3

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

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)

得到的 PSDSt_welch是通过在u大小为 的阵列段中执行几个 FFT 获得的seg_size。因此，我的问题是：

应该St_welch乘以一个因子1/seg_size来给出一个能量一致的PSD？是否应该乘以1/N? 根本不应该倍增吗？

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

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

关于此事的期刊文章

4

1 回答 1

3

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

scaling : { 'density', 'spectrum' }, optional 在计算 Pxx 单位为 V^2/Hz 的功率谱密度（'密度'）和计算Pxx 单位为V^2，如果 x 以 V 为单位，fs 以 Hz 为单位。默认为“密度”</p>

将提供正确的采样频率作为参数fs以检索正确的频率和准确的功率谱密度。为了恢复类似于使用 计算的功率谱，np.multiply(u_fft, np.conj(u_fft))必须分别提供 fft 帧的长度和应用的窗口作为帧的长度和boxcar（相当于根本没有窗口）。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 Nk 的幅度将是 bin k 的幅度的复共轭。因此，该 bin 的能量等于 bin k 的能量，并将与 bin k 的能量相加。因此因子 2。对于幅度为 1 的测试正弦波信号，能量报告为 0.5：它确实是幅度为 1 的平方正弦波的平均值。

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

于 2019-09-10T19:47:59.260 回答