Spectral Analysis of Irregularly Sampled Data

Suppose you have a set of astronomical time-series data which is not uniformly sampled in time. If you expect this data will exhibit some periodicities, how do you do a spectral analysis? If the data is sampled irregularly, Fourier cannot be used directly. Another approach might involve interpolation or fitting the data, and then resampling it uniformly to allow Fourier analysis. Here we will show a simple spectral analysis method which avoids interpolation and resampling.

For each frequency ω of interest, fit α and β in α sin(ω t) + β cos(ω t) to the data by linear least squares. The spectral energy ℰ[ω ] is then the square amplitude integrated over one period, T = 2π/ω:

(((∫_0^(((2 π)/ω)))(((α Sin[ω t] + β Cos[ω t]))^2)
               \[DifferentialD]t)/((∫_0
              ^(((2 π)/ω))) \[DifferentialD]t))//Simplify

To demonstrate this method, we will use the data set pulsar.dat, which consists of timing residuals for the pulsar PSR1257+12 spanning a three-year period (kindly supplied by Alex Wolszczan). The data represent post-fit residuals after fitting a standard timing model assuming a single pulsar. Click here to change to the directory that contains the pulsar data file used in this section.

data = ReadList["pulsar.dat", Number, RecordLists -> True];

Here is the first datum:

Take[data, 1]

Column 2 is the epoch (in decimal years) and column 3 is the timing residual (in milliseconds).

data = Transpose[data][[{2, 3}]];

Here is the timing residual variation in fractional years since 1990:

data = data - {1990, 0} // Transpose;

{(Subscript[t, min]), (Subscript[t, max])} = First/@{First[data], Last[data]}

lp = 
      ListPlot[data, 
        PlotStyle -> {RGBColor[1, 0, 0], AbsolutePointSize[4]}, 
        AspectRatio -> 1/5, AxesOrigin -> {(Subscript[t, min]), 0}];

It is straightforward to compute the spectral energy,

\(ℰ[ω_] := 
      Fit[data, {Sin[ω\ t], Cos[ω\ t]}, t] /. 
		
        α_\ Sin[_] + β_\ Cos[_] :> 
          \((α\^2 + β\^2)\)/2\)

and then plot it as a function of ω.

Plot[ℰ[ω], {ω, 1, 50}, 
      PlotRange -> All];

This graph reveals peaks at a number of frequencies. The two principal frequencies are

(Subscript[ω, 1]) = 
      ω /. 
        Last[FindMinimum[
            -ℰ[ω], {ω, 23, 25}]]

and

(Subscript[ω, 2]) = 
      ω /. 
        Last[FindMinimum[
            -ℰ[ω], {ω, 34, 35}]]

These two frequencies correspond to “planets” orbiting the pulsar with orbital periods of approximately

(Subscript[T, 1]) = ((2  π)/(Subscript[ω, 1])) 365.24 Days

and

(Subscript[T, 2]) = ((2  π)/(Subscript[ω, 2])) 365.24 Days

These values can be compared with those at http://darc.obspm.fr/

Since we now have good estimates for the orbital periods, it is easy to do a nonlinear fit of a double planet model.

\(\(basis[α_, β_, θ_, γ_, δ_, ϕ_, 
        ϵ_] = 
	
      α\ Sin[θ\ t] + β\ Cos[θ\ t] + 
        γ\ Sin[ϕ\ t] + δ\ Cos[ϕ\ t] + ϵ; 
    \)\)

The approximate periods provide good starting values for the nonlinear fitting routine.

Needs["Statistics`NonlinearFit`"]

NonlinearFit[data, 
      basis[α, β, θ, γ, δ, ϕ, 
        ϵ], 
      t, {{α, 1}, {β, 1}, {γ, 1}, {δ, 0}, {
          θ, (Subscript[ω, 1])}, {ϕ, (Subscript[ω, 2])}, {ϵ, 0}}]

This plot shows the equality of the nonlinear fit:

Plot[%, {t, (Subscript[t, min]), (Subscript[t, max])}, Epilog -> lp[[1]]];

For a visualization of the planet data, see http://www.empire.net/~whatmoug/Extrasolar/home_PSR1257.html.

References

A. Wolszczan and D. Frail. A planetary system around the millisecond pulsar PSR1257+12, Nature 255:145, 1992.

A. Wolszczan. Confirmation of Earth-mass planets orbiting the millisecond pulsar PSR B1257+12, Science 264:538, 1994.

Pulsar planets: Department of Astronomy and Astrophysics at Penn State University.