Program Spectra

original version april 1998, this version january 2017.
Tabaré Gallardo (Dpto. Astronomía, Instituto de Física , Facultad de Ciencias)

[ascl:1701.003]



Power Spectrum

The search for periodicities in a time series is a classic and everyday problem in astronomy. The series can be equispaced (case of numerical integrations) or non-equispaced (case of observations). For equispaced series we have the powerful FFT method. However, this method has the important limitation that the power of the spectrum is evaluated in certain fixed frequencies (the Fourier frequencies) that may not be the ones that we are looking for. More precisely, it is a decomposition of the time series into an orthonormal basis defined by a finite number of frequencies, N/2 + 1, where N is the number of data. The FFT method is definitely useless for evaluating low frequencies. Although often used, the FFT should not be applied to unequal spaced observations because the fundamental hypothesis (equispaced data) is not fulfilled. Worse still, the FFT is usually used to evaluate the power in non-Fourier frequencies, arguing that it is working with the continuous transform. All this makes that the power evaluated by FFT is incorrect and therefore the spectral lines are badly identified. The first serious attempt to solve this problem related to unequally spaced data was made by Ferraz-Mello (1981, Astron., Journal 86 (4), 619) followed by Horne & Baliunas (1986, Astroph.J., 302, 757) and by Foster in a succession of memorable works in 1995 (Ast.J. 109 (4), 1889) and 1996 (Ast. J. 111 (1), 541 and Ast. J. 111 (1), 555).

There are many methods for obtaining spectra or periodograms. While the efficiency of methods depends on our purpose and the type of data we have, in general Foster's ideas are at the forefront. Here we present a method based on the ideas of Ferraz-Mello and Foster and whose details can be found in Gallardo & Ferraz-Mello (1997, Astr.J. 113 (2), 863). Basically it is about adjusting the data (t, x) to the model function

Y = C1 + C2*cos(2pi*f*t) + C3*sen(2pi*f*t)

For tentative values of the frequency f, the constants Ci and the Spectral Correlation Coefficient, CCE, (Ferraz-Mello 1981, Foster 1996a) are calculated. The CCE is a number less than 1 and is associated with the probability that the frequency f is real and not product of chance. The maxima in the graph CCE (f) indicate the presence of lines and the corresponding Ci define amplitude, phase and zero level. Logically, the power of the method lies in the definition of the CCE, for which we suggest to see the source program in Fortran language spectra.for and the references. Spectra also allows you to choose between using a Hanning window or a square (all weights are equal). The Hanning window eliminates the secondary lobes of the spectral lines but their width doubles. The square window maintains the natural width of the lines but the secondary lobes are notorious.


Instructions

  1. Compile "spectra.for" or use the compiled version provided.
  2. Save with the name "datos.dat" the file containing the time series in two columns (t,x) sorted chronologically. Put in the same directory the executable "spectra" and "datos.dat".
  3. When you run spectra, the program will ask you for which type of window to use. According to that the program will indicate the width of each spectral line in frequency units (cycles per unit time used).
  4. Then inticate the limits between which you want to calculate the spectrum. Here it is convenient to keep in mind that in the case of equispaced series with interval T between data it will not make sense to evaluate the spectrum at frequencies outside the range (0, 1/(2T)) or Nyquist frequency. The program will also ask the number of points that will have the plot CCE(f). The higher the number of points, the better the definition of the maxima in the spectrum, but if several thousand points are used the program can take an eternity.
  5. Once finalized, the results will be in the file spectra.dat with the format:
    Frequency, CCE, C1, C2, C3, amplitude.
    where amplitude = sqrt(C2^2 + C3^2).
The possible spectral lines will be defined by the CCE maxima. If you want the amplitude look for the amplitude corresponding to the maxima of CCE.

In this way the spectrum is obtained. If the lines are sufficiently separated (several times the width) there will be no leakage between them and the CCE maxima will define them. But if the lines are at distances smaller than the width then there will be coupling and its precise determination requires an extra work that is not contemplated in this program. One thing is obtaining the spectrum, CCE (f), and another is the precise determination of the lines. In order to separate close lines there are several methods, and we strongly recommend Foster (1995).

Low frequencies

Spectra is especially accurate for low frequencies, that is, when we try to determine periods by having observations covering only a fraction of it. For example, a pure sinusoid (no noise) is detected satisfactorily having only 1/6 of its period, something unthinkable in the traditional methods based on the Fourier Transform, where the spectrum is completely spurious in the low frequencies (by non orthogonality of the base on which the data are projected to calculate the power).

Download spectra

The file spectra.zip contains:
spectra.for (F77 source code)
spectra.exe (executable Windows)
spectra (executable Ubuntu)
datos.dat (example)
spectra.html (this webpage)