实测值在对数轴上的功率谱

Question

我已经读过许多关于这个主题的讨论(lomb与fft的比较、绘制蟒蛇的功率谱、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)')

在这里，(无用的)情节产生：

​

Answer 1

下面是上面代码的一个改进版本：

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)。

我不想阻止这个话题，所以我会在这里问:我怎样才能检索到峰值的频率？如果把频率画在峰旁就很好了。