我指的是以下帖子:Using scipy.signal.spectral.lombscargle for period discovery
我知道在某些情况下给出的答案是正确的。在
# 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)
然而,下面的情况似乎出了问题。在
# 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
似乎不正确。我错过了哪一步?在
您理所当然地期望峰值出现在
1 / pi
,但是您测试的最高频率是1 / 2 / pi
。。。尝试以下单个更改:现在输出是预期的:
^{pr2}$不过,请注意,如果您将
periodogram
与freqs / 2 / np.pi
对应,则图形如下所示:所以对于更复杂的信号,你不能仅仅依靠寻找周期图的
max
来找到主频,因为谐波可能会欺骗你。在相关问题 更多 >
编程相关推荐