正确使用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

但是,以下情况似乎出错了.

罪的频率(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,图形看起来是这样的:

因此,对于更复杂的信号,您不能仅仅依靠寻找max周期图来找到主导频率,因为谐波可能会欺骗您.