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 的频率是 angular 频率。

如果您更改测试频率列表，您会看到在 ~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 在另一个 SO 问题中所做的相同的评论:

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