Multivariate Frequency and Amplitude Estimation for Unevenly Sampled Data Using by Extending the Lomb-Scargle Method

1. Introduction

In numerous signal processing applications, the characterization of a physical process via its power spectrum, or its amplitude and phase spectra, is of fundamental importance. Typically, the recorded signal is sampled at uniform intervals over a finite observation period and stored as discrete sampled data. Spectral analysis serves multiple purposes, including the identification of characteristic frequencies for process recognition, the selection or attenuation of specific spectral components, and the precise estimation of amplitude and phase at specified frequencies. Conventional approaches predominantly rely on the discrete Fourier transform (DFT) or its computationally efficient companion, the fast Fourier transform (FFT), both of which decompose the data into Fourier series coefficients, from which the power spectrum can be derived [1,2]. These transformations are inherently invertible, computationally efficient, and straightforward to implement. However, their applicability is constrained by the requirement of equidistant sampling, a limitation that poses challenges for data acquisition, in general. In practical applications, missing or corrupted data is often handled by substituting zero values, a common but problematic approach. In contrast to the intension, zero values provide some sort of information and do not mean “not aquired”. Subsequently, this procedure introduces a systematic error that leads to increased inaccuracies in frequency and amplitude estimation and ultimately affects the reliability of parameter estimation. While the DFT naturally extends to higher dimensions, facilitating applications such as image reconstruction in magnetic resonance imaging [3], image filtering [4], and higher-order spectral filtering in spatiotemporal domains [5], its performance deteriorates in the presence of non-uniform sampling. In ultrafast nuclear magnetic resonance (NMR) spectroscopy, non-uniform sampling enables substantial acceleration of chemical analysis [6,7], albeit at the cost of increased missing data, reduced sensitivity, and signal attenuation. These limitations are inherent to DFT-based and similar techniques, as they exhibit reduced accuracy when applied to incomplete datasets with data gaps [8]. Although alternative methods, such as the "non-uniform Fourier transform" [9,10,11], attempt to address these issues, they remain reliant on classical orthogonal mode decomposition (e. g., FFT, DFT, wavelet transform) and thus inherit systematic errors, as discussed in Section 2.2.

The application of techniques based on Fourier series to signals with non-equidistant sampling is difficult, since it usually requires zero padding (replacing missing values by zeros) on a grid and in the general case resampling or interpolation of the data to a regular grid [12,13]. The effect of a gap filling method is investigated extensively by Munteanu et al. [8]. The conclusion of the latter study is that without a corrective measure, errors such as amplitude minimization, which depends on the total gap density, cannot be avoided when FFT or DFT is applied. On the other hand, interpolation could be utilized to shift non-uniformly sampled data to a regular grid, so that the new interpolated data points are composed of the signal information and a projection of the accompanying noise. In the general case such distortions or noise are not band limited leading to biased estimates which consist of the local and aliased errors (noise, outliers, missing values, etc.) of the surrounding data.

Usually, astrophysical data is often affected by gaps, missing values and unevenly sampling. For example, when ground based radio telescopes are exploring areas in space there might exist time intervals in which the antenna is not pointing towards the object of interest due to the rotation of the Earth. The result is an incomplete data set. Non-uniform (random) sampling emerges for example if the irregular appearance of objects like Sunspots is measured as a binary quality depending on time and location. Section 5.3 provides a working example which analyzes frequencies and periods in latitude and time from the two dimensional data set of Sunspot observations. An other technical example for random sampling (with highly variable sampling frequency) is asynchronous data acquisition in large sensor networks, which are used in Smart Home, Industry 4.0 and automated driving, see Geneva et al. [14], Cadena et al. [15], Sudars [16]. Here the time series provide missing values and data gaps originating from time periods, in which strong noise prevents the measurement or simply the source of the signal to be measured is not in the range of the sensor.

The Lomb-Scargle method (LSM) was developed especially for ground based non-uniformly sampled one dimensional data, from which the amplitude spectrum is calculated [17,18]. The main advantage of this method is to directly estimate the spectrum without iterative optimization of trigonometric models as discussed in Section 2.1. A fast version for one dimensional signals is presented by Townsend [19] and Leroy [20].

An extension of LSM to two- or three-dimensional time series have not been presented yet. Especially for large multivariate data sets such a direct approach would be comparably faster than present iterative procedures as proposed in (Babu and Stoica [21], Chap. 9). The advantages of multivariate LSM will be demonstrated in Section 5 by means of a 3D Ultrasound flow profile measurement and the 2D analysis of Sunspot time series data.

In case of the analysis of flow profile measurements which are obtained in an experiment investigating the magnetorotational instability in liquid metals [22], such method would be highly desirable. Due to the complex experimental setup and the weak signal to noise ratio, the flow profile measurements contain several time intervals, in which distortions are dominant, see Seilmayer et al. [23]. These time intervals had to be rejected leading to a time series with invalid (missing) data points, from which the multivariate version of LSM is able to determine the parameters of the characteristic traveling wave.

In order to demonstrate the method and to motivate the basic idea of the multivariate version of LSM we compare it, in terms of necessary conditions and error (noise) behaviour, with the traditional orthogonal mode decomposition (OMD) with trigonometric basis functions, which is the essence of the classical Fourier transform. It will be shown that LSM fits better to the conditions of an arbitrary finite length of the sampling series in comparison to DFT, because LSM reduces the error in model parameter estimation. In contrast to the traditional approach, which was derived from a statistical point of view (see the appendix of Scargle [18]), the presented method is deduced from a technical point of view and focuses on its application. This leads to a slight change in the scaling of model parameters, which will be discussed in Section 4.2. However, the introduced procedure includes all benefits, like arbitrary sampling, fragmented data and a good noise rejection.

The starting point of this paper is the analysis of a continuous 1D signal $s:\mathbb{R}\to \mathbb{R}$ which is composed of an arbitrary and finite set of individual frequency components ${\omega }_i,0\le i\le M$. Without loss of generality, band limitation is assumed stating that there exists an upper maximum frequency ${\omega }_{max}$ with ${\omega }_i<{\omega }_{max},0\le i\le M$, which mimics an intrinsic low pass filter characteristic of the measurement device. Since the measurement time is finite, the value of s is only known in the time interval $\left[0,T\right]$ with $T\in \mathbb{R}$ and $T>0$. The m-dimensional extension of this signal is $S:{\mathbb{R}}^{\mathit{m}}\to \mathbb{R}$. If such a continuous signal is sampled, the pair $({\widehat{s}}_i,t_i)\in {\mathbb{R}}^2,0\le i\le N-1,$ represents the measured 1D value and the corresponding instant in time whereas the pair $({\widehat{S}}_i,{\overrightarrow{t}}_i)\in {\mathbb{R}}^{m+1},0\le i\le N-1,$ represents a measured value and the m-dimensional location (${\overrightarrow{t}}_i$ space/time) for the i-th sample (N is the number of samples). All the methods presented in this paper are implemented in a package written in R and published on CRAN [24]. The reader will find more elaborate examples in the package and in the supplementary material.

2. Mathematical Model and Comparison of OMD and LSM

In order to delineate the differences between trigonometric OMD and LSM, we start with the basic model of a periodic signal as a sum of signals of different frequencies ${\omega }_k,0\le k\le M,$ with the corresponding amplitude $A_k\in \mathbb{R}$ and phase shift ${\phi }_k\in \mathbb{R}$. The corresponding trigonometric model function

$$\begin{array}{cc}\hfill y(t)& =\sum _{k=0}^M{A}_{k}cos({\omega }_{k}t+{\phi }_k)=\sum _{k=0}^M(a_{k}cos({\omega }_{k}t)+b_{k}sin({\omega }_{k}t))\hfill \end{array}$$

(1)

describes an infinite, stationary and steady process $y:\mathbb{R}\to \mathbb{R}$ with the coefficients $a_k,b_k\in \mathbb{R}$ for the defined frequency ${\omega }_k$ and the identities $A_k=\sqrt{a_k^2+b_k^2}$ as well as ${\phi }_k={tan}^{-1}(b_k/a_k)$. Furthermore, the trigonometric model above consists of an arbitrary number $M\in {\mathbb{N}}_0$ of frequency components. If the given signal $s(t)$ is described by the defined model from Equation (1) the model misfit $ϵ$ is given by

$$ϵ(t)=s(t)-y(t).$$

(2)

Thus, the signal can now be described by inserting Equation (1) into Equation (2) and defining the misfit ${ϵ}_k$ for each discrete frequency k with $ϵ(t)={\sum }_k{ϵ}_k(t)$ as follows:

$$\begin{array}{cc}\hfill s(t)& =\sum _{k=0}^M\left(a_{k}cos({\omega }_{k}t)+b_{k}sin({\omega }_{k}t)+{ϵ}_k(t)\right).\hfill \end{array}$$

(3)

The defined misfit originates form measurement uncertainties or parametric errors from $a_k$, $b_k$. In the general case $ϵ(t)$ can be any function or distribution. The challenge is to precisely determine the model parameters $a_k$ and $b_k$ for a given signal $s(t)$ achieving minimal $ϵ(t)$. This can be accomplished by one of the three methods described in the following sections: (i) Least-Square fit; (ii) orthogonal mode decomposition; (iii) Lomb-Scargle method (LSM).

2.1. Least Square Fit

The optimal fit is reached by least square fitting, resulting in a minimum $ϵ$, which was shown by Barning [25], Mathias et al. [26]. Since such procedures are iterative, the convergence of the algorithm might need a large number of function evaluations of Equation (3). Therefore, a direct version like LSM is preferable. It can be shown that LSM becomes equivalent to a least square fit of a sinusoidal model [17,25].

2.2. Trigonometric OMD

Generally, two functions $f,g:\mathbb{R}\to \mathbb{R}$ are said to be orthogonal on the interval $\left[a,b\right]\subset \mathbb{R},$ if the following condition holds, refer to Weisstein [27]:

$${\int }_a^{b}f(x)g(x)\mathrm{d}x=0.$$

(4)

By selecting $f(x)=sin(x)$ and $g(x)=cos(x)$, the integration leads, by exploiting the identity $cos(x)sin(x)={\displaystyle \frac{1}{2}}sin(2x)$, to ${\displaystyle \frac{1}{4}}(cos(2a)-cos(2b))=0$, which is only zero, if $cos(2a)=cos(2b).$ This is true for any $a\in \mathbb{R}$, if the length of the interval $[a,b]$ is a multiple of the period $2\pi$ so that $2b=2a+2\pi k$ with $k\in {\mathbb{N}}^{+}$. It is interesting to note that this interval can be shorted to one half of the period, if $cos(2a)$ is zero. By exploiting this feature of the trigonometric functions, the individual model coefficients are calculated by multiplying the sine or the cosine to the measured data $s(t)$ and integrating over all times as shown by (Cohen [1], Chap. 15):

$$\begin{array}{cc}\hfill a_k& ={\displaystyle \frac{2}{T}}{\int }_{-\infty }^{\infty }s(t)cos({\omega }_{k}t)\mathrm{d}t\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill b_k& ={\displaystyle \frac{2}{T}}{\int }_{-\infty }^{\infty }s(t)sin({\omega }_{k}t)\mathrm{d}t.\hfill \end{array}$$

(6)

For a measured signal, the integration can only be performed over the interval $\left[0,T\right]$. It is obvious that an error is introduced, if $T\ne 2\pi n$. In order to investigate the properties of the finite integration, we will focus on the cosine term (Equation (5)), since the analysis of the sine term (Equation (6)) is similar. By setting the integral boundaries to the finite time interval and expressing the signal using Equation (3), the integral in Equation (5) can be written as

$$\begin{array}{cc}\hfill & {\int }_0^{T}s(t)cos({\omega }_{k}t)dt=\hfill \\ \hfill & {\int }_0^T(a_k{cos}^2({\omega }_{k}t)+b_{k}sin({\omega }_{k}t)cos({\omega }_{k}t)+{ϵ}_k(t)cos({\omega }_{k}t))\mathrm{d}t.\hfill \end{array}$$

(7)

Using the identities ${cos}^2(x)={\displaystyle \frac{1}{2}}\left(1+cos(2x)\right)$ and $sin(x)cos(x)={\displaystyle \frac{1}{2}}sin(2x)$, Equation (7) becomes

$$\begin{array}{cc}\hfill {\displaystyle \frac{2}{T}}{\int }_0^{T}s(t)cos({\omega }_{k}t)\mathrm{d}t& =a_k(1+\underset{\mathrm{truncation}\mathrm{error}{ϵ}_k^{\mathrm{T}}}{\underbrace{{\displaystyle \frac{1}{T}}{\int }_0^T(cos(2{\omega }_{k}t)+{\displaystyle \frac{b_k}{a_k}}sin(2{\omega }_{k}t))\mathrm{d}t})}\hfill \\ \hfill & +\underset{\mathrm{random}\mathrm{error}{ϵ}_k^{\mathrm{FS}}}{\underbrace{{\displaystyle \frac{2}{T}}{\int }_0^T{ϵ}_k(t)cos({\omega }_{k}t)\mathrm{d}t}}.\hfill \end{array}$$

(8)

The equation above indicates that the coefficient $a_k$ is effected by two errors ${ϵ}_k^{\mathrm{T}}$ (truncation) and ${ϵ}_k^{\mathrm{FS}}$ (random) which might be nonzero for an arbitrary T. Thus, if both errors are neglected, the integral on the left hand side turns into an approximation of $a_k$. Taking the consideration about the orthogonality described in Equation (4) into account, ${ϵ}_k^{\mathrm{T}}$ becomes only zero, if the integration time T is an integer multiple of $\pi /{\omega }_k$. This means for the Fourier series that a given integration time T (observation window) defines the lowest allowed frequency ${\omega }_0$. Furthermore, only multiples of this fundamental frequency ${\omega }_k=k{\omega }_0$ with $k\in \mathbb{N}$ are allowed in the Fourier series because ${ϵ}_k^{\mathrm{T}}=0$ in this case.

The technical realization of OMD is called “quadrature demodulation”, which is the discrete version of equation (8) by replacing the integral over $s(t)$ into a sum over the discrete sampled values ${\widehat{s}}_n,0\le n\le N-1$ and setting the truncation error to zero. For a time discrete signal with N samples and equidistant sampling with a constant sampling period $T_s$, and $t_n=n{T}_s$, the parameter $a_k$ (corresponding to the frequency ${\omega }_k$) is approximated by

$$a_k\approx {\displaystyle \frac{2}{N}}\sum _{n=0}^{N-1}{\widehat{s}}_{n}cos({\omega }_k{t}_n).$$

(9)

According to the previous considerations, the truncation error ${ϵ}_k^{\mathrm{T}}$ is larger than zero if the measurement range $T=N{T}_{\mathrm{s}}$ is not exactly a multiple of $\pi /{\omega }_k$, as shown in Figure 1.

In this example a signal with a time period of $t=2\pi$ is sampled with $N=23$, $T_s=2\pi /20$. The total integration time is $T=2\pi +1/5\pi$, being slightly longer than the period of the signal, which is indicated by the two additional sampling points after $t=2\pi$. The truncation error equals the gray shaded area.

Generally speaking, the sampling error in the discrete version can be estimated by

$${\displaystyle \frac{2}{N}}\sum _{n=0}^{N-1}{\widehat{s}}_{n}cos({\omega }_k{t}_n)\approx a_k(1+\underset{{ϵ}_k^{\mathrm{T}}}{\underbrace{{\displaystyle \frac{\Delta \phi }{T}}}})\pm \underset{{ϵ}_k^{\mathrm{FS}}}{\underbrace{{\Phi }_{1-\alpha }{\displaystyle \frac{2\sigma }{\sqrt{N}}}}},$$

(10)

where $\Delta \phi ={min}_i(T-\pi i/{\omega }_k)$, $i\in \mathbb{N}$. The random error is modeled by the $\alpha$-quantile of the sampling error distribution ${\Phi }_{1-\alpha }$ related to the underlying process with its standard deviation $\sigma$ (e. g. JCGM [28]). This estimates the confidence interval of $a_k$ and $b_k$, respectively (see supplementary material for details).

The above considerations lead to following four statements, which are derived in detail in Section 2 of the supplementary material : (i) the maximum absolute value of the truncation error is bounded $\left|{ϵ}_k^{\mathrm{T}}\right|\le 0.2$ if $T\pi /{\omega }_k$, which is consistent with the experimental findings from Thompson and Tree [29]; (ii) ${ϵ}_k^{\mathrm{T}}$ is independent from the total number of samples N for a constant time T, which makes the OMD a non-consistent estimator for model parameters $a_k$ and $b_k$; (iii) by increasing the sample rate, the random error decreases by $\mathcal{O}(N^{-0.5})$ but scales with twice the standard deviation of noise distribution; (iv) DFT can be derived from Equation (10) by restricting to equidistant sampling with constant $T_{\mathrm{s}}$.

Concluding, the truncation error, which is an intrinsic feature of OMD, produces a systematic deviation from the true value depending on the difference between the sampling time interval and the corresponding period of the frequency of interest. In contrast to that, the random error diminishes with increasing sampling rate $T_{\mathrm{s}}^{-1}=N/T$. A detailed discussion on this topic can be found in Thompson and Tree [29], Jerri [30], Shannon [31]. In case of multivariate and randomly sampled data the recent work of Al-Ani and Tarczynski [32] suggests two calculation schemes estimating the Fourier transform in the general case. The continuous time Fourier transform estimation (similar to Equation (9)) takes samples with arbitrary spacing as the general approach, in contrast to the secondly proposed discrete time Fourier transform ($t_n=n{T}_{\mathrm{s}}$ and ${\omega }_k=k{\omega }_0$) estimation scheme. Here the data is projected onto a regular grid, which now may provide locations of missing data. Both schemes consider the sampling pattern $t_n$ and its power spectrum distribution function $p(t_n)$. However, even for such sophisticated methods, the main conceptual drawbacks (see points (i) and (ii)) remain, which motivates the subsequent Lomb-Scargle method as non OMD method and its extension to multivariate data.

2.3. Lomb Scargle Method

As described in the previous section, the disadvantage of OMD is that cosine and sine are not orthogonal on arbitrary intervals. By introducing an additional parameter $\tau \in \mathbb{R}$ into Equation (4) it will be shown that

$${\int }_a^{b}cos(x-\tau )sin(x-\tau )\mathrm{d}x=0$$

(11)

holds for arbitrary intervals $\left[a,b\right]$. By utilizing the trigonometric identities, ${cos}^2(\varphi )-{sin}^2(\varphi )=cos(2\varphi )$ and $cos(\varphi )sin(\varphi )={\displaystyle \frac{1}{2}}sin(2\varphi )$, to remove the differences in the arguments the integral can be transformed to

$${\int }_a^b\left({\displaystyle \frac{1}{2}}sin\left(2x\right)cos\left(2\tau \right)-{\displaystyle \frac{1}{2}}sin\left(2\tau \right)cos\left(2x\right)\right)\mathrm{d}x=0.$$

If this integral has to be zero, the following condition has to hold

$$cos(2\tau ){\int }_a^{b}sin(2x)\mathrm{d}x=sin(2\tau ){\int }_a^{b}cos(2x)\mathrm{d}x$$

which can be transformed to

$${\displaystyle \frac{{\int }_a^{b}sin\left(2x\right)\mathrm{d}x}{{\int }_a^{b}cos\left(2x\right)\mathrm{d}x}}=tan\left(2\tau \right).$$

(12)

Therefore, from Equation (12) the parameter $\tau$ verifying Equation (11) may be found. In case of equidistant sampling, the value of $\tau$ can be directly calculated from the integration boundaries by:

$$\tau ={\displaystyle \frac{b-a}{2}}.$$

Based on this general consideration and in order to remove the truncation error, for each frequency ${\omega }_k$ the time shifting parameter ${\tau }_k\in \mathbb{R},0\le k\le M,$ was introduced by Lomb and Scargle [17,18] into the model given in Equation (1) so

$$s(t)=\sum _{k=0}^M\left(a_{k}cos\left({\omega }_k(t-{\tau }_k)\right)+b_{k}sin({\omega }_k(t-{\tau }_k))+{ϵ}_k(t)\right).$$

(13)

The parameter ${\tau }_k$ can be calculated by

$$tan(2{\omega }_k{\tau }_k)={\displaystyle \frac{{\int }_a^{b}sin(2{\omega }_{k}t)\mathrm{d}t}{{\int }_a^{b}cos(2{\omega }_{k}t)\mathrm{d}t}}$$

(14)

similar to Equation (12). For the time discrete version the integrals transform into a sum resulting in

$$tan(2{\omega }_k{\tau }_k)={\displaystyle \frac{{\sum }_{n=0}^{N-1}sin(2{\omega }_k{t}_n)}{{\sum }_{n=0}^{N-1}cos(2{\omega }_k{t}_n)}}.$$

(15)

The parameters $a_k$ can be determined beginning with Equation (7) but by factorising them by ${cos}^2(x)$, instead of expanding ${cos}^2(x)$ to ${\displaystyle \frac{1}{2}}\left(1+cos(2x)\right)$ as done for Equation (8). In the following, the method is delineated for the discrete set ${\widehat{s}}_i$ since it is applied to sampled data. This data is multiplied by $cos({\omega }_k(t_n-{\tau }_k))$ resulting in the next equation for a single frequency ${\omega }_k$:

$$\sum _{n=0}^{N-1}{\widehat{s}}_{n}cos({\varphi }_k)=a_k\sum _{n=0}^{N-1}{cos}^2({\varphi }_k)+b_k\underset{\mathrm{truncation}\mathrm{error}{ϵ}_k^{\mathrm{T}}}{\underbrace{\sum _{n=0}^{N-1}sin({\varphi }_k)cos({\varphi }_k)}}+\underset{\mathrm{truncation}\mathrm{error}{ϵ}_k^{\mathrm{T}}}{\underbrace{\sum _{n=0}^{N-1}{ϵ}_k(t_n)cos({\varphi }_k)}},$$

(16)

where ${\varphi }_k={\omega }_k(t_n-{\tau }_k)$. The sum of the terms $cos({\varphi }_k)sin({\varphi }_k)$ vanishes because of the proper selection of ${\tau }_k$ according to Equation (15). The sum of the terms ${ϵ}_k(t_n)cos({\varphi }_k)$ describes a modulated noise distribution function. The value of the parameters $a_k$ can be directly obtained by dividing by ${\sum }_{n=0}^{N-1}{cos}^2({\varphi }_k)$. This gives

$${\displaystyle \frac{{\sum }_{n=0}^{N-1}{\widehat{s}}_{n}cos({\varphi }_k)}{{\sum }_{n=0}^{N-1}{cos}^2({\varphi }_k)}}=a_k+\underset{{ϵ}_k^{\mathrm{LS}}}{\underbrace{{\displaystyle \frac{{\sum }_{n=0}^{N-1}{ϵ}_k(t_n)cos({\varphi }_k)}{{\sum }_{n=0}^{N-1}{cos}^2({\varphi }_k)}}}}.$$

(17)

The residual error term on the right hand side, ${ϵ}_k^{\mathrm{LS}}$, consists of a random distributed part divided by a sum over the square of the cosine. In contrast to OMD, the estimation error of the parameter $a_k$ only depends on noise and is independent from the realization of sampling.

The next step is to determine ${ϵ}_k^{\mathrm{LS}}$ in terms of a confidence interval with respect to $\alpha$-quantile of the sampling error distribution ${\Phi }_{1-\alpha }$, as it was done for OMD in Equation (10). Given a normal distributed error function with zero expectation value $\mu =0$ and a standard deviation $\sigma$ independent from time (${ϵ}_k\leftrightarrow \mathcal{N}(0,\sigma )$), the error can be factorized (Parzen [33], Theorem 4A, p. 90), leading to

$$\begin{array}{cc}\hfill {ϵ}_k^{\mathrm{LS}}& ={\Phi }_{1-\alpha }{\displaystyle \frac{\sigma }{\sqrt{N}}}{\displaystyle \frac{{\sum }_{n=0}^{N-1}cos({\varphi }_k)}{{\sum }_{n=0}^{N-1}{cos}^2({\varphi }_k)}}\le {\Phi }_{1-\alpha }{\displaystyle \frac{\sigma }{\sqrt{N}}}{e}_{max}\hfill \end{array}$$

estimating the most probable limits of ${ϵ}_k^{\mathrm{LS}}$. In the above equation $e_{max}$ neither depend on the frequency nor the shifting parameter and is defined as

$$e_{max}=\underset{\beta \in [0,2\pi ]}{max}\left[{\displaystyle \frac{{\int }_0^{\beta }cos(\varphi )\mathrm{d}\varphi }{{\int }_0^{\beta }{cos}^2(\varphi )\mathrm{d}\varphi }}\right]=\underset{\beta \in [0,2\pi ]}{max}\left[{\displaystyle \frac{4sin(\beta )}{(2\beta +sin(2\beta ))}}\right]={\displaystyle \frac{4}{\pi }}.$$

The final parameter estimation gives

$$a_k={\displaystyle \frac{{\sum }_{n=0}^{N-1}{\widehat{s}}_{n}cos({\omega }_k(t_n-{\tau }_k))}{{\sum }_{n=0}^{N-1}{cos}^2({\omega }_k(t_n-{\tau }_k))}}\pm \underset{\Delta a_k}{\underbrace{{\displaystyle \frac{4}{\pi }}{\Phi }_{1-\alpha }{\displaystyle \frac{\sigma }{\sqrt{N}}}}},$$

(18)

where the confidence interval $\Delta a_k={\displaystyle \frac{4}{\pi }}{\Phi }_{1-\alpha }{\displaystyle \frac{\sigma }{\sqrt{N}}}$ converges to zero by increasing the number of samples N, for a fixed time interval. This property qualifies LSM to be a consistent estimator in amplitude and phase [26]. In addition, the confidence interval $\Delta a_k$ for the LSM is smaller than the value given in Equation (10) for the OMD (this also occurs for the confidence interval of the parameters $b_k$, $\Delta b_k=\Delta a_k$). We notice that the definition of $a_k$ given in Equation (18) is slightly different from the original definition given by Lomb [17] but both definitions converge for large N (see supplementary material).

The confidence interval for the amplitude, $A_k=\sqrt{a_k^2+b_k^2}$, can be deduced by propagating the error $\Delta a_k$:

$$\begin{array}{cc}\hfill \Delta A_k& =\left|{\displaystyle \frac{\partial A_k}{\partial a_k}}\right|\Delta a_k+\left|{\displaystyle \frac{\partial A_k}{\partial b_k}}\right|\Delta b_k={\displaystyle \frac{4}{\pi }}{\Phi }_{1-\alpha }\sqrt{{\displaystyle \frac{2}{N}}}\sigma .\hfill \end{array}$$

In a similar way, the confidence interval for the phase ${\phi }_k$ can be defined by

$$\begin{array}{cc}\hfill {\phi }_k& ={tan}^{-1}\left({\displaystyle \frac{b_k}{a_k}}\right)\pm \Delta \phi ={tan}^{-1}\left({\displaystyle \frac{b_k}{a_k}}\right)\pm {\displaystyle \frac{4}{\pi }}{\Phi }_{1-\alpha }\sqrt{{\displaystyle \frac{2}{N}}}{\displaystyle \frac{\sigma }{A_k}}.\hfill \end{array}$$

It is interesting to note that the confidence interval for the phase ${\phi }_k$ decreases for increasing amplitude.

3. Spectrum from the Lomb-Scargle Method

In order to determine the spectrum with LSM, we assume the recorded signal is described by many frequencies. The most significant frequency is then represented by a peak in the frequency spectrum of a certain width and height. The width is determined by the frequency resolution $\Delta f$ which equals $1/T$. From this point of view the precision of a frequency estimation changes only with observation length T and seams to be independent from the number of samples N and signal quality.

The quality $\Sigma$ is measured by a signal-to-noise ratio like expression

$$\begin{array}{cc}\hfill \Sigma & =\sqrt{{\displaystyle \frac{1}{N}}\sum _{n=0}^{N-1}{\displaystyle \frac{{\left(s_n-y_n\right)}^2}{{\sigma }_n^2}}}\approx {\displaystyle \frac{\sqrt{{\sum }_l{A}_l^2}}{\sqrt{\frac{1}{N}{\sum }_{n=0}^{N-1}{\left(ϵ(t_n)\right)}^2}}}\hfill \end{array}$$

(19)

with $A_l$ counting the significant amplitudes. In Equation (19) $y_n$ is the fitted model, ${\sigma }_n$ is the uncertainty per sample, with $\sigma =\sqrt{{\displaystyle \frac{1}{N}}{\sum }_{n=0}^{N-1}{\sigma }_n^2}$, and $ϵ(t_n)$ is the noise per sample. Following VanderPlas [34], we apply Bayesian statistics and assume that every peak is Gaussian shaped, i. e. ${\mathrm{e}}^{P(f_{max}\pm \Delta f)}\propto {\mathrm{e}}^{-\Delta f^2/\left(2{\sigma }_f^2\right)}$. It follows that a significant peak $A_{max}^2=A^2(f_{max})$ appears at $f_{max}$, in a way that $A_{max}^2/2=A^2(f_{max}\pm \Delta f)$ is valid. Here, $A^2$ is related to the power spectral density $P(f)\propto a_k^2+b_k^2$. The frequency uncertainty (or standard deviation) is then given by

$${\sigma }_f\approx \Delta f\sqrt{{\displaystyle \frac{2}{N{\Sigma }^2}}}$$

so that a significant peak is located in the interval $f_{max}\pm {\sigma }_f$. This approach suggests that increasing the number of samples in a fixed interval T enhances the precision of $f_{max}$ by reducing ${\sigma }_f.$ However, if the original signal contains two frequencies with a distance in the range of $1/T$, it cannot be excluded even for LSM that these peaks merge together in one single peak with small ${\sigma }_f$ according to Kovács [35].

Concluding the properties of LSM: (i) there is no truncation error ${ϵ}_{\mathrm{T}}=0$; (ii) LSM provides a better noise rejection compared to OMD, ${ϵ}^{\mathrm{LS}}<{ϵ}^{\mathrm{FS}}$; (iii) the explicit sampling pattern is not of interest to work with LSM.

4. The Multivariate Lomb-Scargle Method

A signal S depending of m-independent variables represents a function $S(\overrightarrow{t}):{\mathbb{R}}^m\to \mathbb{R}$ with the input vector described by $\overrightarrow{t}=[t_1,t_2,\dots ,t_m]$. The model function for multivariate LSM is gained by replacing the scalar arguments of cosine and sine of the univariate model function of Equation (13) by vectors, resulting in

$$Y(\overrightarrow{t})=\sum _{k=0}^M\left(a_{k}cos\left({\overrightarrow{\omega }}_k·\left(\overrightarrow{t}-{\overrightarrow{\tau }}_k\right)\right)+b_{k}sin\left({\overrightarrow{\omega }}_k·\left(\overrightarrow{t}-{\overrightarrow{\tau }}_k\right)\right)\right).$$

In this case, the shifting parameter ${\overrightarrow{\tau }}_k\in$${\mathbb{R}}^{\mathit{m}},0\le k\le M,$ is a vector and in principle, hard to calculate. However, if the argument of the cosine is expanded, it is obvious that the scalar product ${\overrightarrow{\omega }}_k·{\overrightarrow{\tau }}_k\in \mathbb{R}$ does not depend on time and thus, the cosine argument can be written as ${\overrightarrow{\omega }}_k·\overrightarrow{t}-{\tau }_k^{*}$ with ${\tau }_k^{*}={\overrightarrow{\omega }}_k·{\overrightarrow{\tau }}_k$. The determination of ${\tau }_k^{*}$ is similar as for $\tau$ shown in Equation (14), but with some differences in the equations since the phase, instead of the time coordinate, is shifted now. This determination is described in the following.

4.1. Derivation of the Shifting Parameter

This section shows that shifting the phase, instead of time, does not affect the Lomb-Scargle algorithm. The derivation of the shifting parameter for the multivariate case is delineated for the time discrete signal $\widehat{S}$$=\{({\widehat{S}}_i,{\overrightarrow{t}}_i)\in R^{\mathit{m}+\mathit{1}},0\le i\le N-1\}.$ Starting with the orthogonality condition

$$\sum _{n=0}^{N-1}sin\left({\overrightarrow{\omega }}_k·{\overrightarrow{t}}_n-{\tau }_k^{*}\right)cos\left({\overrightarrow{\omega }}_k·{\overrightarrow{t}}_n-{\tau }_k^{*}\right)=0$$

and applying the trigonometric identities to remove the differences in the arguments:

$$\begin{array}{c}\hfill \sum _{n=0}^{N-1}\left[\left(\underset{{\mathrm{c}}_n}{\underbrace{cos\left({\overrightarrow{\omega }}_k·{\overrightarrow{t}}_n\right)}}\underset{{\mathrm{c}}_{Ø}}{\underbrace{cos\left({\tau }_k^{*}\right)}}+\underset{{\mathrm{s}}_n}{\underbrace{sin\left({\overrightarrow{\omega }}_k·{\overrightarrow{t}}_n\right)}}\underset{{\mathrm{s}}_{Ø}}{\underbrace{sin\left({\tau }_k^{*}\right)}}\right)\right.\\ \hfill \left.\left(\underset{{\mathrm{s}}_n}{\underbrace{sin\left({\overrightarrow{\omega }}_k·{\overrightarrow{t}}_n\right)}}\underset{{\mathrm{c}}_{Ø}}{\underbrace{cos\left({\tau }_k^{*}\right)}}-\underset{{\mathrm{c}}_n}{\underbrace{cos\left({\overrightarrow{\omega }}_k·{\overrightarrow{t}}_n\right)}}\underset{{\mathrm{s}}_{Ø}}{\underbrace{sin\left({\tau }_k^{*}\right)}}\right)\right]& =0\hfill \end{array}$$

By using the defined abbreviations and rearranging the summation (the terms ${\mathrm{c}}_{Ø}$ and ${\mathrm{s}}_{Ø}$ do not depend on n), the following equation is easily deduced:

$${\displaystyle \frac{{\mathrm{c}}_{Ø}{\mathrm{s}}_{Ø}}{{\mathrm{c}}_{Ø}^2-{\mathrm{s}}_{Ø}^2}}={\displaystyle \frac{{\sum }_{n=0}^{N-1}{\mathrm{c}}_n{\mathrm{s}}_n}{{\sum }_{n=0}^{N-1}({\mathrm{c}}_n^2-{\mathrm{s}}_n^2)}}$$

The fraction on both sides can be further simplified by applying ${cos}^2(x)-{sin}^2(x)=cos(2x)$, $cos(x)sin(x)={\displaystyle \frac{1}{2}}sin(2x)$, and the definition of the tangent. This results in

$$\begin{array}{cc}\hfill tan(2{\tau }_k^{*})& ={\displaystyle \frac{{\sum }_{n=0}^{N-1}sin(2{\overrightarrow{\omega }}_k·{\overrightarrow{t}}_n)}{{\sum }_{n=0}^{N-1}cos(2{\overrightarrow{\omega }}_k·{\overrightarrow{t}}_n)}}.\hfill \end{array}$$

(20)

Notice that in comparison with one dimensional LSM, the frequency ${\omega }_k$ is missing on the left side in Equation (20).

4.2. Parameter Estimation

Similar to the procedure for the LSM for a discrete signal in one dimension (see Equation (16)), a discrete multivariate signal $({\widehat{S}}_i,{\overrightarrow{t}}_i)$ with N samples is multiplied by $cos({\overrightarrow{\omega }}_k·\overrightarrow{t_i}-{\tau }_k^{*})$ resulting in

$$\begin{array}{c}\sum _{n=0}^{N-1}\widehat{S}({\overrightarrow{t}}_n)cos({\overrightarrow{\omega }}_k·{\overrightarrow{t}}_n-{\tau }_k^{*})=a_k\sum _{n=0}^{N-1}{cos}^2({\overrightarrow{\omega }}_k·{\overrightarrow{t}}_n-{\tau }_k^{*})\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+b_k\sum _{n=0}^{N-1}\underset{=0,\mathrm{if}\mathrm{orthogonal}}{\underbrace{sin({\overrightarrow{\omega }}_k·{\overrightarrow{t}}_n-{\tau }_k^{*})cos({\overrightarrow{\omega }}_k·{\overrightarrow{t}}_n-{\tau }_k^{*})}}+\sum _{n=0}^{N-1}\underset{\mathrm{random}\mathrm{error}}{\underbrace{{ϵ}_k({\overrightarrow{t}}_n)cos({\overrightarrow{\omega }}_k·{\overrightarrow{t}}_n-{\tau }_k^{*})}}.\end{array}$$

The determination of the parameter $a_k$, and its confidence interval $\Delta a_k$, are similar to the one-dimensional case shown in Equation (18):

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\sum }_{n=0}^{N-1}y({\overrightarrow{t}}_n)cos({\overrightarrow{\omega }}_k·{\overrightarrow{t}}_n-{\tau }_k^{*})}{{\sum }_{n=0}^{N-1}{cos}^2({\overrightarrow{\omega }}_k·{\overrightarrow{t}}_n-{\tau }_k^{*})}}& =a_k,\Delta a_k={\displaystyle \frac{4}{\pi }}{\Phi }_{1-\alpha }{\displaystyle \frac{\sigma }{\sqrt{N}}}.\hfill \end{array}$$

The parameter $b_k$, and its confidence interval $\Delta b_k$, are calculated analogously:

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\sum }_{n=0}^{N-1}y({\overrightarrow{t}}_n)sin({\overrightarrow{\omega }}_k·{\overrightarrow{t}}_n-{\tau }_k^{*})}{{\sum }_{n=0}^{N-1}{sin}^2({\overrightarrow{\omega }}_k·{\overrightarrow{t}}_n-{\tau }_k^{*})}}& =b_k,\Delta b_k=\Delta a_k.\hfill \end{array}$$

We recall that the $\alpha$-quantile of the sampling error distribution ${\Phi }_{1-\alpha }$ (with standard deviation $\sigma$) has been used to model the random error. The power spectral density, as well as the false alarm probability, are calculated in the same way compared to the one dimensional LSM method (see supplementary material).

An implementation of the multivariate Lomb-Scargle method is available, in the spectral (R) package published on CRAN [24], to give the user access to a multi-dimensional analysis with an easy to use interface.

5. Application

The selected applications are ordered by increasing complexity of the sampling, and the data itself. First, a regular grid with missing values for a synthetic test data is considered. Second, Ultrasound Doppler Velocimetry (UDV) measurement data, which contains jitter and missing values is analized. Finally, a quite general data scenario given by an astrophysical 2D set of Sunspots, which appear freely in time and space, is investigated.

5.1. Synthetic Test Data

We start with a simple test case, where the input signal is a simple two dimensional plain wave $z=cos\left(2\pi (x{f}_x+y{f}_y)+{\displaystyle \frac{\pi }{4}}\right)$, with $f_x=3.25$ and $f_y=6.32$ representing the dimensionless frequencies, specifically selected in such a way that $f_x,f_y\ne k·{\omega }_0$, with ${\omega }_0=2\pi /T_{x,y}$ and $T_{x,y}=1$. The 2D sampling of $x,y\in [0,1]$ is performed with $\delta x,\delta y=0.025$. With respect to the chosen frequencies $f_x$ and $f_y$ it becomes evident, that both frequencies do not fit to the data range by an integer fraction. Data gaps are introduced by removing randomly distributed and uncorrelated grid points covering 60% of the total data set. Figure 2(a) shows the nonuniform distributed input data. Gray areas indicate missing values (non available numbers “$\mathrm{nan}$”). As LSM requires an appropriate input vector of frequencies we choose for both variables $f_x,f_y\in [-10,10]$ with a resolution of $\delta f_x,\delta f_y=0.025$.

The corresponding power spectral density (PSD) from the LSM (see supplementary material), shown in Figure 2(b), displays the maxima at the 1st and 3rd quadrant, which represents a wave traveling upwards. The figure illustrates that even with a huge amount of missing data, it is possible to properly detect the periodic signal. In the present case, the value of $\mathrm{psd}(f_x,f_y)\approx 1$ indicates a perfect fit to the corresponding sinusoidal model. For comparison, Figure 2(c) depicts the standardized PSD

$${\mathrm{psd}}_{\mathrm{FT}}(\overrightarrow{\omega })={\displaystyle \frac{N}{N-1}}·{\displaystyle \frac{1}{2{\sigma }_0^2}}{\left({\displaystyle \frac{2\left|\mathcal{F}\left(z(x,y)\right)(\overrightarrow{\omega })\right|}{1-\frac{N_{\mathrm{zero}}}{N}}}\right)}^2$$

calculated from a discrete Fourier transform DFT symbolized by the operator $\mathcal{F}(z(x,y))(\overrightarrow{\omega })=\int z(x,y){\mathrm{e}}^{\mathrm{i}({\omega }_{x}x+{\omega }_{y}y)}\mathrm{d}x\mathrm{d}y$. Here, missing values ($\mathrm{nan}$) are handled by zero padding which means that the introduced gaps are filled with zero. Since the DFT represents a conservative transformation (mapping) into the spectral domain, zero padding leads to a significant reduction of the average signal energy. A proper rescaling to the number of zeros $N_{\mathrm{zero}}$ (missing values) is required, to ensure approximated amplitude estimations. Zero padding always changes the character of the input signal, in a way that the corresponding DFT threats the zeros as if they where part of the original signal. The resulting PSD becomes then finally different from the original one. Since the signal frequencies $f_x,f_y$ do not fit into an integer spaced scheme of the data range, the DFT suffers from its drawbacks, as is briefly discussed in Section 2.2. The effect of leakage and misfit of frequencies can be recognized in Figure 2(c), in terms of a rather coarse resolution and a lower PSD value.

5.2. 3D UDV Flow Measurement

The second example is taken from previous experimental flow measurements [22] devoted to investigate the azimuthal magneto-rotational instability AMRI. The experiment consists of a cylindrical annulus containing a liquid metal between the inner and the outer wall. The inner wall rotates at a frequency ${\omega }_{\mathrm{in}}=2\pi ·0.05\mathrm{H}\mathrm{z}$ whereas the outer wall rotates at a frequency ${\omega }_{\mathrm{out}}=2\pi ·0.013\mathrm{H}\mathrm{z}$. The flow is driven by the shear, since ${\omega }_{\mathrm{in}}-{\omega }_{\mathrm{out}}>0$. Two ultrasound sensors mounted at the outer cylinder, on opposite side of each other and at the same radius, measure the axial velocity component $v_z$ along the line of sight parallel to the rotation axis of the cylinder. Additionally, the liquid metal flow is exposed to a magnetic field $B_{\phi }\propto r^{-1}$ (r is the radial coordinate), originated from a current $I_{\mathrm{axis}}$ on the axis of the cylinder. A typical time series of the axial velocity (along the measuring line) is displayed in Figure 3. The existence of periodic patterns (large blue inclined stripes) in this figure indicates a travelling wave propagating in the fluid. In cylindrical geometry, a travelling wave is described in terms of the drift frequency $\omega =2\pi f$, the vertical wave numbers k, and the azimuthal wave numbers m. A close look to Figure 3 evidences that this wave is not axially symmetric, and that the leading azimuthal wave number is $m=1$. This is indeed proven in the following with the LSM analysis.

The measured time series for $v_z$ depends on time $t_n$, depth $d_n$ and the angular positions of the sensors ${\phi }_n^{{\mathrm{S}}_1}={\omega }_{\mathrm{out}}{t}_n$ and ${\phi }_n^{{\mathrm{S}}_2}={\phi }_n^{{\mathrm{S}}_1}+\pi$, where the subscript n indicates the sampling. The azimuthal angles depend on time, since the sensors are attached to the outer wall. This induces the mapping $v_z=f(t_n,d_n,{\phi }_n)$, where ${\phi }_n=\{{\phi }_n^{{\mathrm{S}}_1},{\phi }_n^{{\mathrm{S}}_2}\}$ is meant to be one out of the set of two sensors at each sampling instance. Finally the measured velocity component is a mapping $v_z:{\mathbb{R}}^3\mapsto \mathbb{R}$, which fits perfectly to the proposed multivariate LSM.

Figure 4 illustrates the result of the LSM decomposition into several azimuthal m-modes with the corresponding amplitudes. The $m=0$ mode contains a stationary structure, at a frequency $f\approx 0$, which originates from sensor miss-alignments and (thermal) side effects in the flow. The minor non-stationary components, here the two point-symmetric peaks, originate from cross-talk (or alias projections, see labels on Figure 4) of the $m=1$ mode. This is mainly for two reasons: First, because the sensors behave not exactly identical (e. g. misalignment or different sensitivity) and therefore respond slightly different to the same flow signal. This means that one of the sensors projects a little bit more energy into the data than the other, which leads to a “leakage”-effect and a weak signal in $m=0$ mode. Second, the rotation of the sensors acts like an additional sampling frequency $f_{\mathrm{out}}$. In consequence, the spectrum of this regular sampling function folds with the spectrum of the observed process, leading to alias images (copies) of the original process spectrum into other frequency ranges (i. e. m’s). In the present multivariate analysis, these alias structures exhibit point symmetry in $m=0$ panel, indicating corresponding traveling waves. It is evident that these aliases originate as a systematic error from the LSM method and, due to their symmetry, effectively cancel each other out. One could also argue, that multivariate LSM enables a robust identification of alias amplitudes, which do not belong to the observed process.

The AMRI wave itself is located in the $m=1$ panel, with a characteristic frequency in time (f) and in the vertical direction (k). Here, two components can be identified: the dominant wave at $f\approx 9m\mathrm{H}\mathrm{z}$ and $k\approx 20m^{-1}$ and a minor counterpart at $f\approx 4m\mathrm{H}\mathrm{z}$ and $k\approx -20m^{-1}$. The signals observed in the $m=3$ and $m=5$ panels might be due to aliases, because of the rotation of the inner cylinder at a frequency $f_{\mathrm{i}}=0.05\mathrm{H}\mathrm{z}$.

The advantage of “high” dimensional spectral decomposition is the improved noise rejection. With respect to the raw data, given in Figure 3, it is obvious that high frequency noise is present in the data. Depending on the exact distribution of this noise, its energies spread over a certain range of frequencies. If we would select a representative depth and angle so that $d_n,{\phi }_n=\mathrm{const}.$, the velocity $v_z$ only depends on time. In the subsequent one dimensional analysis the noise would accumulate along the single frequency ordinate probably hiding the signal of interest. Taking the higher order analysis distributes the noise energies over multiple domain variables. For the present example, this means that the distortions are projected into higher frequencies f, higher m’s and larger k’s. Since the signal of interest remains in the same spectral corridor, its signal amplitude becomes more clear, because of the “reduced” local noise.

5.3. Analyzing 2D s Unspot Data

Since the beginning of the 17th century systematic visual observations of Sunspots are available. This allows to investigate dynamic processes taking place in the Sun. The appearance of Sunspots on the Sun’s surface depends on the level of solar activity, and therefore it renders some fundamental features of the underlying solar dynamo, such as the $11\mathrm{y}\mathrm{r}\mathrm{s}$ solar cycle. Since the beginning of the 1820’s observational data is available, in terms of a two dimensional time series displayed in Figure 5. The data exhibits a periodic wing-like pattern essentially symmetric with respect to the Sun’s equator, the Sunspot butterfly diagram. This diagram summarizes the individual Sunspot groups appearing, at a certain time and latitude on the Sun, for the past 190 years. Additionally, Figure 5 depicts the assigned field polarity order indicated by the Sunspot color (gray or black). Given a year Y and a mean latitude L (degrees), the field polarity $P(Y,L)\in \{-1,1\}$ gives the arrangement of leading North/South polarity of a Sunspot group. This polarity is changing approximately each 11 yrs, which is called Schwabe’s cycle. The full period of $22\mathrm{y}\mathrm{r}\mathrm{s}$ is called the Hale cycle. The segmentation of Figure 5 is obtained following suggestions from Leussu et al. [36], but with some simplifications leading only to minor miss-assignments for individual Sunspot groups.

The Sunspot data is taken as-is from Leussu et al. [37], which originates from several sources with different qualities, i. e. the Royal Greenwich Observatory – USAF/NOAA(SOON)1[38,39], Schwabe2[40] and Spoerer3[43] data sets. A detailed discussion of the data collection is given by Leussu et al. [36], Leussu et al. [37] and the referenced literature.

In this section, the LSM spectral decomposition from this unevenly sampled binary data set is obtained, in order to identify typical periods. In contrast to the two examples analysed in Section 5.1 and Section 5.2, the Sunspot data set consists of real arbitrary sampling in both, time and space, providing the most general scenario input data for the LSM. The starting point for the subsequent analysis is the simplified segmentation of the butterfly diagram reassigning the field polarity. The segmentation takes place by optimizing the distance d (in the $(L,Y)$-plane) between each Sunspot and the segmentation line

$$L(Y)=\left\{\begin{array}{cc}{m}_1\left(Y-T_0\right)\hfill & m_1<0,L>0\hfill \\ m_2\left(Y-T_0\right)\hfill & m_2>0,L\le 0\hfill \end{array}\right.,$$

which is a piece wise linear function with intersection point $(T_0,0)$. The individual slopes, with $|m_1|,|m_2|<12deg/yrs$, are assigned on the northern and southern hemispheres, respectively. The final segmentation of Figure 5 consist of 18 individual cycles subdivided by a > - shaped separation area.

The LSM spectral decomposition is based on the definition of a two dimensional wave model

$$P(Y,L)\sim Acos({\omega }_{\mathrm{Y}}\mathit{Y}+{\omega }_{\mathrm{L}}\mathit{L}+{\tau }^{*})+Bsin({\omega }_{\mathrm{Y}}\mathit{Y}+{\omega }_{\mathrm{L}}\mathit{L}+{\tau }^{*})+ϵ,$$

(21)

where P is given by the assigned patch polarity. To show the advantage and robustness of the LSM the raw data is employed without any further preprocessing. The selected frequencies ${\omega }_{\mathrm{Y},n},{\omega }_{\mathrm{L},n}\propto n^{-1}$ are given on a rectangular grid with inverse distance so that the periods are distributed uniformly alongside the consecutive counter n. Furthermore, the frequency resolution close to $|{\omega }_{\mathrm{L}}|=0$ is selected finer to better resolve this region.

Figure 6 gives the resulting LSM spectral decomposition in a reciprocal log-scale plot to analyse the different periods. Each of the selected peaks (numbered black dots) provides a FAP value of $p<{10}^{-10}$, meaning that these periods are significantly different from noise level. Each of these points originate from a local frequency refinement to achieve the local maximum amplitude. The binary order information $\{-1,1\}$ induced by the Sunspots, introduces higher harmonics into the spectrum, which do also appear with a significantly low FAP value.

Figure 6 provides good agreement between the peaks found and the common known periods. Table 1 summarizes the major outcome in comparison with literature. For example the $22\mathrm{y}\mathrm{r}\mathrm{s}$ Hale cycle varies from $18\dots 28\mathrm{y}\mathrm{r}\mathrm{s}$ (Usoskin [44]) which is identified by the two main peaks. Since the solar cycle is modulated, refer to Hathaway [45], it is natural that a broad spectrum with many harmonics becomes present. These subsequent patterns are related to a set of local maxima mainly collapsing on a horizontal line at $\left|f_{\mathrm{Lat}}\right|\approx 1.4·{10}^{-2}\mathrm{d}\mathrm{e}{\mathrm{g}}^{-1}$, indicating that the complex structure of the individual wing-shapes of the butterfly diagram with different widths, heights and orientation, corresponds to a main period $P_{\mathrm{Hale}}=21.634\mathrm{y}\mathrm{r}\mathrm{s}$. Next to that we can identify typical periods which are related to the Gleissberg process (see Table 1) or the short periods in the range of $P\approx 7.2\mathrm{y}\mathrm{r}\mathrm{s}$, which are consistent with the data provided by Prestes et al. [46], Kane [47] and partly with Deng et al. [48].

The latitudinal dimension of the spectrum spreads the peaks vertically which is a great advantage compared to a purely 1D analysis, where all signal would be projected on the $f_{\mathrm{Lat}}=0$ line. In addition, complex patterns such as lines of constant phase velocities (dashed curve of Figure 6) may be used to analyse modulations of the dynamics of the Sunspot motion.

6. Conclusions

In the present work a multidimensional extension of the Lomb-Scargle method is developed. The key aspect is the redefinition of phase argument to ${\overrightarrow{\varphi }}_{\mathrm{new}}=\overrightarrow{\omega }·\overrightarrow{t}-{\tau }^{*}$. We suggest using a modified shifting parameter ${\tau }^{*}$ in contrast to the traditional approach which shifts the ordinate ${\varphi }_{\mathrm{orig}}=\omega ·\left(t-\tau \right)$ instead of the phase. This enables multivariate modeling with one single scalar value ${\tau }^{*}$ for all independent variables as there is always a shifting parameter ${\tau }^{*}$ for which ${\int }_a^{b}sin(\varphi (t))cos(\varphi (t))\mathrm{d}t$ vanishes on any interval $[a,b]$.

With respect to the evaluation of measurement results, the noise rejection ${ϵ}_{\mathrm{LS}}\propto 4N^{-0.5}/\pi$ and confidence intervals ($\Delta a_k$ and $\Delta b_k$) of model parameters $a_k$ and $b_k$ are delineated. To emphasize the advantages of the LSM the traditional Fourier mode decomposition is compared. It turned out that the systematic error ${ϵ}_{\mathrm{T}}$ does not vanish for FT based methods with increasing number of samples on a fixed interval T, which finally lead to the common leakage effect. Here the signal amplitudes are distributed on an area of neighboring spectral components, i. e. individual bins or pixels. We conclude that the standard orthogonal mode related procedures do not represent a consistent estimator for model parameters $a_k$ and $b_k$. Whereas the introduction of ${\tau }^{*}$ in LSM leads to a consistent estimator with ${ϵ}_{\mathrm{T}}=0$ even for higher dimensions. Finally, it was shown that LSM converges to the true model parameters with increasing number of samples and provides a better noise rejection (${ϵ}_{\mathrm{LS}}<{ϵ}_{\mathrm{FS}}$) as well.

The examples from Section 5 underline the strengths of the developed procedure in a consecutive way. First the application on ideal two dimensional test data shows the ability to analyze fragmented time series. Here the sampling remains regular meaning $t=n{T}_0$ with $n\in \mathbb{N}$, but with $n_i-n_{i+1}\ne 1$ as an incomplete set of locations to describe missing values. A second quasi similar situation is given in the experimental data set of UDV measurements, expect a certain jitter. The dimensionality in this example was extended to ${\mathbb{R}}^3$ to decompose the wave parameters in frequency f, symmetry m and spatial frequency k.

As the third analysis, the Sunspot data sets represent the most general case of non-evenly sampling. Taking only the time series of positions into account, LSM is able calculate the spectrum on a individual frequency grid pronouncing the low frequencies close to zero. The assignment of $\{-1,1\}$ to each data point is a minor modification so the data can be assumed as binary raw data. The spectral decomposition with LSM shows the characteristic footprint of waves and its dynamics present in the butterfly diagram. The time frequency spectrum itself provides the commonly known frequencies, i. e. Hale-Cycle or Gleissberg process. The second frequency $f_{\mathrm{Lat}}$ provides information about the minor peaks which are spread out into the latitudinal domain. Moreover, from these values the Sunspot drift motion may be calculated using a characteristic v-line with its side bands.