Algorithms for Spectral Analysis of Irregularly Sampled Time Series

In this paper, we present a spectral analysis method based upon least square approximation. Our method deals with nonuniform sampling. It provides meaningful phase information that varies in a predictable way as the samples are shifted in time. We compare least square approximations of real and complex series, analyze their properties for sample count towards infinity as well as estimator behaviour, and show the equivalence to the discrete Fourier transform applied onto uniformly sampled data as a special case. We propose a way to deal with the undesirable side effects of nonuniform sampling in the presence of constant offsets. By using weighted least square approximation, we introduce an analogue to the Morlet wavelet transform for nonuniformly sampled data. Asymptotically fast divide-and-conquer schemes for the computation of the variants of the proposed method are presented. The usefulness is demonstrated in some relevant applications.


Introduction
Despite the vast number of methods for time series analysis, there is still a lack in manageable methods dealing with nonuniformly sampled data. In this paper, we present a spectral analysis method based upon least square approximation, applied to different types of data. Our method is based upon assumptions similar to those of the well-known Lomb-Scargle periodogram (Lomb 1976;Scargle 1976;Press, Teukolsky, Vetterling, and Flannery 1992;Bretthorst 2000), and deals with nonuniform sampling. Beyond that, it has the advantage of a meaningful phase that varies in a predictable way as the samples are shifted in time, while the Lomb-Scargle periodogram was designed to compute the power spectrum. We study estimator properties and the convergence of our results for the sample counts going towards infinity. We compare least square approximations of real and complex series and show the equivalence to the discrete Fourier transform applied onto uniformly sampled data as a special case. We analyze how constant offsets affect the obtained spectral coefficients and propose a way to deal with the undesirable side effects of nonuniform sampling.
We then reanalyze the approach with focus on weighted least square approximation and introduce an analogue to the Morlet wavelet transform (Chui 1992;Daubechies 1987;Holschneider 1995) for nonuniformly sampled data. Next, we construct a family of asymptotically fast divide-and-conquer schemes for the evaluation of the presented analysis methods. We demonstrate the functioning of our methods by applying them onto the paleoclimatic CO 2 and Deuterium concentrations covering a time span of approx. 420000 years that were obtained from the ice core (Petit JR et al. 1999) extracted from the ice cover of Lake Vostok, Antarctica. The occupation with the nonuniformly sampled data presented here were the main reason for the work that lead to this publication. The paper and the data are available from http://www.nature.com/cgi-taf/DynaPage.taf?file=/nature/journal/ v399/n6735/abs/399429a0_fs.html. Finally, we apply the methods to tick data obtained on one day at the London stock exchange. Since the trade occurs at irregular times, one also has to deal with nonuniformly sampled data.
The coefficients c ω and s ω can be obtained by multiplying Eq. 1 with the Moore-Penrose pseudo-inverse (Penrose 1955), (c ω s ω ) = x Ω T t,ω Ω t,ω Ω T (c ω s ω ) = x k cos ωt k sin ωt k T cos 2 ωt k cos ωt k sin ωt k cos ωt k sin ωt k sin 2 ωt k −1 = 2 x k cos ωt k sin ωt k T n − cos 2ωt k − sin 2ωt k − sin 2ωt k n + cos 2ωt k n 2 − ( cos 2ωt k ) 2 − ( sin 2ωt k ) 2 ( 2) and, based on that, we introduce P ω = c 2 ω + s 2 ω , a value that will be used later in considerations about power spectra. Lomb (1976) made similar assumptions about the nature of the analyzed data and introduced a time shift θ of the series that makes the vectors {cos (ω (t k − θ))} k∈K and {sin (ω (t k − θ))} k∈K orthogonal, dependent on ω: The representation using arctan is the common way to express θ. Note that we silently assume it to be equivalent to the representation using arg, i.e. taking into account the signs of numerator and denominator separately and dealing with 0 in the denominator.
Applying Eq. 2 onto the original time series with changed sample times t k − θ, k ∈ K, yields The absolute square value is which, apart from the squared denominators, is equivalent to the Lomb periodogram,  Figure 1: The solid line shows the absolute square result of Eq. 3 applied to sin(4.8t) + X 0,2 , sampled at 400 values of t that are uniformly randomly distributed over [0, 2π). X 0,2 denotes a normally distributed random variable with mean 0 and variance 2. The dashed line shows the Lomb periodogram applied to the same data, scaled with 1 182 . This factor has to be adapted for different sample counts.
Using the mean x and the variance σ 2 of {x k }, the Lomb normalized periodogram is defined as The replacement of the x k with x k − x is one way to deal with time series with constant offset, for which we propose an alternative solution below. This replacement is also commonly used for the unnormalized Lomb periodogram.
The amplitudes of both versions of the Lomb periodogram appear to depend on the sampling density, as can be seen in in Figure 2.
In the case of a complex time series ξ = {ξ k } k∈K and an approximation by complex functions e iωt , the result is even simpler. We use the matrix Ψ t,ω , Ψ t,ω = e iωt 0 e iωt 1 . . . e iωt n−1 to formulate the expression   analogously to Eq. 1 with a complex coefficient ζ ω and a vector of complex zero average error terms.
The pseudo-inverse solution yields in this case Using z = 1 2 (z + z) and z = 1 2i (z − z), we rewrite the real approximation of Eq. 2 in terms of the complex approximation. Here, ζ ω denotes the complex approximation of the real series {x k } k∈K , ι 2ω denotes the complex approximation at frequency 2ω of the unit time series {1} k∈K over {t k } k∈K , and c and c mean real and imaginary parts of a complex number c, respectively. We then obtain In the uniform case of t k = k∆t, ∆t = 0, k ∈ K = {0, . . . , n − 1} and ω = lπ n∆t , l ∈ K \ {0}, the sums cos 2ωt k and sin 2ωt k vanish in the real approximation given by Eq.
(2), and the result is equivalent to the complex conjugate of one spectral coefficient of the discrete Fourier transform. For ω = 0, the average c ω = x is obtained. Each term of the sum in Eq. 6 can be written as a matrix multiplication, Approximating such a uniform real-valued series with the complex approximation of Eq. 7 leaves the imaginary part undetermined, whereas the strictly real approximation of Eq. 2 can be interpreted as a complex equation with imaginary part set to zero, x k = 0. Eq. 2 can be used to approximate both time series x k and x k of Eq. 9, whereby x k = 0, as mentioned before, which obviously leads to c ω = s ω = 0. The complex approach of Eq. 9 chooses the average between the results from Eq. 2 and c ω = s ω = 0, which is 1 2 (c ω s ω ). This explains the factor 2 in Eq. 2.
In the general case of nonuniformly sampled time series, one usually estimates power spectra by means of periodograms, as for example the frequently used Lomb-Scargle periodogram (Lomb 1976;Scargle 1976;Press et al. 1992). When however phase information is desired, e.g. for the reconstruction of time series out of spectral coefficients and for comparison of signals, Eq. (2) does the job. In the sequel, we also use it for the extension of wavelet analysis onto nonuniform time series.
The amplitude of our result is invariant with respect to time shift, whereas the phase is covariant, as can be verified easily. By introducing a time shift ∆t, we have Ω t+∆t,ω = cos ω (t 0 + ∆t) cos ω (t 1 + ∆t) · · · cos ω (t n−1 + ∆t) sin ω (t 0 + ∆t) sin ω (t 1 + ∆t) · · · sin ω (t n−1 + ∆t) The pseudo-inverse of Ω t+∆t,ω becomes then −1 R T ω∆t and the coefficients after time shift Obviously, for the complex approximation (Eq. 7) the time shift yields Both the real and the complex approximation's amplitudes are thus shift invariant, and a meaningful phase is obtained. The well-known advantage of periodograms applied to nonuniformly sampled time series is also given here, namely their insensitivity to aliasing phenomena when the differences between the t i obey a normal distribution.
We note in passing that the reconstruction of a signal from the spectral analysis of a time series requires a full matrix inverse as a consequence of the nonuniform sampling.

Convergence
We base our following considerations on a fixed real infinite sequence {t k } k∈N 0 with a generalized function describing the empirical density that the t k obey, and a real function x(t) where x(t)ϕ(t) is absolutely integrable. For sequences with noncontinuous empirical density, e.g. t k = k mod m for a fixed m ∈ N, ϕ(t) may become infinite for some values of t, a case which is handled by the theory of generalized functions or distributions. In most cases it suffices however to consider the limits of ordinary functions, e.g.
for the previously mentioned series.
Under these circumstances, we can write We next consider a complex function ξ(t), where ξ(t)ϕ(t) is absolutely integrable, and a measurement error ρ k obeying an independent complex symmetric distribution with zero mean. The values ξ (t k ) are identified with the ξ k . The complex approximation of Eq. 7 applied to ξ(t k ) + ρ k becomes then where (F t f (t)) (ω) denotes the Fourier transform of f (t) with respect to the variable t at frequency ω. We note that when F t ξ(t) exists, the complex approximation converges to F t ξ(t) * F t ϕ(t). The nonconstant empirical sample density ϕ(t) thus makes the complex approximation result a biased spectral estimator. The bias converges to 0 as F t ϕ(t) approaches δ 0,t , i.e. components of ϕ(t) with non-zero frequency vanish. We note that which means that the complex approximation is a consistent estimator.
The rewritten real approximation in Eq. 8 becomes We consider ω = 0. As F t ϕ(t) approaches δ 0,t , (F t ϕ(t)) (2ω) vanishes, leaving a result that, except for the factor 2 explained above, is equivalent to the conjugate of the complex approximation. Under the assumed circumstances, the real approximation thus converges to something different from the complex approximation. If ϕ is symmetric, then F t ϕ(t) is real, and the real and imaginary parts of F t x(t)ϕ(t) are scaled individually. In the general case, the real approximation converges to a complex linear combination of the complex approximation's limit and its complex conjugate.
Since the real approximation is based on the complex approximation, it is, for F t ϕ(t) → δ 0,t , an unbiased estimator. Assuming a real symmetric distribution for the error term ρ k in Eq. 12, we also obtain zero variance for ζ ω which means that the real approximation is a consistent estimator as well.
Concluding, we can state that when the Fourier transform of the randomly sampled signal exists, i.e. the signal is a superposition of sinusoids, both the complex and the real approximation converge and provide consistent and for F t ϕ(t) → δ 0,t unbiased estimators.
In order to study the convergence of the Lomb periodogram, we first introduce Rewriting Eq. 5 and applying Eq. 11, we obtain

Wavelet transforms
We now combine the approximations derived above with a weighted least square approximation to formulate a windowed transform. In our case, a weighted least square term that should be minimized with respect to (c ω s ω ) reads where β 2 k ≡ w(t k ) is an appropriate weight function. We obtain the following pseudo-inverse solution for Eq. 14: The complex version of this is again much simpler, For continuous functions as well as for evenly sampled data, the Morlet wavelet transform is an instance of a weighted Fourier spectrum, whereby the window width scales with ω. It uses the Gaussian e −(σωt) 2 with an additional scaling parameter σ as weight function.
We may now define Morlet-like wavelet transforms for unevenly sampled data by using the windowed approximations from Eq. (15) and Eq. (16). The convolution 1 a ψ * ( b−τ a )f (τ )dτ uses the mother wavelet ψ(t), which we choose to be ψ(t) = w(σt)e −iωt in accordance with our approach. It remains to introduce a time shift. The complex case then becomes, In order to accelerate the computation, we use a window function with compact support in time that does not employ transcendental functions, This window function has a shape that is very similar to the popular Hanning window (see Figure 5) Its Fourier transform isŵ (ω) = 12 2 − 2 cos(ω) − ω sin(ω) ω 4 , and it produces comparable sideband ripples (see Figure 5). The first sideband ripple is even lower than the Hanning window's one and not substantially worse than the sideband falloff of the Gaussian's spectrum. All three power spectra are shown in Figure 5. An appropriate choice for the value of σ is 0.05. This means that 2 σ radians or 1 πσ periods of e −iωt fit into the whole support of the weight function.
In order to estimate the spectral phase differences between two complex time series ξ k , t k , k ∈ 0, . . . , K − 1 and η l , u l , l ∈ 0, . . . , L − 1, we use the complex conjugate product of their approximations according to Eq. 17, Since e −iωt k ξ k and e iωu l η l do not vary with t, they can be precomputed, which accelerates the computation in comparison to Eq. 17 when the result is computed with increasingly fine time resolution.

Fast approximation
We use the power series approximation of e it to formulate a divide-and-conquer scheme for the evaluation of Eq. 7 consisting of 4 fundamental steps described below. A fast algorithm for the real approximation that makes use of the complex approximation is given thereafter.

Subdivision
In a first step, the index set K = {0, . . . , n − 1} is subdivided into subsets K h with associated τ h such that for a given highest analysis frequency ω max , when e −iωmax(t k −τ h ) is approximated by a power series, the error is sufficiently small.

Precomputation
Next, for each K h , the vector s h of the P + 1 first terms of the power series expansion of The factor µ plays a role in the improvement of numeric accuracy. For now, we can assume µ = 1.

Summation for a given ω
Using the s h , we now compute the spectral coefficients for 1 2 ω max < ω ≤ ω max .
Using the fact that the τ h are evenly spaced, the e −iωτ h do not have to be calculated explicitly but can be obtained from e −iωτ 0 by repetitive multiplication of e −2iω∆τ .
Merging of s 2h and s 2h+1 into s h For the upper limit of the next frequency interval ω max := 1 2 ω max , the power series approximation has the required exactness over an interval of double length. It follows that the s 2h and s 2h+1 can be combined to s h , using the offset τ h := 1 2 (τ 2h + τ 2h+1 ) and size ∆τ := 2∆τ belonging to the index set K h with −∆τ ≤ t k − τ h < ∆τ for all k ∈ K h . The index h now runs from 0 to H := H 2 . The vectors of power series expansion terms s 2h and s 2h+1 are combined into s h using reexpansion: With the introduction of ∆τ = τ 2h+1 − τ h = − (τ 2h − τ h ), this can be further simplified: This process is repeated until either ω max passes the desired minimal analysis frequency ω 0 , or the merging of adjacent subsets leaves just one subset behind (which means that at frequencies below the ω (d) max where this occurs, the time range of the time series covers merely a fraction of a full period of e −iωt , which is of questionable value, except for ω = 0). When the number of subsets in one level is odd, the last remaining vector s

Real approximation
Based on the scheme for complex approximation and the complex representation of the real approximation in Eq. 8, a computation scheme for the real approximation can be constructed, As in Eq. 7, ζ ω denotes the complex approximation result for the considered time series, whereas ι 2ω denotes the complex approximation at frequency 2ω of the unit time series (1 . . . 1) over (t 0 . . . t n−1 ) that can be computed in a separate subdivision scheme with a few obvious simplifications stemming from the fact that the sample values equal 1.
We note that in the real case, the s h,p are real when p is even and imaginary otherwise. In Eq. 23, it can be seen that most of the computation steps for the merging process can be performed with purely real arithmetic, which further speeds up the computation.

Complexity considerations
We analyze the application of the fast scheme onto logarithmic and linear ranges of frequencies.
In the logarithmic case, the algorithm provides the highest acceleration compared to the straighforward implementation of 7. Assuming that n oct octaves of n coeff spectral coefficients per octave shall be computed, the following steps, each shown with their complexity, have to be performed: • Subdivision and precomputation: O(nP ), • Summation for n coeff spectral coefficients in each of the n oct octaves with ω • Merging of the s h and s h+1 : Summing everything up and applying the sum of geometric series, we obtain a total complexity of When a linear frequency range ω l = lωmax L , l ∈ {0, . . . , L} is desired, the number of coefficients that can be computed with subdivision level d is For ω = 0, the coefficient is computed using the 0 th term of the s h . Summation of the above and substitution of geometric series yields

Weighted complex approximation
Here, we use the Hanning window from Eq. 18. We build upon subdivision and precomputation steps similar to the ones described above, with where Q can be chosen smaller than P .
We then rewrite Eq. 17 as follows, using the abbreviations ω − = ω(1 − σ) and ω + = ω(1 + σ). K(t) denotes the range k a . . . k b of indices of samples that lie within the translated and scaled Hanning window, i.e. σω (t ka − t) < π and σω (t k b − t) < π. H(t) denotes the range of indices of subdivision intervals that cover the Hanning window. The factors e −iω − (τ h −t) , e −iω(τ h −t) , e −iω + (τ h −t) and e iσω(τ h −t) can be multiplied up iteratively from their respective start value for τ ha . The pairwise merging of the s h and j h is performed similarly to the complex approximation. It has to be noted that the process of covering the periodic extension of the Hanning window used with precomputed sample ranges introduces small errors since the sample range boundaries almost never coincide with the window's boundary.

Approximation for time series with constant offset
When analyzing time series with a constant offset, e.g. x + d(1 . . . 1), d ∈ R, using Eq. 2, we have d(1 . . . 1)Ω T t Ω t Ω T t −1 = 0 in the general case, due to the nonorthogonality of the row vectors of Ω t and the unit vector. This leads to increasingly wrong spectral coefficients for low frequencies.
We therefore compute (c ω s ω d ω )   cos ωt 0 cos ωt 1 . . . cos ωt n−1 sin ωt 0 sin ωt 1 . . . sin ωt n−1 1 1 by means of the pseudoinverse, n + cos 2ωt k sin 2ωt k 2 cos ωt k sin 2ωt k n − cos 2ωt k 2 sin ωt k 2 cos ωt k 2 sin ωt k 2n  Figure 6: The absolute value from Eq. 25 applied to the Vostok CO 2 series (solid) compared with the absolute value from Eq. 2 applied to the CO 2 series (densely dashed) and to the CO 2 series with subtracted mean value (widely dashed, nearly coincident with the solid curve).
using the symbols ζ ω and ι ω as above to denote complex approximations for the time series and the unit series. The ζ ω , ι ω and ι 2ω can be computed in divide-and-conquer schemes, whereby the results for ι ω can be reused as ι 2ω in subsequent octaves.
When looking at Eqs. 1 and 24, it becomes obvious that at low frequencies, the approximation becomes problematic. In Figure 7 which displays results down to a frequency corresponding to a third of one period over the full sample time span of 422000 years, it can be seen that the low frequency coefficients start to diverge. As a rule of thumb, we note that approximation of time series over a time span of less than half a period becomes undefined.
Looking for periodic components in a signal, we did some promising experiments with the simultaneous approximation of sinusoids with m integer multiple frequencies of a given ω 0 , using the pseudoinverse solution of (d ω c ω s ω . . . c mω s mω ) It is possible to construct a divide-and-conquer scheme for these approximations as well. It becomes impracticable however to perform the matrix inversion algebraically.

Application samples
The software accompanying this paper is written as an R module, and was developed under Linux. Besides the functions written in native R in order to provide a comprehensive sample, the real and complex spectral transform algorithms and the weighted transform, as well as their non-accelerated counterparts have been implemented in C as a loadable library. We give an overview and a couple of examples using the implemented functions here; please refer to the online documentation in R for details concerning their application and further examples.
The R package is distributed as a compressed installation directory. Please use the shell command tar xzvf nuspectral.tgz; R CMD install nuspectral to install the package into the appropriate directories of your R installation. You may have to login as super-user first before performing the installation.
After having started R, please use library("nuspectral") to load the library. For the function names mentioned below and the data sets provided with the package, help pages are available that can be viewed using the R help command, e.g.

help("fastnureal")
The following functions, each operating on the time series defined by the abscissa vector X and the ordinate vector Y , are provided: • lombcoeff(X, Y, o) computes a coefficient of the Lomb periodogram for the given circular frequency.
• lombnormcoeff(X, Y, o) computes a coefficient of the Lomb normalized periodogram for the given circular frequency.
• nurealcoeff(X, Y, omegamax, ncoeff, noctave)nurealcoeff • nuwaveletcoeff(X, Y, t, o, wgt = cubicwgt, wgtrad = 1)computes a coefficient of the nonuniform wavelet transform in Eq. 17 at the given time and circular frequency. The weight function is by default the cubic polynomial weight function shown in Figure 5, and the radius beyond which the weight function is truncated to 0 is by default 1.
The following functions have been implemented in C.
• nucomplex(X, Y, omegamax, ncoeff, noctave)computes the logarithmically scaled complex spectral approximation of Eq. 7 for the given maximum circular frequency, number of coefficients and coefficients per octave.
• nureal(X, Y, omegamax, ncoeff, noctave)computes the logarithmically scaled real spectral approximation of Eq. 3 for the given maximum circular frequency, number of coefficients and coefficients per octave.
The data we will be using stem from • the Vostok Lake ice drilling core (Petit JR et al. 1999); the ice column extracted from a more than 3 km deep hole drilled into the antarctic ice over the huge frozen lake near the Russian Vostok research station provides invaluable information about climate, chemical composition of the atmosphere and other factors within the last 422000 years. The time series are all nonuniform, and almost none of the sample times of different series match.
• tick data from the 11. July 2001 at the London stock exchange. Given are the rates of all stocks sampled at the transaction times. Since the rates change upon trade transactions that occur more or less randomly, these time series are also sampled nonuniformly.
Distributed in R format with the package are the Vostok data series • co2 -the CO 2 content of the ice.
• ch4 -the Methane content of the ice.
• deut -the Deuterium content of the ice; the ratio of Deuterium to total Hydrogen is a proxy for the atmospheric temperature.
• o18 -the 18 O (a heavy Oxygen isotope) content of the air enclosed in the ice, a proxy for the solar irradiation.
• dust -the dust content of the ice.
The stock data example has been downloaded from the internet, and has not been included into the R package because of copyright issues.
The wavelet transform diagrams below haven't been generated with R. Some exemplary invocations of R functions that perform the equivalent transformations are presented, however, in the respective figure captions.

Summary
We have presented a way to compute spectra for nonuniformly sampled data based on least square approximation and the Moore-Penrose pseudoinverse, both for real and complex approximations. The presented methods are shift-covariant in time and provide phase information. We have analyzed the correspondence of our methods with the discrete Fourier transform and the impact of constant signal offsets onto the spectral coefficients, and studied estimator properties and the convergence of our results for the sample counts going towards infinity. We further have introduced weighted versions of these spectral transforms that provide wavelet transforms for nonuniformly sampled data.
As a most remarkable result, we have presented asymptotically fast algorithms for the computation of the real and complex spectral analysis, which makes our methods comparably as powerful as the fast Fourier transform. These fast and robust methods to evaluate the spectral signal in time series open up new possibilities for real-time applications of increasing importance in many fields, e.g. in medicine, physiology as well as video and audio processing. The methods have also been applied in an interactive work of art where the user can explore animated 3D renderings in image and sound of a spectral analysis of the Vostok data. Our methods also promise to be useful for reconstruction and interpolation of missing and distorted data. The potentiality of our methods has been hinted by the presented relevant applications.
Future work will include an improvement of the numerical basis of our fast algorithms, and the development of other fast algorithms based on the presented scheme.    Karlsruhe; Medialab, Madrid). The installation allows the visitor to navigate through the Vostok data material (here the CO 2 data set) in a way that is inspired by a sound sampler, using loops and variable playback speed. The data is analyzed using Eq. 17, and the absolute values are displayed as elevation and the phase as vector field. Further, an animated graphic representation of the equation acting on the data is shown. Beyond graphics, the original data and the spectral analysis result are used to directly control sound parameters of a software synthesizer, a process called audification.