使用 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 频率。

如果更改测试频率列表，您会看到在 angular ~138 频率附近有一个峰值，对应于 22.28 的频率（=angular 频率 /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.

使用 scipy.signal.lombscargle

Use of scipy.signal.lombscargle

python

signal-processing

numpy

astronomy

scipy