0

この優れた回答を読んだ後、実数値データと scipy.signal.welch を介して計算されたデータの FFT ベースのパワー スペクトル密度 (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