我已经读了很多关于这个话题的讨论(comparison between lomb-scargle and fft ,Plotting power spectrum in python,Scipy/Numpy FFT Frequency Analysis,还有很多其他的),但是仍然无法处理,所以我需要一些提示。
我有一个光子事件列表(探测与时间),数据可用here。这些列分别是time
、counts
、errors
,并在不同的能带中计数(您可以忽略它们)。我知道震源在8.9 days = 1.3*10^-6 Hz
附近有一个周期。
我想绘制功率谱密度,在这个频率上显示一个峰值(可能在log 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)')
在这里(无用的)情节产生了:
最佳答案
下面是上述代码的改进版本:
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()
在这里产生的情节(正确地):
最高峰是我试图重现的高峰。第二个峰值也在预料之中,但功率较小(事实上是这样)。
如果使用
rfft
而不是fft
(和rfftfreq
而不是fftfreq
),则复制相同的图(在这种情况下,可以使用频率值而不是模块numpy.fft.rfft)我不想屏蔽这个话题,所以我要问:我如何才能检索到峰值的频率?把频率和峰值放在一起很好。
关于python - 使用fftpack在对数轴上的实际数据功率谱,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/21542194/