python - scipy.signal.lombscargle 的使用

Question

对于类(class)，我们试图使用 scipy 中的嵌入式包来证明 Lomb-Scargle 周期图的简单示例。关于如何使用此功能的文档很少，我在网上也找不到任何帮助。当我运行代码时，我得到的周期图主峰值为 ~6.3，而不是预期的 ~23.3。我们从中提取的数据是一个简单的 .dat 文件，其中包含数字列表。这是代码，对发生的事情有什么想法吗？

import scipy as sp
import math as m
import numpy as np
from scipy.signal import lombscargle
import pylab as plt
from numpy import shape


x=[]
y=[]
nout = 10000
file=open("hv878_1945.dat", 'r')

for pair in file:
    xandy=pair.split()
    x.append(float(xandy[0]))
    y.append(float(xandy[2]))


x=np.asarray(x)
y=np.asarray(y)
f = np.linspace(0.1, 50, nout)

periodogram=sp.signal.spectral.lombscargle(x,y,f)
normval = x.shape[0]
plt.plot(f, np.sqrt(4*(periodogram/normval)))
plt.show()

这是文件，如果有人想运行它:http://www.mediafire.com/download/9a0aw9nc40r4c73/hv878_1945.dat

如有任何帮助，我们将不胜感激!

Answer 1

您似乎遇到了描述的相同问题 here 。基本上，您没有测试感兴趣的频率，也没有考虑到传递给 lombscargle 的频率是角频率。

如果更改测试频率列表，您会看到角频率 ~138 附近有一个峰值，对应的频率为 22.28(=角频率/2/pi):

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import lombscargle

plt.ion()
data = np.loadtxt('hv878_1945.dat')

ang_freq = np.linspace(.1, 25*6, 10000)
periodogram = lombscargle(data[:,0], data[:,2], ang_freq)
print(ang_freq[np.argmax(periodogram)]/2/np.pi)
plt.plot(ang_freq, periodogram)

我将做出与 Jaime 在其他问题中相同的评论:

you cannot rely on just looking for the max of periodogram to find the dominant frequency, because the harmonics may fool you.