2

Python3 を使用してパワー スペクトルを計算したいと思います。このトピックに関する別のスレッドから、基本的な材料を入手しました。私はそれが次のようなものであるべきだと思います:

ps = np.abs(np.fft.fft(x))**2
timeres = t[1]-t[0]
freqs = np.fft.fftfreq(x.size, timeres)
idx = np.argsort(freqs)
plt.plot(freqs[idx], ps[idx])
plt.show()

ここtに時間とx光子数があります。私も試しました：

W = fftfreq(x.size, timeres=t[1]-t[0])
f_x = rfft(x)
plt.plot(W,f_x)
plt.show()

しかし、どちらもほとんどゼロ付近のピークを示します (ただし、同じではありません)。これからパワースペクトルを計算しようとしています：

データ

これにより、580Hz付近の信号が得られるはずです。ここで何が間違っていますか？

4

2 に答える 2

5

@kwinkunksの回答に欠けていると感じることがいくつかあります。

  1. あなたはゼロで大きなピークを見たと言いました。上記のコメントで述べたように、これは、入力信号の平均がゼロでない場合に予想されます。DC 成分を取り除きたい場合は、DFT を実行する前に、たとえば平均を減算することによって、信号のトレンドを除去する必要があります。

  2. スペクトル漏れの問題を回避するために、DFT を実行する前に、常にウィンドウ関数を信号に適用する必要があります。

  3. DFT のモジュラスの 2 乗を取ると、スペクトル密度の大まかな推定値が得られますが、これは信号のノイズに非常に敏感です。ノイズの多いデータに対するより堅牢な方法は、信号の複数の小さなセグメントのピリオドグラムを計算し、セグメント間でこれらを平均化することです。これは、ロバスト性の向上と引き換えに、周波数ドメインの分解能をいくらか犠牲にします。ウェルチの方法はこの原理を利用しています。

個人的にはscipy.signal.welch、上記のすべてのポイントに対応する を使用します。

from scipy.signal import welch

f, psd = welch(x,
               fs=1./(t[1]-t[0]),  # sample rate
               window='hanning',   # apply a Hanning window before taking the DFT
               nperseg=256,        # compute periodograms of 256-long segments of x
               detrend='constant') # detrend x by subtracting the mean
于 2015-11-28T22:07:20.113 に答える