2

我已经阅读了很多关于这个主题的讨论（lomb-scargle 和 fft 之间的比较、在 python 中绘制功率谱、Scipy/Numpy FFT 频率分析等），但仍然无法管理它，所以我需要一些提示。我有一个光子事件列表（检测与时间），数据可在此处获得。列是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)')

这里产生了（无用的）情节： 密度泛函

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 最高峰是我试图重现的高峰。预计还会出现第二个峰值，但功率较小（确实如此）。如果rfft使用代替fft（和rfftfreq代替fftfreq）相同的图被再现（在这种情况下，频率值，而不是模块，可以使用numpy.fft.rfft）

我不想堵住话题，所以在这里问一下：如何检索峰值的频率？在峰值旁边绘制频率会很棒。

于 2014-02-04T08:04:57.920 回答