7

私は次の投稿を参照しています:期間の発見に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

間違っているようです。私が逃したステップはありますか?

4

1 に答える 1

10

当然のことながら、ピークがに現れることを期待してい1 / piますが、テストしている最高周波数は1 / 2 / pi...次の単一の変更を試してください。

freqs = linspace(0.01, 3, 3000)

そして今、出力は期待されています:

1/pi = 0.318309886184
Frequency = 0.318311478264

ただし、periodogramに対してプロットするとfreqs / 2 / np.pi、グラフは次のようになります。

ここに画像の説明を入力してください

maxしたがって、より複雑な信号の場合、高調波があなたを騙す可能性があるため、ピリオドグラムのを探して支配的な周波数を見つけるだけに頼ることはできません。

于 2013-01-25T17:23:48.743 に答える