Correct way to use scipy.signal.spectral.lombscargle
numpy, python, scientific-computing, scipy, signal-processing
Solution
You are rightfully expecting the peak to show up at `1 / pi`, but the highest frequency you are testing is `1 / 2 / pi`... Try the following single change :
freqs = linspace(0.01, 3, 3000)
and now the output is the expected:
1/pi = 0.318309886184
Frequency = 0.318311478264
Note, though, that if you plot `periodogram` against `freqs / 2 / np.pi`, the graph looks like this:
So for a more complicated signal, you cannot rely on just looking for the `max` of periodogram to find the dominant frequency, because the harmonics may fool you.
Problem
I'm refering to the following post : Using scipy.signal.spectral.lombscargle for period discovery I realize the answer given correct for certain case. Frequency for sin(x), which is 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) ``` The following is printed. Is fine. I guess. The reason we divide the `lombscargle` result with `2pi` is that, we need to convert radian to frequency. (f = radian / 2pi) ``` 1/2pi = 0.159154943092 Frequency = 0.159154943092 ``` However, thing seems goes wrong for the following case. Frequency for sin(2x), which is 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) ``` The following is being printed. ``` 1/pi = 0.318309886184 Frequency = 0.0780862900972 ``` Seem incorrect. Any step that I had missed out?