Correct way to use scipy.signal.spectral.lombscargle

Question

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?

Crellen

Answer 1

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.

Greyhen

Thanks for the info. By the way, I generate more frequencies between 0.01 to 3 freqs = np.linspace(0.01, 3, 6000). I expect the result will be much more closer to 1/pi = 0.318309886184. However, when I run, that is not the case. I get Frequency = 0.417415475576. Is there any rule of thumb to follow? Thanks – Crellen

It turns out that I do not provide sampling points which is closed enough to each others : time = np.linspace(1, 1000, 100) – Crellen

It also helps to spread the frequencies geometrically if the span is a couple of orders of magnitude, e.g 0.25 Hz to 10kHz: ang_freq = np.geomspace(start_af, end_af, 10000) – Afc

The testing signal has been created with evenly spaced samples while lambscargle has been developed to manage unevenly spaced samples. If you look to scipy documentation you can find that by sampling the same signal with random periods the harmonics will disappear. You may consider FFT if you are dealing with evenly-spaced sampled signals. – Vacla