正确使用scipy.signal.spectral.lombscargle的方法

Question

我想说的是以下这篇帖子：使用scipy.signal.spectral.lombscargle进行周期发现

我发现这个答案在某些情况下是正确的。

sin(x)的频率是1/(2* pi)

# imports the numerical array and scientific computing packages
import numpy as np
import scipy as sp
from scipy.signal import spectral

# generates 100 evenly spaced points between 1 and 1000
time = np.linspace(1, 1000, 100)

# computes the sine value of each of those points
mags = np.sin(time)

# scales the sine values so that the mean is 0 and the variance is 1 (the documentation specifies that this must be done)
scaled_mags = (mags-mags.mean())/mags.std()

# generates 1000 frequencies between 0.01 and 1
freqs = np.linspace(0.01, 1, 1000)

# computes the Lomb Scargle Periodogram of the time and scaled magnitudes using each frequency as a guess
periodogram = spectral.lombscargle(time, scaled_mags, freqs)

# returns the inverse of the frequence (i.e. the period) of the largest periodogram value
print "1/2pi = " + str(1/(2*np.pi))
print "Frequency = " + str(freqs[np.argmax(periodogram)] / 2.0 / np.pi)

下面的结果被打印出来了。看起来没问题。我猜我们之所以要把lombscargle的结果除以2pi，是因为我们需要把弧度转换成频率。（f = 弧度 / 2pi）

1/2pi = 0.159154943092
Frequency = 0.159154943092

但是，接下来的情况似乎就不对了。

sin(2x)的频率是1/(pi)

# imports the numerical array and scientific computing packages
import numpy as np
import scipy as sp
from scipy.signal import spectral

# generates 100 evenly spaced points between 1 and 1000
time = np.linspace(1, 1000, 100)

# computes the sine value of each of those points
mags = np.sin(2 * time)

# scales the sine values so that the mean is 0 and the variance is 1 (the documentation specifies that this must be done)
scaled_mags = (mags-mags.mean())/mags.std()

# generates 1000 frequencies between 0.01 and 1
freqs = np.linspace(0.01, 1, 1000)

# computes the Lomb Scargle Periodogram of the time and scaled magnitudes using each frequency as a guess
periodogram = spectral.lombscargle(time, scaled_mags, freqs)

# returns the inverse of the frequence (i.e. the period) of the largest periodogram value
print "1/pi = " + str(1/(np.pi))
print "Frequency = " + str(freqs[np.argmax(periodogram)] / 2.0 / np.pi)

下面的结果被打印出来了。

1/pi = 0.318309886184
Frequency = 0.0780862900972

看起来不正确。我是不是漏掉了什么步骤？

Answer 1

你期待看到的峰值应该出现在 1 / pi 这个位置，但你测试的最高频率其实是 1 / 2 / pi。试试下面这个简单的修改：

freqs = linspace(0.01, 3, 3000)

现在输出结果就符合预期了：

1/pi = 0.318309886184
Frequency = 0.318311478264

不过要注意，如果你把 periodogram 和 freqs / 2 / np.pi 画成图，结果会是这样的：

所以对于更复杂的信号，你不能仅仅依靠查看 periodogram 的最大值来找到主频率，因为谐波可能会让你产生误导。