Least-squares spectral analysis

In the Vaníček method, a discrete data set is approximated by a weighted sum of sinusoids of progressively determined frequencies, using a standard linear regression, or least-squares fit.[16] The frequencies are chosen using a method similar to Barning's, but going further in optimizing the choice of each successive new frequency by picking the frequency that minimizes the residual after least-squares fitting (equivalent to the fitting technique now known as matching pursuit with pre-backfitting[11]). The number of sinusoids must be less than or equal to the number of data samples (counting sines and cosines of the same frequency as separate sinusoids).

A data vector Φ is represented as a weighted sum of sinusoidal basis functions, tabulated in a matrix A by evaluating each function at the sample times, with weight vector x:

${\displaystyle \phi \approx {\textbf {A}}x}$

where the weight vector x is chosen to minimize the sum of squared errors in approximating Φ. The solution for x is closed-form, using standard linear regression:[17]

${\displaystyle x=({\textbf {A}}^{\mathrm {T} }{\textbf {A}})^{-1}{\textbf {A}}^{\mathrm {T} }\phi .}$

Here the matrix A can be based on any set of functions that are mutually independent (not necessarily orthogonal) when evaluated at the sample times; for spectral analysis, the functions used are typically sines and cosines evenly distributed over the frequency range of interest. If too many frequencies are chosen in a too-narrow frequency range, the functions will not be sufficiently independent, the matrix will be badly conditioned, and the resulting spectrum will not be meaningful.[17]

When the basis functions in A are orthogonal (that is, not correlated, meaning the columns have zero pair-wise dot products), the matrix A^TA is a diagonal matrix; when the columns all have the same power (sum of squares of elements), then that matrix is an identity matrix times a constant, so the inversion is trivial. The latter is the case when the sample times are equally spaced and the sinusoids are chosen to be sines and cosines equally spaced in pairs on the frequency interval 0 to a half cycle per sample (spaced by 1/N cycle per sample, omitting the sine phases at 0 and maximum frequency where they are identically zero). This particular case is known as the discrete Fourier transform, slightly rewritten in terms of real data and coefficients.[17]

${\displaystyle x={\textbf {A}}^{\mathrm {T} }\phi }$     (DFT case for N equally spaced samples and frequencies, within a scalar factor)

Trying to lower the computational burden of the Vaníček method in 1976 [7] (no longer an issue), Lomb proposed using the above simplification in general, except for pair-wise correlations between sine and cosine bases of the same frequency, since the correlations between pairs of sinusoids are often small, at least when they are not too closely spaced. This is essentially the traditional periodogram formulation, but now adopted for use with unevenly spaced samples. The vector x is a good estimate of an underlying spectrum, but since correlations are ignored, Ax is no longer a good approximation to the signal, and the method is no longer a least-squares method – yet it has continued to be referred to as such.

Rather than just taking dot products of the data with sine and cosine waveforms directly, Scargle modified the standard periodogram formula to first find a time delay ${\displaystyle \tau }$ such that this pair of sinusoids would be mutually orthogonal at sample times ${\displaystyle t_{j}}$, and also adjusted for the potentially unequal powers of these two basis functions, to obtain a better estimate of the power at a frequency,[3][8] which made his modified periodogram method exactly equivalent to Lomb's method. The time delay ${\displaystyle \tau }$ is defined by the formula

${\displaystyle \tan {2\omega \tau }={\frac {\sum _{j}\sin 2\omega t_{j}}{\sum _{j}\cos 2\omega t_{j}}}.}$

The periodogram at frequency ${\displaystyle \omega }$ is then estimated as:

${\displaystyle P_{x}(\omega )={\frac {1}{2}}\left({\frac {\left[\sum _{j}X_{j}\cos \omega (t_{j}-\tau )\right]^{2}}{\sum _{j}\cos ^{2}\omega (t_{j}-\tau )}}+{\frac {\left[\sum _{j}X_{j}\sin \omega (t_{j}-\tau )\right]^{2}}{\sum _{j}\sin ^{2}\omega (t_{j}-\tau )}}\right)}$

which Scargle reports then has the same statistical distribution as the periodogram in the evenly sampled case.[8]

At any individual frequency ${\displaystyle \omega }$, this method gives the same power as does a least-squares fit to sinusoids of that frequency, of the form

${\displaystyle \phi (t)=A\sin \omega t+B\cos \omega t.}$[18]

The standard Lomb–Scargle periodogram is valid for a model with zero mean. Commonly, this is approximated by subtracting the mean of the data before calculating the periodogram. However, this is an inaccurate assumption when the mean of the model (the fitted sinusoids) is non-zero. The generalized Lomb–Scargle periodogram removes this assumption, and explicitly solves for the mean. In this case, the function fitted is

${\displaystyle \phi (t)=A\sin \omega t+B\cos \omega t+C.}$[19]

The generalized Lomb–Scargle periodogram has also been referred to as a floating mean periodogram.[20]

Michael Korenberg of Queen's University in Kingston, Ontario, developed a method for choosing a sparse set of components from an over-complete set, such as sinusoidal components for spectral analysis, called fast orthogonal search (FOS). Mathematically, FOS uses a slightly modified Cholesky decomposition in a mean-square error reduction (MSER) process, implemented as a sparse matrix inversion.[15][21] As with the other LSSA methods, FOS avoids the major shortcoming of discrete Fourier analysis, and can achieve highly accurate identifications of embedded periodicities and excels with unequally spaced data; the fast orthogonal search method has also been applied to other problems such as nonlinear system identification.

Palmer has developed a method for finding the best-fit function to any chosen number of harmonics, allowing more freedom to find non-sinusoidal harmonic functions.[22] This method is a fast technique (FFT-based) for doing weighted least-squares analysis on arbitrarily spaced data with non-uniform standard errors. Source code that implements this technique is available.[23] Because data are often not sampled at uniformly spaced discrete times, this method "grids" the data by sparsely filling a time series array at the sample times. All intervening grid points receive zero statistical weight, equivalent to having infinite error bars at times between samples.