Спектр мощности реальных данных с помощью fftpack на оси журнала

Question

Я уже много читал об этой теме (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. Я хотел бы построить плотность спектра мощности, показывающую пик на этой частоте (возможно, на логарифмической оси х). Было бы неплохо, если бы я мог избежать половины части сюжета (симметричный). Это мой код до сих пор, не так далеко, но все-таки что-то:

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)

Я не хочу, чтобы блокировать тему, поэтому я спрошу здесь: И как я могу получить частоты пиков? Было бы здорово построить частотные частоты на вершинах.