## Abstract

An ellipse describes the polarized part of a partially polarized quasi-monochromatic plane wave field. Its azimuth angle and aspect ratio are functions of the elements of the covariance matrix associated with the polarized part at a particular time instant. Given an ensemble of *K* independent samples at that time, the distributions of the estimators of these parameters are derived. The estimation is thus based on a sample ensemble at any time, and does not assume temporal stationarity. Additionally, the azimuth angle estimator has an angular distribution so that non-standard statistical methods are needed when deriving its mean and standard deviation.

## 1. Introduction

When a quasi-monochromatic (narrow-band) plane wave propagates in the *z*-direction of a Cartesian coordinate system, it is found that, in the *x*–*y* plane perpendicular to the direction of travel, the endpoint of the vector traces out an instantaneous ellipse whose shape changes continuously with time (e.g. Brosseau 1998); the wave is partially polarized. A wave whose ellipse maintains a constant azimuth angle and aspect ratio, but whose size varies, is called completely polarized, and if no regular pattern is exhibited the wave is unpolarized.

As well as being of fundamental importance in optics, polarization is of great interest in many other fields such as astronomy (Simmons & Stewart 1985), atmospheric science (Hayashi 1979), oceanography (Emery & Thomson 1998), geophysics (Park *et al*. 1987) and elsewhere (Schreier submitted). Parameters describing the state of polarization can be calculated from observations but these are only estimates of the true parameter values (e.g. Simmons & Stewart 1985). Modern data-acquisition systems can often record multiple (*K*, say) independent views of a propagating wave in the form of a set of discrete-time finite-length data series. For example, the presence of an ultra-low-frequency (ULF) wave in the solar magnetic field was captured simultaneously by *K*=4 spacecraft in the Cluster mission, an international solar physics experiment, in February 2003 (Archer *et al*. 2005); see figure 1 and §6. Assuming spatial homogeneity, these multiple samples can be used to reduce variance in the estimation of instantaneous polarization ellipse parameters, a procedure known as ensemble averaging. Considering the instantaneous ellipse describing the polarized portion of the signal, we shall look at estimators of its *azimuth* angle, the angle which the major axis of the ellipse makes with the *x*-direction, and its *aspect ratio*, the ratio of length of the minor axis to that of the major axis (signed according to right- or left-handedness). New statistical results are required for the estimation of polarization parameters from such data and are the subject of this paper.

A commonly used model for a partially polarized plane wave is the complex representation whereby *Z*_{1}(*t*) and *Z*_{2}(*t*) say, are the analytical signal representations of the *x* and *y* components of the field (Barakat 1985): , *j*=1, 2, where _{j}(*t*) is the real field component and is its Hilbert transform. The real-valued stochastic processes _{1}(*t*) and _{2}(*t*) are taken to be zero mean and Gaussian-distributed, justified by the central limit theorem (Brosseau 1998, §3.3). Let , *k*=0, …, *K*−1 denote the *K* independent complex-valued observations (‘random samples’) of at an arbitrary time, where the superscript ‘T’ denotes transpose. Associated with these observations is an instantaneous sample Hermitian covariance matrix of the form(1.1)where superscript ‘*’ denotes complex conjugate; superscript ‘H’ denotes Hermitian (complex-conjugate) transpose; and the samples have zero mean and positive-definite covariance matrix . Since is constructed from *K* independent complex-valued vectors, the parameter *K* is also the number of complex degrees of freedom in , and .

Considering the partially polarized quasi-monochromatic wave field as an incoherent superposition of a fully polarized monochromatic wave field and a completely unpolarized wave field, can be uniquely decomposed as (Born & Wolf 1970)(1.2)respectively, where has a determinant of 0. The estimators are based on the sample covariance matrix of (1.1) at any arbitrary time and could thus be described as providing instantaneous estimates of the azimuth angle and aspect ratio. The results by Barakat (1985) are a special case of our results, since his results correspond to taking *K*=1; also, the azimuth angle estimator has an angular distribution so that we provide non-standard statistical methods for deriving its mean and standard deviation.

In §2 we carefully define the azimuth angle and aspect ratio of the ellipse describing the polarized portion of the signal. Some fundamental statistical background is provided in §3, and the forms of the parameter estimators are derived from the *K*-sample covariance matrix of (1.1). In §4 we derive the probability density function (PDF), mean and standard deviation of the azimuth estimator, taking into account that it has an angular distribution, and a simulation experiment is carried out showing that the standard definitions of these quantities are, by comparison, quite unsatisfactory. The PDF of the aspect ratio is discussed in §5. The ULF wave in the solar magnetic field is analysed in §6. Our results are summarized in §7.

## 2. Stokes and ellipse parameters

We write the matrices in (1.2) aswhere , and , where ‘det’ denotes determinant. The elements of these matrices are given by (Born & Wolf 1970)If , then , and the polarized part of is associated with elliptical motion in the original coordinates. If, further, , and then , and the polarization component is circular (Hayashi 1979), a special case of elliptical polarization.

Associated with the *partial polarization* covariance matrix, , are the usual Stokes parameters (Eliyahu 1994) given byFor notational simplicity, we denote , , , by *s*_{0},*s*_{1}, *s*_{2} and *s*_{3}.

The Stokes parameters associated with the *full polarization* covariance matrix, , areTherefore, while the other Stokes parameters are the same for and .

The azimuth angle, *ψ*, satisfies (Born & Wolf 1970, pp. 31, 555)(2.1)Given *s*_{2}/*s*_{1} no distinction can be made between a solution *ψ*_{0} and a solution *ψ*_{0}±*π*/2; there is not enough information to know which is which. Hence, as done by Eliyahu (1993) we in fact take *ψ* to be the angle of the major *or minor* axis of the ellipse from the *x*-direction, with −*π*/4≤*ψ*<*π*/4; see figure 2.

For circular polarization, both the numerator and the denominator of the ratio in (2.1) are 0, and the angle is undefined, as it should be since the idea of the orientation of an ellipse axis is meaningless in the case of a circle.

Let 2*a* and 2*b*(*a*≥*b*) be the lengths of the major and minor axes of the ellipse, as in figure 2, and *Χ*, (−*π*/4≤*Χ*<*π*/4) be the angle which characterizes the ellipticity and the sense in which the ellipse is being described. Then (Born & Wolf 1970, p. 555), tan *Χ*=±*b*/*a* (*Χ*≷0 according as the polarization is right- or left-handed) and in terms of the Stokes parameters(2.2)But for (−*π*/4≤*Χ*<*π*/4). If we define *y*=tan 2*Χ* and *ϵ*=tan *Χ*=±*b*/*a*, then(2.3)(2.4)The only solution for *ϵ* in (2.4) for which −1≤*ϵ*≤1 is(2.5)a sigmoidal curve. Using (2.3) and (2.5), we see that the (signed) aspect ratio, *ϵ*, is given by(2.6)

The degree of polarization is (Born & Wolf 1970, p. 555)(2.7)

## 3. Statistical fundamentals

### (a) Complex bivariate Gaussian distribution

is said to have the bivariate proper (or circular) complex Gaussian distribution with zero mean and covariance matrix , written if *Z*_{0} has PDF (Goodman 1963; Picinbono 1993, p. 122)(3.1)Under the Gaussian assumption on the real-valued stochastic processes _{1}(*t*) and _{2}(*t*) the random vector is obviously a bivariate complex Gaussian vector, but it is only proper and hence has probability density (3.1) if (Picinbono 1993, p. 120) . For this to be true, we need to show that and . These follow immediately from the results on Hilbert transform relationships given by Bendat & Piersol (1986, table 13.2, p. 499) and thus we conclude that has the bivariate proper (or circular) complex Gaussian distribution with PDF (3.1).

### (b) Complex Wishart distribution

Let where are *K* independent bivariate proper complex Gaussian random vectors, each with the distribution. When *K*≥2 the random 2×2 matrix is full rank and(3.2)that is, ** W** has the non-singular two-dimensional complex central Wishart distribution with

*K*complex degrees of freedom and mean . Goodman (1963) pointed out that in (1.1) is the maximum likelihood estimator of the covariance matrix .

### (c) Parameter estimators

From the sample covariance matrix in (1.1) can be defined by the random variables(3.3)From (2.1) the estimator for the azimuth *ψ* is the random variable(3.4)The estimator of the aspect ratio, *ϵ*, is, from (2.6), the random variable given by(3.5)Hence, we seek the PDFs of (3.4) and (3.5).

## 4. Statistical properties of azimuth estimator

### (a) Distribution of *Ψ*

The PDF can be derived via the unitary transformation matrixApplying it to *Z*_{0} gives where the Hermitian covariance matrix is given by with elementsThen, , while We also note that *ρ*^{2} is the squared correlation coefficient between the variables and .

Now let be *K* independent bivariate proper complex Gaussian random vectors, each with the distribution. The maximum likelihood estimator of the true covariance matrix is given by . The random variable , has PDF (Miller 1980, p. 91)(4.1)where , , *K*≥2. Here is the hypergeometric function with 2 and 1 parameters, *α*_{1}, *α*_{2} and *β*_{1}, and scalar argument *z*, which may be written explicitly as(4.2)(It is a special case of the generalized hypergeometric series, defined by (Gradshteyn & Ryzhik 1980, p. 1045)(4.3)where (*y*)_{m} is defined in terms of the Gamma function as with (*y*)_{0}=1).

On transforming (4.1) with and *ψ*=*ω*/2 we get(4.4)(4.5)The PDF in (4.4) is valid for *K*=1 as well as *K*≥2. For *K*=1 , where denotes Barakat's PDF, (Barakat 1985, eqn (5.27)), defined over −*π*/2≤*ψ*<*π*/2; this is exactly what we should obtain (Eliyahu 1993, p. 2885). Of course the importance and novelty of (4.4) is its validity for *K*≥2. Note that the PDF in (4.4) depends on *ρ*^{2} and the true value, *ψ*, of the azimuth. If we define(4.6)thenso that the degree of polarization, *P*, and the aspect ratio, *ϵ*, together determine *ρ*^{2}, and so the PDF could alternatively, but not so neatly, be directly parameterized in terms of these quantities.

The PDF is plotted in figure 3*a*–*c* for three values of *K*, *ψ* and *ρ*^{2}, respectively. We see that as *K* increases, the distribution becomes concentrated around the true value of *ψ*=0.2 in figure 3*a*, as expected. The same effect can be seen when *ρ*^{2} increases for a fixed *K* in figure 3*c*. It is also easy to see that when *ψ*=*π*/5 in figure 3*b*, the distribution shows its circular nature, clearly wrapping around in a circular manner. This is an important observation when thinking about the mean and variance of such a distribution, and we look at two specific examples to illustrate the more complicated approach required.

We consider the two covariance matrices, and where(4.7)The PDF is shown in figure 4 for both the matrices.

### (b) Mean and standard deviation for the angular distribution

Since the PDF corresponds to an *angular* distribution in the range [−*π*/4,*π*/4). As a result, special techniques are required to sensibly define the mean and variance for *Ψ*; in particular, care must be taken with the fact that the distribution is defined on a quarter circle.

We look first at the mean direction. Let . The first trigonometric moment about 0 is . Now write , where *μ*_{0} is the mean direction. The first central trigonometric moment (about the mean direction) is defined as(4.8)(4.9)so that, from (4.8)(4.10)and the mean direction is hence given by the value *μ*_{0} satisfying (4.10), agreeing with Mardia (1972, p. 70).

Using (4.2) we can write in (4.4) as(4.11)But and sowhere _{m} is given by(4.12)But _{m}=0 if *μ*_{0}=*ψ* since the overall function to be integrated is odd. Hence the solution to (4.10) is *μ*_{0}=*ψ*, and we see that the mean direction is precisely *ψ*, so that is an unbiased estimator of the exact azimuth *ψ*.

The circular variance is a measure of circular dispersion for points on the unit circle (Mardia 1972, p. 21); for *Ψ* it is defined as (Mardia 1972, pp. 45, 71)(4.13)After considerable algebra, it is found that *V*_{0} may be written as(4.14)Now (4.13) means that *V*_{0} takes values between 0 and 1, which is very different to the usual variance of a random variable which can take any positive value. To relate *V*_{0} to the standard deviation on the positive real line, the circular standard deviation is defined as (Mardia 1972, p. 74)(4.15)where the divisor of 4 takes into account that . *σ*_{0} is plotted as a function of *K* and *ρ*^{2} in figure 5; it increases with decreasing *K* and *ρ*^{2} and is maximized for *ρ*^{2}=0, which corresponds to circular polarization.

Consider again the two covariance matrices and . The columns labelled *μ*_{0} and *σ*_{0} of table 1 give the exact values of *ψ* and the circular standard deviations for these two matrices.

For the PDF appears like a periodic Gaussian (figure 4). The interval *μ*_{0}±1.96*σ*_{0} is , and numerical integration of the PDF shows that this interval covers a probability of 0.95, as in the Gaussian case.

For the interval *μ*_{0}±1.96*σ*_{0} is (−0.623, 0.345). The corresponding probability coverage is 0.93. It appears then that will be an approximate 95% CI for *ψ*=*μ*_{0}; in practice, *σ*_{0} will need to be estimated as say, so that(4.16)defines a rough 95% CI.

### (c) Simulation

For each of the model covariance matrices and , we simulated in (1.1) a number *N* of times. (For each of these repetitions, was simulated from with replaced by or , using the technique given in Medkour & Walden (2007, §V). *K*=5 was used. Each repetition allows the computation of a sample azimuth . Therefore, we obtain with . Since covers only a quarter circle we next create , *j*=1, …, *N* (Mardia 1972, p. 26). Then, let and and calculate where . The four-quadrant inverse tangent is then obtained from(4.17)(The result produced from (4.17) is the same as that obtained with the mathematical function atan2, with .) Finally, the estimated mean direction is given by . For *N*=10 000 repetitions, was obtained as −0.647 for and –0.140 for very close to the exact azimuths of –0.646 and –0.139, respectively.

The sample circular variance is given by (Mardia 1972, p. 22)and then this is converted to a measure on the positive real line as in (4.15), . Using the *N*=10 000 repetitions, was obtained as 0.078 for , and 0.244 for , compared with the theoretical values of 0.078 and 0.247, respectively, given by (4.15).

### (d) Comparison with usual mean and standard deviation

It is interesting to compare our results properly formulated for an angular distribution with what we would obtain using standard methods. We took our simulated values and computed the sample mean, , and standard deviation, , in the usual way, obtaining, for , and 0.273, respectively, very different to the values of and The standard method does not recognize that the distribution is angular, finding the mean as the centre of gravity of the PDF exactly as shown in figure 4*a*, and then inflating the standard deviation due to the probability density near to the right of the incorrect mean. By way of contrast, the properly formulated approach basically *circularizes* the PDF of figure 4*a* before determining the centre of gravity and spread.

For the sample mean and standard deviation of are given by –0.122 and 0.257 compared with the values of and . The mean is affected by the increasing tail of the PDF near *π*/4 (figure 4*b*), but the overall effect of this on the standard deviation is relatively small for this particular example.

## 5. Aspect ratio

Using the estimators in (3.3), equations (2.3) and (2.4) give the random variable(5.1)(5.2)say, where is given by (3.5). First, the PDF of *Y* in (5.1) may be derived, and then the PDF of found by the transformation in (5.2). The resulting PDF, is found in appendix A. For *K*=1,(5.3)which agrees with Barakat (1985, eqn (5.9)) as corrected in Eliyahu (1993, p. 2885). For *K*≥2,(5.4)with *α*′ and *β*′ given in (A 7). This PDF is plotted in figure 6*a*–*c* for three values of *K*, *ϵ* and *P*, respectively. We see that as *K* increases, the distribution becomes concentrated around the true value of *ϵ*=0.5 in figure 6*a*, as expected. The same effect can be seen when the degree of polarization increases for a fixed *K* in figure 6*c*.

The mean-squared error (MSE) of estimation of *ϵ* is given by and comprises the bias (squared) plus variance; it is a more useful measure of the accuracy of estimation than either the bias or the variance alone. Owing to the complexity of the PDF the *r*th moments, , are most easily computed by numerical integration. The MSE is shown in figure 7; it falls with increasing *K* or with the increasing degree of polarization *P*. We also see that the MSE increases with the absolute value of the aspect ratio, so that, as would be expected, the estimate worsens the closer the polarization approaches circular.

## 6. ULF wave example

The two-component narrowband series for the ULF wave in the solar magnetic field recorded by each of the four Cluster mission spacecraft are shown in figure 1. The measurement unit is nanoTeslas, each series has 70 values, the sample interval is 2 s and the bandpass interval 0.02–0.05 Hz. Archer *et al*. (2005) found that the wave could be assumed uniform and planar over the separations of the craft. Figure 8 shows plots of the instantaneous estimated values of (i) *ρ*^{2} defined in (4.5), estimated using , (ii) the azimuth *ψ* estimated using (3.4), (iii) the standard deviation of the azimuth estimator *σ*_{0} given in (4.15) estimated using (4.14) with *ρ*^{2} replaced by its estimate , (iv) the degree of polarization *P*, defined in (2.7), estimated using and (v) the aspect ratio *ϵ* estimated using (3.5). The degree of polarization is just a little less than unity at all times.

No assumption of temporal stationarity was necessary to produce the estimates of figure 8. However, the PDFs for the azimuth and aspect ratio estimators, which assume Gaussianity, could change with time if any of the governing parameters in (4.4) and (5.4) are time-varying.

But if a temporally stationary vector process is assumed, then the governing parameters will be time-invariant. Suppose we treat the estimated azimuth values in figure 8*b* as a random sample of azimuth estimates from the distribution in (4.4) where we take *K*=4, the now time-invariant *ψ* to be the estimated mean direction equal to 0.4, calculated as in §4*c*, and the now time-invariant *ρ*^{2} to be the median of the values in figure 8*a*. A histogram of the sample values is compared with the resulting theoretical PDF in figure 9*a*. The fit appears quite good, despite the fact that the estimated azimuth values are of course correlated in time, rather than independent as for a random sample. Our distributional results thus appear consistent with the temporal stationarity holding here.

Again, under the assumption of temporal stationarity, suppose we treat the estimated aspect ratio values in figure 8*e* as a random sample of azimuth estimates from the distribution in (5.4) where we take *K*=4, the now time-invariant *ϵ* to be the median of the values in figure 8*e*, and the now time-invariant *P* to be the median of the values in figure 8*d*. A histogram of the sample values is compared with the theoretical PDF in figure 9*b*; again, the fit is reasonable considering the necessarily crude approximations, and appears consistent with temporal stationarity.

The standard deviation of the azimuth estimator in figure 8*c* defines, via (4.16), an approximate 95% CI of the form about the instantaneous estimate. With averaging approximately 0.3, we can see that the value of *ψ* is barely constrained within its range. The difficulty of estimating the azimuth is consistent with the fact that here polarization is near circular—the aspect ratio has a median of 0.91.

## 7. Summary

We have derived the statistical distributions for instantaneous estimators of the azimuth and aspect ratio of the polarization ellipse derived from a sample covariance matrix formed from a number *K* of independent samples under the Gaussian assumption. We showed that the estimation accuracy improves as (i) *K* increases, as expected or (ii) as *ρ*^{2} increases (azimuth) or *P* increases (aspect ratio) for a fixed *K*. The distribution of the azimuth is angular and special techniques were used to derive its mean and standard deviation. An assumption of temporal stationarity of the vector process is not necessary for the validity of the theoretical results.

For the ULF wave data analysis, the empirical temporal distributions of azimuth and aspect ratio matched quite well the theoretical distributions with fixed parameters, suggesting that in this example at least temporal stationarity might hold.

## Acknowledgments

T.M. thanks the government of the People's Democratic Republic of Algeria for financial support. The authors are very grateful to Dr Tim Horbury for supplying the Cluster mission time series.

## Footnotes

- Received June 11, 2007.
- Accepted September 11, 2007.

- © 2007 The Royal Society