scipy库 signal 导入_scipy.signal.lombscargle的使用

For class, we are trying to prove a simple example of the Lomb-Scargle Periodogram using the embedded package in scipy. There is little documentation on how to use this function, and I have not been able to find any help online. When I run the code, I get a value of ~6.3 for the main peak of the periodogram, as opposed to the expected ~23.3. The data that we are pulling from is a simple .dat file with lists of numbers. Here is the code, any ideas on what's happening?

import scipy as sp

import math as m

import numpy as np

from scipy.signal import lombscargle

import pylab as plt

from numpy import shape

x=[]

y=[]

nout = 10000

file=open("hv878_1945.dat", 'r')

for pair in file:

xandy=pair.split()

x.append(float(xandy[0]))

y.append(float(xandy[2]))

x=np.asarray(x)

y=np.asarray(y)

f = np.linspace(0.1, 50, nout)

periodogram=sp.signal.spectral.lombscargle(x,y,f)

normval = x.shape[0]

plt.plot(f, np.sqrt(4*(periodogram/normval)))

plt.show()

Any help would be appreciated!

解决方案

It would appear you have the same problem described here. Basically you are not testing the frequency of interest and did not take into account that the frequencies passed to lombscargle are angular frequencies.

If you change your list of tested frequencies, you'll see there's a peak near the angular frequency of ~138, which corresponds to a frequency of 22.28 (=angular frequency /2/pi):

import numpy as np

import matplotlib.pyplot as plt

from scipy.signal import lombscargle

plt.ion()

data = np.loadtxt('hv878_1945.dat')

ang_freq = np.linspace(.1, 25*6, 10000)

periodogram = lombscargle(data[:,0], data[:,2], ang_freq)

print(ang_freq[np.argmax(periodogram)]/2/np.pi)

plt.plot(ang_freq, periodogram)

I will make the same remark that Jaime made in the other SO question:

you cannot rely on just looking for the max of periodogram to find the dominant frequency, because the harmonics may fool you.