IvS-KULeuven / IvSPythonRepository

Python Repository of the Institute of Astronomy @ KU Leuven
14 stars 22 forks source link

scargle and fasper returning insensible results #10

Open JorisDeRidder opened 5 years ago

JorisDeRidder commented 5 years ago

time = np.arange(0.,30.,1.) signal = np.sin(2.*np.pi*2.15*time) freq, spec = scargle(time, signal, f0=0., fn=20., df=0.0001, norm='amplitude') print(freq)

gives

[3.44827586e-04 4.44827586e-04 5.44827586e-04 ... 4.99744828e-01 4.99844828e-01 4.99944828e-01]

Why is the frequency array only going to 0.5 and not to 20. as specified? I have the same problem with fasper().

I'm using Python 3.6. Compilation of the Fortran routines after installing the repository went fine without any error messages.

mike-ivs commented 5 years ago

Just looked into this and it appears to be intentional behaviour for both scargle and fasper.

Using your variables and input in conjunction with the 'pergrams.chec_input()' function gave the following results:

>>> perg.check_input(time, signal, f0=0., fn=20., df=0.0001, norm='amplitude') OK: inputs are arrays
OK: inputs are 1D and have same length
OK: no infs or nans
OK: time array is sorted
No inconsistencies found or inconsistencies are fixed
Default Nyquist frequency: 0.5
Final frequency 'fn' is larger than the Nyquist frequency

(array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27., 28., 29.]), array([ 0.00000000e+00, 8.09016994e-01, 9.51056516e-01, 3.09016994e-01, -5.87785252e-01, -1.00000000e+00, -5.87785252e-01, 3.09016994e-01, 9.51056516e-01, 8.09016994e-01, 8.81869492e-15, -8.09016994e-01, -9.51056516e-01, -3.09016994e-01, 5.87785252e-01, 1.00000000e+00, 5.87785252e-01, -3.09016994e-01, -9.51056516e-01, -8.09016994e-01, -1.76373898e-14, 8.09016994e-01, 9.51056516e-01, 3.09016994e-01, -5.87785252e-01, -1.00000000e+00, -5.87785252e-01, 3.09016994e-01, 9.51056516e-01, 8.09016994e-01]))

According to the Pergrams.py documentation the output is capped at the Nyquist frequency, in this case a cap of 0.5. You should be able to change this if you wish.

Section 2. Nyquist frequency

The periodogram functions are written such that they never exceed the value of the Nyquist frequency. This behaviour can be changed.

By default, the Nyquist frequency is defined as half of the inverse of the smallest time step in the data. That is, the C{nyq_stat} variable is set to the function C{np.min}. If you prefer the Nyquist to be defined as the median, set the C{nyq_stat} variable for any periodogram to C{np.median}. If you want a more complex or self-defined function, that is also acceptable. If you give a number to C{nyq_stat}, nothing will be computed but that value will be considered the nyquist frequency.

JorisDeRidder commented 5 years ago

The Scargle documentation mentions that (f0, fn, df) are the start, stop, and step frequency. So fn is only relevant when it's below the Nyquist frequency?

Scargle, Fasper, and the other functions were originally designed to work for non-equidistant time series where the Nyquist frequency may be ill-defined. How is the stopping frequency then determined?

JorisDeRidder commented 5 years ago

For fasper() I just found the answer on my previous question by looking into the Python code: the average time sampling is taken. This is pretty horrible. The routines were probably written for analysis of CoRoT and Kepler data which is almost equidistant. For ground-based, Gaia, or LSST data this would not make sense.

For scargle() I don't see any reason in the code why fn would be capped to the Nyquist frequency.

mike-ivs commented 5 years ago

Hmm, I also can't find a reason for the Scargle 'fn=Nyquist frequency' in the Python code, nor in the Fortran code pyscargle.f (though I am neither an expert in fortran nor frequency analysis, maybe I missed something).

Perhaps the ivs.timeseries is an area that can be improved by using Astropy...?