0

在阅读了这个出色的答案后，它演示了实值数据的基于 FFT 的功率谱密度 (PSD) 与通过 scipy.signal.welch 计算的功率谱密度 (PSD) 的比例关系，我想知道如果输入时域数据是复杂的。

我直接复制了上面答案中的代码，包括缩放因子（“比例”），并简单地在输入正弦曲线中添加了一个随机虚部：

import matplotlib.pyplot as plt
from scipy import signal
import random

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

framelength=1.0
N=1000
x=np.linspace(0,framelength,N,endpoint=False)

y=np.sin(44*2*np.pi*x)
y = np.asarray([x + 1.j*random.random() for x in y]) # add the imaginary component
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=signal.welch(y, fs=len(y)/framelength,window='boxcar',nperseg=len(y),scaling='spectrum', axis=-1, average='mean')

print('Total power in welch psd: ', np.sum(power2), ' total power in fft psd: ', np.sum(power))

plt.figure()
plt.loglog(freq[0:len(y)//2],power[0:len(y)//2],label='(2/N^2) |np.fft.fft()|^2', alpha=0.6)
plt.loglog(freq2[0:len(y)//2],power2[0:len(y)//2],label='scipy.signal.welch()', alpha=0.6)
plt.ylim(1e-11, 10)
plt.legend()

这输出：

Total power in welch psd:  0.5847839312439229 , total power in fft psd:  1.6502770748462239

和示例数据图：

由 welch 和 fft 计算的叠加 PSD。

Welch PSD 与 FFT PSD 中的样本之间的比率为 0.5。

我认为问题在于，一旦离散傅里叶变换被赋予复数值而不是纯实数值，缩放就不再相同。但是，我不确定如何计算比例因子应该是多少。我注意到，如果我将比例因子设置为 1/(len(y)^2)，则 PSD 中的样本之间的比率为 1，但光谱中的总功率不相等。

有人可以帮助我更好地理解这一点吗？

4

0 回答 0