мои спектры мощности правдоподобны? сравнение между lomb-scargle и fft (scipy.signal и numpy.fft)

Может ли кто-нибудь любезно указать, почему я получаю совсем другие результаты?
Есть много пиков, которые не должны появляться. На самом деле, должен быть только один пик.
Я новичок в Python, и все комментарии о моем коде ниже приветствуются.

Тестовые данные здесь. введите описание ссылки здесь
Вы можете напрямую wget https://clbin.com/YJkwr.
Его первый столбец - время прибытия серии фотонов. Источник света испускает фотоны случайным образом. Общее время 55736 секунд и 67398 фотонов. Я пытаюсь обнаружить некоторую периодичность интенсивности света.
Мы можем сопоставить время и интенсивность света пропорционально количеству фотонов в каждом интервале времени.

Я попробовал numpy.fft и lomb-scargle из scipy.signal, чтобы получить спектр мощности, но получил очень разные результаты.

FFT

 import pylab as pl
 import numpy as np

 timepoints=np.loadtxt('timesequence',usecols=(0,),unpack=True,delimiter=",")
 binshu=50000
 interd=54425./binshu  
 t=np.histogram(timepoints,bins=binshu)[0]
 sp = np.fft.fft(t)
 freq = np.fft.fftfreq(len(t),d=interd)  
 freqnum = np.fft.fftfreq(len(t),d=interd).argsort()  
 pl.xlabel("frequency(Hz)")
 pl.plot(freq[freqnum],np.abs(sp)[freqnum])

введите описание изображения здесь

Lomb-scargle

 timepoints=np.loadtxt('timesequence',usecols=(0,),unpack=True,delimiter=",")
 binshu=50000
 intensity=np.histogram(timepoints,bins=binshu)[0].astype('float64') 
 middletime=np.histogram(timepoints,bins=binshu)[1][1:binshu+1]-np.histogram(timepoints,bins=binshu)[1][3]*0.5
 freq1=1./(timepoints.max()-timepoints.min())
 freq2=freq1*len(timepoints)
 freqs=np.linspace(freq1,freq2,1000)
 pm1=spectral.lombscargle(middletime,intensity,freqs)

 pl.xlabel("frequency(Hz)")
 pl.plot(freqs,pm1)
 pl.xlim(0.5*freq1,2*freq2)
 pl.ylim(0,250)
 pl.show()

введите описание изображения здесь

*********************************
Спасибо, Уоррен, я ценю это. Спасибо за подробный ответ.

Вы правы, есть время интеграции, здесь оно составляет около 1,7 с.
Есть много разовых воздействий. Каждая экспозиция стоит 1.7s
В одной экспозиции мы не можем точно сказать время ее прибытия.

Если временные ряды похожи на:
0 1,7 3,4 8,5 8,5

Время интеграции двух последних фотонов: 1.7s, а не (8.5-3.4)s. Так что я пересмотрю часть вашего кода.

ОДНАКО мой вопрос все еще остается. Вы настраиваете несколько параметров, чтобы получить в некоторой степени пик 0.024Hz в пике ломаного ниспадающего пятна. И вы используете это, чтобы вести ваши параметры в FFT.

Если вы не знаете число 0.024, возможно, вы можете использовать разные параметры для получения другого максимального пика?

Как гарантировать каждый раз, когда мы можем получить право num_ls_freqs? Вы можете видеть, выбираем ли мы num_ls_freqs, самые высокие пиковые сдвиги.

Если у меня много временных рядов, каждый раз мне приходится указывать разные параметры? И как их получить?

4 голоса | спросил questionhang 18 WedEurope/Moscow2013-12-18T16:32:55+04:00Europe/Moscow12bEurope/MoscowWed, 18 Dec 2013 16:32:55 +0400 2013, 16:32:55

1 ответ


0

Спектр мощности - это абсолютное значение преобразования Фурье в квадрате . Если вдобавок к этому вы избавитесь от значения DC в вашем БПФ, нанесите на график только положительную половину ваших результатов и вычислите ширину ячейки по выходным данным np.histogram вместо какого-то закодированного магического числа, вот что вы получите:

timepoints=np.loadtxt('timesequence',usecols=(0,),unpack=True,delimiter=",")
binshu=50000
t, _=np.histogram(timepoints,bins=binshu)
interd = _[1] - _[0]
sp = np.fft.fft(t)
freq = np.fft.fftfreq(len(t),d=interd)
n = len(freq)
pl.xlabel("frequency(Hz)")
pl.plot(freq[1:n//2],(sp*sp.conj())[1:n//2])

введите описание изображения здесь

Это выглядит очень похоже на вашу периодограмму. Существует несоответствие в положении большого пика, который составляет около 0,32 в спектре, основанном на БПФ, и больше, чем 0,4x в выходном сигнале ломбарда. Не уверен, что там происходит, но я не совсем понимаю ваши расчеты частоты, особенно в случае с ломбардами ...

ответил Jaime 18 WedEurope/Moscow2013-12-18T18:33:21+04:00Europe/Moscow12bEurope/MoscowWed, 18 Dec 2013 18:33:21 +0400 2013, 18:33:21

Похожие вопросы

Популярные теги

security × 330linux × 316macos × 2827 × 268performance × 244command-line × 241sql-server × 235joomla-3.x × 222java × 189c++ × 186windows × 180cisco × 168bash × 158c# × 142gmail × 139arduino-uno × 139javascript × 134ssh × 133seo × 132mysql × 132