## Abstract

Rotary analysis decomposes vector motions on the plane into counter-rotating components, which have proved particularly useful in the study of geophysical flows influenced by the rotation of the Earth. For stationary random signals, the motion at any frequency takes the form of a random ellipse. Although there are numerous applications of rotary analysis, relatively little attention has been paid to the statistical properties of the random ellipses or to the estimated rotary coefficient, which measures the tendency to rotate counterclockwise or clockwise. The precise statistical structure of the ellipses is reviewed, including the random behaviour of the ellipse orientation, aspect ratio and intensity. Special attention is then paid to spectral matrix estimation from physical data and to hypothesis testing and confidence intervals computed using the estimated matrices.

## 1. Introduction

‘As theories attempt to explain more details of the motion of the atmosphere and ocean, it becomes more important to make rigorous comparisons between theory and observations. Because physical phenomena vary as a function of frequency, measurements are usually described in the frequency domain where different theories apply in different frequency ranges’ [1], p. 627. Rotary component analysis, discussed from a statistical perspective in this paper, is a tool for aiding such studies. In oceanographic studies, rotary components have proved useful for investigating currents over topography, inertial motions and shelf waves [2]. Frequency-dependent quantities such as the rotary coefficient—which measures the tendency to rotate clockwise or counterclockwise—may be extracted, and moreover many of the quantities discussed are invariant to rotation of the coordinate system.

To study vector time series such as current and wind, the data are resolved into orthogonal components. Oceanographic currents are resolved into an eastward (zonal) component {*X*_{t}}, and a northward (meridional) component {*Y* _{t}}[2], and wind vectors are constructed similarly in meteorology [3]. (The eastward and northward plane representation of the variation of a velocity vector has often been examined in the scientific literature under the name hodograph or hodogram [1].) We can of course construct the discrete complex-valued random process where *X*_{t}=Re(*Z*_{t}) and *Y* _{t}=Im(*Z*_{t}), so we have *Z*_{t}=*X*_{t}+i*Y* _{t}. For such vector time series, the relationships between the two real-valued components lead to structures on the complex plane, and geometrical information is encoded in these structures.

The rotary analysis method [2,4–8] decomposes motions on the complex plane into counter-rotating components, which have proved particularly useful in the study of geophysical flows influenced by the rotation of the Earth. For example, inertial waves [5] in oceanography have clockwise/counterclockwise polarization in the Northern/Southern Hemisphere, and because the spectrum of the complex-valued time series is asymmetric about zero frequency, it is possible to readily distinguish between clockwise and counterclockwise effects. This is illustrated in figure 1 for current data (centimetres per second) recorded in the Labrador Sea [9–11]. (The Labrador Sea observations analysed here are available in the electronic supplementary material.) For the real-valued series {*X*_{t}} and {*Y* _{t}}, figure 1*a* gives the necessarily symmetric power spectra; because of the symmetry about frequency *f*=0, only positive frequencies are shown. For the complex-valued series {*Z*_{t}}, figure 1*b* shows the clockwise (negative frequencies, overplotted on the positive frequency axis) and counterclockwise (positive frequencies) spectra. While the spectra for the real series are extremely similar, there are easily discernible differences between the clockwise and counterclockwise spectra, e.g. the vertical dotted line marks the inertial frequency strongly represented in the clockwise spectrum, but weakly so in the counterclockwise one. Figure 1*c* shows the rotary coefficient (discussed in detail in §3*f*) corresponding to figure 1*b*, and a value of almost exactly −1 is seen at the inertial frequency, in line with the expected clockwise rotation at the inertial frequency in the Northern Hemisphere.

At any frequency, the motion on the complex plane takes the form of an ellipse. Under the assumption that {*Z*_{t}} is a (second-order stationarity) random process, this ellipse is random in the sense that properties such as its orientation, aspect ratio and size will vary from realization to realization. The precise statistical structure of the ellipses is reviewed, including the random behaviour of the ellipse orientation, aspect ratio and intensity. (If, for every realization, the ellipse has a fixed orientation, aspect ratio and sense in which it is described, the process is completely polarized at that frequency.) Special attention is then paid to spectral matrix estimation from physical data and to hypothesis testing and confidence intervals (CIs) computed using the estimated matrices. A (non-exhaustive) glossary of notation and terminology is provided in table 1.

## 2. Rotary components model

Consider two deterministic series and with *a*_{x},*a*_{y} being the amplitudes, *f*_{0}>0 the frequency, *t* time and *θ*_{x},*θ*_{y} the phases. Then the complex quantity *z*_{t}=*x*_{t}+i*y*_{t} can be written as
i.e. as an ellipse that is the sum of a counterclockwise phasor, and a clockwise rotating phasor, where a rotating phasor is of the form *a*e^{i(2πft+θ)}. (Practically, a negative-frequency complex exponential can be seen as a counterclockwise rotation evolving in the opposite time direction—or just a clockwise rotation.) The phasors are called the *rotary components* [2].

We now consider the case of discrete stationary random signals. Without loss of generality, we shall assume all series have zero mean. The sample interval is Δ_{t} and the Nyquist frequency is *f*_{𝒩}=1/(2Δ_{t}). A discrete complex-valued process {*Z*_{t}} is said to be second-order stationary (SOS) if the real-valued series {*X*_{t}} and {*Y* _{t}} are jointly stationary stochastic processes, so that {*X*_{t}} and {*Y* _{t}} are individually SOS processes, and is a function of *τ* only for all *t*. The covariance of *Z*_{t+τ} and *Z*_{t} is We also define the relation between *Z*_{t+τ} and *Z*_{t} as Now,
so that SOS means that both and are functions of *τ* only [12] whence we obtain the autocovariance sequence (ACVS) with and the autorelation sequence (ARS) with (Schreier & Scharf [13] use the terminology ‘wide-sense stationary’ rather than SOS.) We note that and *r*_{Z,τ}=*r*_{Z,−τ}, i.e. the ACVS is complex Hermitian, while the ARS is complex symmetric. If i.e. *Z*_{t+τ} is uncorrelated with for all *τ*, then the process {*Z*_{t}} is said to be *proper* [13].

An SOS signal *Z*_{t}=*X*_{t}+i*Y* _{t} is then decomposed in terms of random ellipses at different frequencies, each random ellipse again being considered the sum of a counterclockwise and a clockwise rotating phasor. To be precise, because {*Z*_{t}} is SOS with zero mean, it has the spectral representation
2.1for all integers *t*, where the above equality is in the mean square sense. It has been shown that {*Z*(*f*)} has the following properties [14], p. 180:

—

*E*{d*Z*(*f*)}=0 for all |*f*|≤*f*_{𝒩};— for frequencies

*f*and*f*′ contained in the interval 2.2where*F*_{1}(*f*) is a real non-decreasing function called the spectral distribution function (SDF), and*δ*(⋅) is the Dirac delta function; so the covariance is zero unless*f*=*f*′; and— in general, d

*Z*(*f*)≠d*Z**(−*f*).

At frequencies *f* and −*f* (i.e. at frequency magnitude |*f*|), the contribution to {*Z*_{t}} is
2.3Thus *Z*_{t}(*f*) is formed from random increments, d*Z*(*f*) and d*Z*(−*f*), and is the parametric equation of a *random* ellipse, comprising the addition of two oppositely rotating motions with random amplitudes and phases, the rotary components.

Figure 2 illustrates, for Δ_{t}=1, the random ellipse that is centred at the origin, with a full rotation being performed with a period of 1/*f*. Because *t* is discrete, *Z*_{t}(*f*) only describes a subset of the points of the ellipse, 𝒵_{t}(*f*) say (shown by the elliptical curve), that would be traced if *t* were continuous. This subset is finite if and only if *f* is rational; for example, if then *Z*_{t}(*f*) takes only 10 values (points in the complex plane) as *t* traverses the integers. These points are shown by the black dots in figure 2, which also illustrates the geometrical properties of the random ellipse 𝒵_{t}(*f*). Here *A*(*f*)≡|d*Z*(*f*)|+|d*Z*(−*f*)| and *B*(*f*)≡||d*Z*(*f*)|−|d*Z*(−*f*)|| are the lengths of the semi-major and semi-minor axes, respectively; and *Θ*(*f*) is the azimuth (orientation), the angle that the major axis of the ellipse makes with the *x*-direction,
2.4The (signed) aspect ratio is ℰ(*f*)=±*B*(*f*)/*A*(*f*), where the sign is sgn{|d*Z*(*f*)|−|d*Z*(−*f*)|}, and is positive/negative for counterclockwise/clockwise rotation, i.e.
2.5A measure of the size of the ellipse is its ‘intensity’ [15], p. 78
2.6The statistics (2.4)–(2.6), which describe the geometrical properties of the ellipse, depend on the random variables (RVs) d*Z*(*f*) and d*Z*(−*f*).

Given a realization of the process {*Z*_{t}}, the six types of motion resulting from the realized values of the RVs d*Z*(*f*) and d*Z*(−*f*) are given in table 2.

## 3. Complex-valued second-order stationary processes and spectral matrices

### (a) Covariances

It will be useful to define an augmented process for which, under SOS, is given by
and the sequence of matrices captures all the second-order statistics of {*Z*_{t}}. In terms of the constituent real-valued jointly stationary processes, we define a process *V*_{t}=[*X*_{t},*Y* _{t}]^{T}, for which, under SOS, is given by
where , , , Now *U*_{t}=*T**V*_{t} where ** T**=

*M*^{−1}with [3] and so and

*s*_{V,τ}=

*M*

*s*_{U,τ}

*M*^{H}.

### (b) Spectra

Yaglom [14], p. 311 considers the covariance of stationary complex-valued vector sequences; it follows that {*s*_{Z,τ}} and {*r*_{Z,τ}} have Fourier–Stieltjes representations
where *F*_{1}(⋅) is as in (2.2), and *F*_{2}(⋅) is a complex-valued function of bounded variation, called the cross-spectral distribution function of {*Z*_{t}} and If {*s*_{Z,τ}} and {*r*_{Z,τ}} are absolutely summable, then they can be represented by ordinary Fourier integrals,
3.1where *S*_{Z}(*f*)d*f*=d*F*_{1}(*f*) and *R*_{Z}(*f*)d*f*=d*F*_{2}(*f*). Because *R*_{Z}(*f*) is complex-valued, we can write *R*_{Z}(*f*)=|*R*_{Z}(*f*)|e^{iω(f)}, where . We have that *R*_{Z}(*f*) is symmetric about *f*=0 because the ARS is symmetric.

The ACVS can be written as
In order that the right-hand side depends only on *τ*, we require *E*{d*Z*(*f*) d*Z**(*f*′)} to be non-zero only when *f*=*f*′. Combining with the first part of (3.1) gives
3.2

The ARS can be written as
In order that the right-hand side depends only on *τ*, we need *E*{d*Z*(*f*) d*Z*(*f*′)} to be non-zero only when *f*=−*f*′. Combining with the second part of (3.1) gives
3.3(Notice that if so that the process is proper, then *R*_{Z}(⋅)=0 and in this case )

Now let d** U**(

*f*)=[d

*Z*(

*f*),d

*Z**(−

*f*)]

^{T}for |

*f*|>0. Then from (3.2) and (3.3), we obtain 3.4where 3.5Now

*E*{d

**(**

*U**f*)d

*U*^{T}(

*f*)}=

*E*{d

**(**

*U**f*)(d

***(**

*U**f*))

^{H}} so that the left of (3.4) means that d

**(**

*U**f*) is uncorrelated with its complex conjugate d

***(**

*U**f*). A complex-valued

*vector*that is uncorrelated with its complex conjugate is said to be a

*proper complex-valued vector*—so d

**(**

*U**f*) is proper. (Notice that propriety of the

*vector*d

**(**

*U**f*) only involves the single frequency

*f*.) If {

*Z*

_{t}} is Gaussian (normal), then the fact that d

**(**

*U**f*) is proper means [16] that d

**(**

*U**f*) is jointly complex Gaussian with mean 0 and covariance matrix

*S*_{U}(

*f*) d

*f*, which we write as The probability density function (PDF) for this distribution is given in [17] (see also (4.5) later); without propriety, the situation is much more complicated [18].

We now consider the relationships between *S*_{Z}(⋅),*R*_{Z}(⋅) and the spectra and cross-spectra, *S*_{X}(⋅),*S*_{Y}(⋅),*S*_{XY}(⋅), of the real-valued components. For SOS complex-valued processes, the real-valued series {*X*_{t}} and {*Y* _{t}} are jointly stationary stochastic processes and the corresponding spectral matrix is given by
Now because {*s*_{U,τ}} and *S*_{U}(⋅) are a Fourier transform pair. This gives us two useful results:
3.6and
3.7

### (c) Conjugate coherence

Now d** U**(

*f*) is proper so that if we consider estimation of one component, d

*Z**(−

*f*), from the other, d

*Z*(

*f*), linear estimation is appropriate and the linear minimum mean square estimator takes the form [19]. We can define a form of (complex-valued) coherence as and the corresponding real-valued magnitude-squared coherence as 3.8If

*S*

_{Z}(

*f*)

*S*

_{Z}(−

*f*)=0, the coherence (3.8) does not exist, but note from table 2 that if

*S*

_{Z}(−

*f*)=0 then |d

*Z*(−

*f*)|≡0 with probability one, implying that

*all realizations*will be counterclockwise circular. Likewise, if

*S*

_{Z}(

*f*)=0, for all realizations the motion is clockwise circular. In both cases, the process may be said to be

*circularly polarized*at

*f*[20].

It is readily shown that *γ*^{2}_{*}(*f*) is bounded between 0 and 1, and measures the accuracy with which d*Z**(−*f*) can be linearly estimated from d*Z*(*f*). It also measures the magnitude of the linear correlation between the components of {*Z*_{t}} and its conjugate process, {*Z**_{t}}, at frequency *f*. Hence we adopt the nomenclature *conjugate coherence* for . When *γ*^{2}_{*}(*f*)=1, the mean square estimation error is zero, and there is convergence in probability, so that, with probability one,
3.9

A statistical study of the random nature of the elliptical paths (2.3) was carried out in Rubin-Delanchy & Walden [16], who showed that the conjugate coherence plays a key role. When *γ*^{2}_{*}(*f*)=1, the use of (2.4) and (3.9) gives for the azimuth (orientation)
3.10while for the aspect ratio of (2.5)
but because so that
3.11We note for future reference that *S*_{U}(*f*) becomes a singular matrix, given by
3.12Hence both the azimuth and aspect ratio are fixed, with values *θ*(*f*) and *ϵ*(*f*), respectively, for the particular frequency *f*; the only aspect of the random ellipse that is random is its size. The intensity of the ellipse is given by the RV *I*(*f*)=*A*^{2}(*f*)+*B*^{2}(*f*)=2[|d*Z*(*f*)|^{2}+|d*Z*(−*f*)|^{2}]; so the mean of the intensity is *E*{*I*(*f*)}=2[*S*_{Z}(*f*)+*S*_{Z}(−*f*)]d*f*.

Taking moduli in (3.9), |d*Z*(−*f*)|=|*R*_{Z}(*f*)/*S*_{Z}(*f*)||d*Z*(*f*)|, and from table 2, we see that motion is counterclockwise elliptical if |*R*_{Z}(*f*)|<*S*_{Z}(*f*) and clockwise elliptical if |*R*_{Z}(*f*)|>*S*_{Z}(*f*). The motion is rectilinear if |*R*_{Z}(*f*)|=*S*_{Z}(*f*), and because this implies *S*_{Z}(*f*)=*S*_{Z}(−*f*). So when *γ*^{2}_{*}(*f*)=1, for every realization the ellipse has a fixed orientation, ellipticity and sense in which the ellipse is described.

### (d) Ellipse decompositions

In the general case, when is not necessarily unity, it was shown in [16] that we can write the ellipse *Z*_{t}(*f*) in the form *Z*_{t}(*f*)=*P*_{t}(*f*)+*Q*_{t}(*f*), where *P*_{t}(*f*) and *Q*_{t}(*f*) are both parametric equations for random ellipses, with a full rotation having a period 1/*f*, centred at the origin. Moreover, *S*_{U}(*f*) in (3.5) can be written as
3.13where and are Hermitian matrices with zero determinant, and
where *C*_{1}=[1+*γ*_{*}(*f*)]/2 and *C*_{2}=[1−*γ*_{*}(*f*)]/2, and *γ*_{*}(*f*) is the positive square root of *γ*^{2}_{*}(*f*). The parametric equation *P*_{t}(*f*) is called the *singular* ellipse, its azimuth is fixed with value *θ*(*f*) and its aspect ratio is fixed with value *ϵ*(*f*). The parametric equation *Q*_{t}(*f*) also has fixed azimuth and aspect ratio. While its aspect ratio again has value *ϵ*(*f*), its azimuth has value *θ*(*f*)+*π*/2 and so *Q*_{t} is called the *orthogonal* ellipse. We note that the mean intensity of *Z*_{t}(*f*) is the sum of the mean intensities of *P*_{t}(*f*) and *Q*_{t}(*f*) and is *E*{*I*(*f*)}=2[*S*_{Z}(*f*)+*S*_{Z}(−*f*)]d*f*, just the same as when *γ*^{2}_{*}(*f*)=1.

The conjugate coherence *γ*^{2}_{*}(*f*) controls the relative sizes of the two ellipses. When it is unity, *Z*_{t}(*f*) reduces to *P*_{t}(*f*), but as *Q*_{t}(*f*) grows towards making the two ellipses have equal influence. The ratio of the mean intensity of *P*_{t}(*f*) to that of *Z*_{t}(*f*) is
3.14

The oceanographer Mooers [7] noted that some of the rotary analysis concepts ‘have analogs in the study of polarized and partially coherent light. Though Born and Wolf [21] did not exploit the concept of negative frequencies … some of the analysis methods used in optics may aid the oceanographic analyst’. In particular, there exists another decomposition for *S*_{U}(*f*) as [21]
3.15where
with and Now let [22]
where and are two independent proper vectors with covariance matrices and respectively. Then the ellipse that arises from summing and *E*_{t}(*f*) say, has the form
3.16and [*E*_{1}(*f*),*E*_{2}(*f*)]^{T} is proper with covariance *S*_{U}(*f*), whereas [d*Z*(*f*),d*Z**(−*f*)]^{T} has covariance *S*_{U}(*f*)d*f*. Under Gaussianity, first- and second-order moments define the entire statistical structure; so under Gaussianity, the statistical properties of *E*_{t}(*f*) and *Z*_{t}(*f*) are identical up to a scaling factor, i.e. the statistical properties of *E*_{t}(*f*) and *Z*_{t}(*f*) [d*f*]^{−1/2} are the same.

The elements of and may be written [21], p. 551 as Because is singular, is an ellipse with fixed azimuth and aspect ratio, and by substituting the appropriate elements in (3.10) and (3.11), we get By contrast, is made up of two oppositely rotating circular motions—independent under Gaussianity—with the same mean amplitude; so its azimuth and aspect ratio should be maximally variable.

The ratio of the mean intensity of to that of *E*_{t}(*f*) is
3.17i.e. *P*(*f*) gives a measure of the departure of *E*_{t}(*f*) (and hence *Z*_{t}(*f*) [d*f*]^{−1/2}) from the singular ellipse. In the study of Born & Wolf [21], the quantity corresponding to the right-hand side of (3.17), i.e. *P*(*f*), is called the degree of polarization, 0≤*P*(*f*)≤1.

The matrix is unitary, so from (3.13),
3.18i.e. where and Because is unitarily similar to *S*_{U}(*f*) it has the same trace, determinant and eigenvalues. Analogous results hold for the two terms on the right of (3.18). So and each have determinant zero and the unitary similarity transform preserves the Hermitian nature of the matrices. The ratio of intensities in (3.14) can also be written as
3.19Likewise we get from (3.15)
3.20where and The unitary similarity transform means that will have zero determinant and Hermitian form like and will be diagonal like The ratio of intensities in (3.17) may therefore be written as
3.21

Comparing (3.19) and (3.21), we see that *γ*^{2}_{*}(*f*)=1 corresponds to *P*(*f*)=1 [3]. Because when the former is true, every realization of the ellipse has a fixed orientation, ellipticity and sense in which the ellipse is described, this behaviour corresponds to ‘complete polarization’ at *f*.

In [15,21], the underlying model is that of a partially polarized quasi-monochromatic plane wave field, and relevant statistical discussion is given in e.g. [15,23]. In [15], analogous decompositions to those presented here are given in terms of the Stokes parameters.

### (e) Rotational invariants

The intensity ratios (3.19) and (3.21) are defined in terms of spectral rotational invariants. To see this more clearly, consider rotating the coordinate system clockwise through a deterministic angle *φ*. The time-series components, referred to the new axes, are given by
where ** G**(

*φ*) is the planar rotation matrix. Then

*S*_{V}(

*f*;

*φ*)=

**(**

*G**φ*)

*S*_{V}(

*f*)

*G*^{T}(

*φ*), where

*S*_{V}(

*f*;

*φ*) is the spectral matrix for [

*X*

_{t}(

*φ*),

*Y*

_{t}(

*φ*)]

^{T}. But

*S*_{V}(

*f*;

*φ*) and

*S*_{V}(

*f*) are unitarily similar matrices; so the determinant and trace and hence

*P*(

*f*) in (3.21) are unchanged. Then the spectra of the rotated time series are where

*S*

_{XX}(

*f*;

*φ*) and

*S*

_{YY}(

*f*;

*φ*) are the SDFs corresponding to

*X*

_{t}(

*φ*) and

*Y*

_{t}(

*φ*), respectively. Because tr{

*S*_{V}(

*f*;

*φ*)}=

*S*

_{XX}(

*f*;

*φ*)+

*S*

_{YY}(

*f*;

*φ*)=

*S*

_{XX}(

*f*)+

*S*

_{YY}(

*f*), we see that indeed the trace is invariant to rotation. Note however that the power spectra are not invariant, and consequently the magnitude-squared coherence of {

*X*

_{t}} and {

*Y*

_{t}} is not invariant, i.e. for all

*φ*. If we define

*Z*

_{t}(

*φ*)=

*X*

_{t}(

*φ*)+i

*Y*

_{t}(

*φ*), then

*Z*

_{t}(

*φ*)=

*Z*

_{t}e

^{−iφ}, and so The ACVS is unchanged by the rotation,

*S*

_{Z}(

*f*;

*φ*)=

*S*

_{Z}(

*f*), where

*S*

_{Z}(

*f*;

*φ*) is the Fourier transform of {

*s*

_{Z,τ}(

*φ*)}. Furthermore,

*R*

_{Z}(

*f*;

*φ*)=

*R*

_{Z}(

*f*)e

^{−2iφ}, where

*R*

_{Z}(

*f*;

*φ*) is the Fourier transform of {

*r*

_{Z,τ}(

*φ*)}. Consequently, for all

*φ*, i.e. unlike the ordinary coherence, the conjugate coherence, and hence (3.19), is invariant to a rotation.

### (f) Rotary coefficient

The rotary coefficient is defined as [2], p. 431 *ρ*(*f*)=[*S*_{Z}(*f*)−*S*_{Z}(−*f*)]/[*S*_{Z}(*f*)+*S*_{Z}(−*f*)]. From (3.6), *ρ*(*f*) can also be written as
3.22where we have used the fact that The rotary coefficient satisfies −1≤*ρ*(*f*)≤1 and measures the tendency to rotate in a counterclockwise or clockwise manner. It provides an objective means of quantifying the rotation associated with the asymmetry of the spectrum for {*Z*_{t}}. Traditionally, a positive frequency is associated with counterclockwise rotation. Let *f*>0. Then if *ρ*(*f*)=+1 (i.e. *S*_{Z}(−*f*)=0), motion is all counterclockwise at that frequency, whereas if *ρ*(*f*)=−1 (i.e. *S*_{Z}(*f*)=0), motion is all clockwise at that frequency, and *ρ*(*f*)=0 for rectilinear motion (unidirectional flow).

The degree of polarization, rotary coefficient and conjugate coherence are all simply related. From (3.8) and (3.17), we get
3.23see e.g. [9]. Because *P*^{2}(*f*) and *γ*^{2}_{*}(*f*) are invariant to coordinate rotation, so must be *ρ*^{2}(*f*). Furthermore, we see that *P*^{2}(*f*)≥*ρ*^{2}(*f*) and *P*^{2}(*f*)≥*γ*^{2}_{*}(*f*).

## 4. Statistical properties

Here we assume that {*Z*_{t}} is complex Gaussian (normal), i.e. for all *n*≥1 and any the joint distribution of the real and imaginary components *X*_{t0},*Y* _{t0},*X*_{t1},*Y* _{t1},…,*X*_{tn−1},*Y* _{tn−1} is multivariate Gaussian.

### (a) Ellipse parameters

With {*Z*_{t}} being complex Gaussian, we know that [d*Z*(*f*),d*Z**(−*f*)]^{T} is jointly complex Gaussian [16]. Using this, the distribution of the aspect ratio—the statistic ℰ(*f*) in (2.5)—was given by [16], parametrized by the degree of polarization and the rotary coefficient; these are affected by *γ*^{2}_{*}(*f*) through (3.23). For , its PDF takes the form
4.1(where we have not explicitly included the dependence on frequency *f*). When *ρ*(*f*) is positive/negative the mass of the PDF is to the right/left, and the distribution is skewed to the left/right, with skewness diminishing as The distribution is symmetric when *ρ*(*f*)=0, and then the mean aspect ratio is zero, corresponding to rectilinear motion.

The azimuth, the statistic *Θ*(*f*) in (2.4), has an angular distribution. For its PDF was given for general *γ*^{2}_{*}(*f*) in [16] as
4.2where and *γ*^{2}_{*}<1. (In (4.2) _{2}*F*_{1}(*α*_{1},*α*_{2};*α*_{3};*z*) is the hypergeometric function with two and one parameters, *α*_{1},*α*_{2} and *α*_{3}, and scalar argument *z* [24].) In [16], it was verified that *θ*(*f*) in (3.10) gives the mean orientation. By combining (3.7) and (3.10), we see that *θ*(*f*) satisfies
4.3where we use the four-quadrant inverse tangent. This equation is the form often quoted in the literature, without formal justification prior to [16], for the mean orientation [2], p. 496.

Suppose we wish to illustrate possible elliptical paths, i.e. simulate a realization {*z*_{t}(*f*)} of the contribution {*Z*_{t}(*f*)}, from a realization {*z*_{t}} of {*Z*_{t}}. This was considered in [16] where it was pointed out that, since *Z*_{t}(*f*)=*O*_{p}([d*f*]^{1/2}), where *O*_{p} denotes order in probability, a rescaling will be necessary for physical implementation. One approach is to work with *Z*_{t}(*f*)[d*f*]^{−1/2} instead of *Z*_{t}(*f*) itself. We can write
4.4where and Of course (4.4) is of the same form as (3.16). The intensity for the ellipse defined by (4.4) is *I*′(*f*), say, given by Now, and ; so *E*{*I*′(*f*)}=2[*S*_{Z}(*f*)+*S*_{Z}(−*f*)]=2 tr{*S*_{U}(*f*)}.

Let [*z*_{1},*z*_{2}]^{T} be a realization of Because is jointly complex Gaussian with mean zero and covariance matrix *S*_{U}(*f*), its PDF is by definition [17]
4.5where Consequently, the PDF of *I*′(*f*)/2 follows from [25], eqn 5.19, and adjusting for the factor of 2 we find, for the PDF of *I*′(*f*) to be
where *C*_{0}=1/[2*P*(*f*)tr{*S*_{U}(*f*)}].

For {*Z*_{t}}, a complex-valued autoregressive process of order 2 (CAR(2)) was used; full details of simulating from *Z*_{t}(*f*)[d*f*]^{−1/2} are given in [16]. The first row of figure 3 gives three independent realizations of the elliptical paths *Z*_{t}(*f*)[d*f*]^{−1/2} at when *γ*^{2}_{*}(*f*_{0})=1. In this case, the only aspect of the random ellipse that is random is its size. This is what we see: the azimuth (orientation) and aspect ratio are the same in all three realizations. The orientation is given by (3.10) (dashed and dotted lines overlaying). By way of contrast, the second row of figure 3 gives three independent realizations when *γ*^{2}_{*}(*f*_{0})=0.5. Here the size, aspect ratio and orientation are different in all three cases. For *γ*^{2}_{*}(*f*_{0})=0.5, the sample mean of intensity realizations 2[|*z*_{1}|^{2}+|*z*_{2}|^{2}] over 2000 simulated ellipses was found to be 881.01 compared with the theoretical value 2 tr{*S*_{U}(*f*_{0})}=880.97, a remarkably good match. (To determine the theoretical value 2 tr{*S*_{U}(*f*_{0})}, the values of *S*_{Z}(*f*_{0}) and *S*_{Z}(−*f*_{0}) for the CAR(2) process were used in (3.5).)

### (b) Spectral estimators and Wishart distributions

We have discussed the distributions of ellipse statistics, *Θ*(*f*),ℰ(*f*) (azimuth, aspect ratio), expressed in terms of d*Z*(*f*) and d*Z*(−*f*), and of *I*′(*f*) (rescaled intensity) expressed in terms of and It is important to also draw statistical inferences on quantities such as *ρ*(*f*). Here we choose to use multitaper estimators for the spectral matrices.

We denote by {*h*_{k,t},*t*=0,…,*N*−1} the *k*th (*k*=0,…,*K*−1) real-valued orthonormal taper, perhaps the sine tapers [26]. Suppose we are given a sample of size *N*, i.e. *U*_{0},…,*U*_{N−1}. We form the product {*h*_{k,t}*U*_{t}} of the *k*th real-valued taper with *U*_{t}, and then compute its Fourier transform
4.6where As an estimator of *S*_{U}(*f*), we take
4.7an estimator with *K* complex degrees of freedom. (The use of in (4.6) ensures, for example, that has the correct form for a power spectrum estimator [27].)

Now *U*_{t}=*T**V*_{t}, so that and an alternative form for the spectral matrix estimator is
where
4.8Under Gaussianity, {*J*_{V,k}(*f*),*k*=0,…,*K*−1} are independently and identically distributed with a complex bivariate Gaussian distribution with mean **0** and covariance matrix *S*_{V}(*f*) [9]
4.9Here 2*W*_{N} is the width of the spectral window induced by tapering. For sine tapers *W*_{N}=(*K*+1)/[2(*N*+1)Δ_{t}] [26], which decreases to zero as for a fixed *K*. Also, the frequency band within which the overall spectral window due to tapering [26] is concentrated must be narrow enough that the components of *S*_{V}(*f*) are essentially constant across it. Likewise, since *J*_{U,k}(*f*)=*T**J*_{V,k}(*f*), and
4.10Finally, and and so provided *K*≥2 (to avoid matrix singularity), the distribution of is two-dimensional complex Wishart with *K* complex degrees of freedom and mean *K**S*_{V}(*f*) [17], which we denote by
4.11Similarly,
4.12Equations (4.9)–(4.12) provide the framework for the statistical results in the following subsections.

### (c) Test for null degree of polarization

If *S*_{U}(*f*) in (3.15) is of the form then *P*(*f*)=0 so that [28], p. 411
4.13Here is an unknown scalar quantity, and *I*_{2} is the 2×2 identity matrix. So a test for perfectly unpolarized components at frequency *f* (*P*(*f*)=0) is equivalent to a test for sphericity for *S*_{U}(*f*), i.e. that d*Z*(*f*) and d*Z**(−*f*) have equal variance and are uncorrelated. The hypothesis test is versus Let *Ω* be the class of all 2-vector-valued complex Gaussian populations, and *ω* be the subclass satisfying the hypothesis of sphericity. The likelihood ratio for the test is the quotient of the maximum of the likelihood for variation of parameters specifying *ω*, say, to the maximum of the likelihood for variation of parameters specifying *Ω*, i.e. Because the quotient is composed from non-negative functions, and *ω*⊂*Ω*, it follows that and hence 0≤*λ*≤1. Clearly, a too *small* value of *λ* would lead to rejection of *H*_{0}. The critical region for the likelihood ratio test is the set of points for which *λ*≤*c*, where *c* is selected so that the test has a specified size *α*. By taking a power of *λ*, a simpler test statistic is obtained [29]. Under *H*_{0} and (4.12), the distribution of *T*_{s}(*f*) is beta with parameters *K*−1 and . For 0<*x*<1, the PDF is
4.14where *B*(*a*,*b*) is the beta function. (The PDF may be deduced from [29], eqns (3.8) by noting that our Wishart distribution has *K* degrees of freedom equivalent to the *n* used in [29].) The hypothesis is rejected for too small values of *T*_{s}(*f*). It is interesting to see how small *T*_{s}(*f*) should be before we reject the null hypothesis when *K* is small. For example, when *K*=5, and the size of the test is *α*=0.05, the hypothesis is rejected if *T*_{s}(*f*) is less than the critical value of 0.40. Critical values derived from *g*(⋅) in (4.14) are shown in figure 4*a*.

### (d) Test for null rotary coefficient

From (3.22),
4.15Then the mean aspect ratio is zero, corresponding to rectilinear motion. Since (4.15) defines *S*_{V}(*f*) to be real-valued, we see that a test for *ρ*(*f*)=0 is equivalent to a test for *real structure* for *S*_{V}(*f*), i.e. whether *S*_{V}(*f*)=Re{*S*_{V}(*f*)}+i Im{*S*_{V}(*f*)} is real-valued [9]. For Gaussian processes, the likelihood ratio criterion for testing *H*_{0}:Im{*S*_{V}(*f*)}=**0** against *H*_{1}:Im{*S*_{V}(*f*)}≠**0** was derived in [30]. By taking a power of the likelihood ratio statistic, the test statistic is obtained [31]. Under *H*_{0} and (4.11), the distribution of *T*_{r}(*f*) is beta with parameters *K*−1 and For 0<*x*<1, the PDF is
4.16Again, and for the same reason as in §4*c*, the hypothesis is rejected for too *small* values of *T*_{r}(*f*). Critical values derived from *g*(⋅) in (4.16) are shown in figure 4*b*.

To illustrate this test, we consider magnetic-field data recorded by a spacecraft in the *Cluster* mission, an international solar physics experiment to collect data on various aspects of the Sun. Three components were recorded: *W* (oriented towards the Sun), *X* (parallel to the ecliptic plane) and *Y* (perpendicular to the ecliptic and *W* in a northward sense). A dominant feature of our February 2003 dataset is a quasi-periodic ultra-low-frequency wave in the solar magnetic field with period about 30 s [32], i.e. a dominant frequency of *f*_{0}≈0.033 Hz. Two complex-valued series were formed, *Z*_{1,t}=*W*_{t}+i*Y* _{t} and *Z*_{2,t}=*X*_{t}+i*Y* _{t}. The measurement unit is nanoteslas (nT), and the sample interval is Δ_{t}=2s, giving a Nyquist frequency of 1/(2Δ_{t})=0.25 Hz. When the data are narrow-band filtered around *f*_{0}, the time evolution of the two complex series is as in figure 5. The filtered {*Z*_{1,t}} series shows rectilinear motion and {*Z*_{2,t}} shows counterclockwise elliptical motion. The spectral matrix estimate was computed from the unfiltered data using *K*=4 tapers, and in this case, the 5 per cent critical point of *T*_{r} is 0.501 and the 1 per cent point is 0.304. Series {*Z*_{1,t}} gives a value for *T*_{r}(*f*_{0}) of 0.9926, which is much larger than either critical value; so *H*_{0} (*ρ*(*f*)=0) is accepted at *f*=*f*_{0}. This is consistent with the rectilinear motion seen in figure 5*a*,*c*. For {*Z*_{2,t}}, the statistic *T*_{r}(*f*_{0}) is 0.0939, much smaller than the 1 per cent point, and *H*_{0} is rejected, consistent with the elliptical motion of figure 5*b*,*d*.

### (e) Confidence intervals for the rotary coefficient

Figure 1 made use of ocean current speed and direction time series for one depth from a set of six (110, 760, 1260, 1760, 2510 and 3476 m) recorded at a mooring in the Labrador Sea [10,11]. We used *N*=1600 observations with a sampling interval of Δ_{t}=1 h. In the spectral analysis, *K*=10 tapers were applied for which *W*_{N}=0.0034 c h^{−1}. At *f*, denote the 100(1−*α*) per cent CI as [9]. The interval is random: it depends on the estimate of the rotary coefficient, with and as in (4.7). It also depends on the chosen quantities *K* and *α*, and which is unknown. As shown in [9], coverage results are still good if we replace by the approximately unbiased estimate where comprising the components of in (4.7). The right end of the interval, , is thus approximated by the value of *ρ* such that (where *F*(⋅) is the distribution function for ), which can be found simply using any standard zero-finding algorithm. The distribution function *F*(⋅) itself can be found by numerical integration from the PDF *g*(⋅) for
where *y*_{ρx}=[(1−*ρ*)(1+*x*)]/[(1+*ρ*)(1−*x*)] and |*x*|<1 [9], eqn (22). Likewise is approximated by finding the value of *ρ* for which

Figure 6 shows the estimated rotary coefficients (solid dots) and corresponding 95 per cent CIs (solid horizontal bars) for the six observation depths at the four given frequencies, leading up to the inertial frequency of approximately 0.07 c h^{−1} [4] for the mooring at 56.75^{°}N. We see that, at the inertial frequency, the rotary coefficient is very close to −1 at all depths, i.e. close to the ideal *theoretical* outcome at the local inertial frequency, and the CIs are all very narrow.

## Acknowledgements

The author is very grateful to Jon Lilly for the Labrador Sea data, Tim Horbury for the *Cluster* mission data, and the Guest Editor and a referee for many helpful suggestions that improved the manuscript.

## Footnotes

One contribution of 17 to a Discussion Meeting Issue ‘Signal processing and inference for the physical sciences’.

- © 2012 The Author(s) Published by the Royal Society. All rights reserved.