2

このトピックに関する多くの議論（lomb-scargle と fft の比較、 Python でのパワー スペクトルのプロット、Scipy /Numpy FFT Frequency Analysis、および他の多く）をすでに読みましたが、まだそれを管理できないため、いくつかのヒントが必要です。光子イベント (検出と時間) のリストがあります。データはこちらから入手できます。列はtime、counts、errors、および異なるエネルギー バンドのカウントです (これらは無視できます)。ソースには の周りに周期性があることがわかってい8.9 days = 1.3*10^-6 Hzます。この周波数でピークを示すパワースペクトル密度をプロットしたいと思います（おそらく対数x軸上）。プロットの半分（対称）を避けることができればそれもいいでしょう。これは今までの私のコードですが、まだ何かがあります:

import numpy as np
from scipy.fftpack import fft, rfft, fftfreq
import pylab as plt

x,y = np.loadtxt('datafile.txt', usecols = (0,1), unpack=True)
y = y - y.mean() # Removes the large value at the 0 frequency that we don't care about

f_range = np.linspace(10**(-7), 10**(-5), 1000)
W = fftfreq(y.size, d=x[1]-x[0])

plt.subplot(2,1,1)
plt.plot(x,y)
plt.xlabel('Time (days)')

f_signal = fft(y)
plt.subplot(2,1,2)
plt.plot(W, abs(f_signal))
plt.xlabel('Frequency (Hz)')

ここで（役に立たない）プロットが作成されました： DFT

4

1 に答える 1

4

上記のコードの改良版は次のとおりです。

import pyfits
import numpy as np
from scipy.fftpack import fft, rfft, fftfreq
import pylab as plt

x,y = np.loadtxt('data.txt', usecols = (0,1), unpack=True)
y = y - y.mean()

W = fftfreq(y.size, d=(x[1]-x[0])*86400)
plt.subplot(2,1,1)
plt.plot(x,y)
plt.xlabel('Time (days)')

f_signal = fft(y)
plt.subplot(2,1,2)
plt.plot(W, abs(f_signal)**2)
plt.xlabel('Frequency (Hz)')

plt.xscale('log')
plt.xlim(10**(-6), 10**(-5))
plt.show()

そして、ここでプロットが生成されました (正しく): fft 最高のピークは、私が再現しようとしていたピークです。2 番目のピークも予想されますが、電力は少なくなります (実際にはそうです)。rfftの代わりにfft(およびrfftfreqの代わりに) を使用すると、同じプロットが再現されます (その場合fftfreq、モジュールの代わりに周波数値を使用できますnumpy.fft.rfft )

トピックをブロックしたくないので、ここで質問します: ピークの周波数を取得するにはどうすればよいですか? ピークの横に周波数をプロットすると素晴らしいでしょう。

于 2014-02-04T08:04:57.920 に答える