Spectral analysis of stationary random bivariate signals

TL;DR

This paper introduces a quaternion-based spectral analysis method for stationary bivariate signals, enabling polarization-aware spectral density estimation and decomposition into polarized and unpolarized components.

Contribution

It presents a novel quaternion Fourier transform approach for spectral analysis of bivariate signals, including a new spectral density definition and polarization periodogram.

Findings

Effective spectral density estimation demonstrated on synthetic data

Decomposition into polarized and unpolarized components achieved

The approach enhances understanding of frequency-dependent polarization attributes

Abstract

A novel approach towards the spectral analysis of stationary random bivariate signals is proposed. Using the Quaternion Fourier Transform, we introduce a quaternion-valued spectral representation of random bivariate signals seen as complex-valued sequences. This makes possible the definition of a scalar quaternion-valued spectral density for bivariate signals. This spectral density can be meaningfully interpreted in terms of frequency-dependent polarization attributes. A natural decomposition of any random bivariate signal in terms of unpolarized and polarized components is introduced. Nonparametric spectral density estimation is investigated, and we introduce the polarization periodogram of a random bivariate signal. Numerical experiments support our theoretical analysis, illustrating the relevance of the approach on synthetic data.

Equations

$$q=a+b\mathbf{i}+c\mathbf{j}+d\mathbf{k},$$

$${\mathbf{i}}^2={\mathbf{j}}^2={\mathbf{k}}^2=\mathbf{i}\mathbf{j}\mathbf{k}=-1.$$

$$q=|q|\exp [\mathbf{i}\theta ]\exp [-\mathbf{k}\chi ]\exp [\mathbf{j}\phi ],$$

$$X(\nu ):=\sum _{t=-\infty }^{+\infty }x[t]\exp (-\mathbf{j}2\pi \nu t).$$

$$x[t]={\int }_{-1/2}^{+1/2}X(\nu )\exp (\mathbf{j}2\pi \nu t)\mathrm{d}\nu .$$

$$X(\nu )=U(\nu )+\mathbf{i}V(\nu ),U(\nu ),V(\nu )\in {\mathbb{C}}_{\mathbf{j}},$$

$$X(-\nu )={\overline{X(\nu )}}^{\mathbf{i}}.$$

$$\mathbf{E}\left\{x[t]\right\}$$

$$R_{uu}[t,\tau ]$$

$$R_{vv}[t,\tau ]$$

$$R_{uv}[t,\tau ]$$

$$x[t]={\int }_{-1/2}^{+1/2}\mathrm{d}X(\nu )\exp (\mathbf{j}2\pi \nu t),$$

$$\mathrm{d}X(\nu )=\mathrm{d}U(\nu )+\mathbf{i}\mathrm{d}V(\nu ).$$

$$x[t]={\int }_{-1/2}^{+1/2}\mathrm{d}X(\nu )\exp (\mathbf{j}2\pi \nu t),$$

$$\mathbf{E}\left\{\mathrm{d}X(\nu )\overline{\mathrm{d}X({\nu }^{\prime })}\right\}=\mathbf{E}\left\{\mathrm{d}X(\nu )\mathrm{d}X({\nu }^{\prime }{)}^{\star \mathbf{j}}\right\}=0,$$

$${\int }_{-1/2}^{+1/2}{\Gamma }_{xx}(\nu )\mathrm{d}\nu =\mathbf{E}\left\{|x[t]{|}^2\right\}+\mathbf{E}\left\{{\left|x[t]\right|}_{\mathbf{j}}^2\right\}\mathbf{j}\in \mathbb{H}$$

$$\mathrm{d}X(-\nu )={\overline{\mathrm{d}X(\nu )}}^{\mathbf{i}}.$$

$${\Gamma }_{xx}(-\nu )={\Gamma }_{xx}(\nu {)}^{\star \mathbf{i}}.$$

$$\mathbf{E}\left\{\mathrm{d}X(\nu )\mathrm{d}X(\nu {)}^{\star \mathbf{i}}\right\}=\mathbf{E}\left\{\mathrm{d}X(\nu )\mathrm{d}X(\nu {)}^{\star \mathbf{k}}\right\}=0.$$

$${\gamma }_{xx}[\tau ]={\int }_{-1/2}^{+1/2}{\Gamma }_{xx}(\nu )\exp (\mathbf{j}2\pi \nu \tau )\mathrm{d}\nu ,$$

$$\begin{array}{rl}{\gamma }_{xx}[\tau ]& =R_{uu}[\tau ]+R_{vv}[\tau ]\\ & +(R_{uu}[\tau ]-R_{vv}[\tau ])\mathbf{j}+2R_{uv}[\tau ]\mathbf{k}.\end{array}$$

$${\Gamma }_{xy}(\nu )\mathrm{d}\nu =\mathbf{E}\left\{\mathrm{d}X(\nu )\overline{\mathrm{d}Y(\nu )}\right\}+\mathbf{E}\left\{\mathrm{d}X(\nu )\mathrm{d}Y(\nu {)}^{\star \mathbf{j}}\right\}\mathbf{j}$$

$${\gamma }_{xy}[\tau ]=R_{u_x{u}_y}[\tau ]+R_{v_y{v}_x}[\tau ]+(R_{u_y{v}_x}[\tau ]-R_{u_x{v}_y}[\tau ])\mathbf{i}$$

$$+(R_{u_x{u}_y}[\tau ]-R_{v_y{v}_x}[\tau ])\mathbf{j}+(R_{u_y{v}_x}[\tau ]+R_{u_x{v}_y}[\tau ])\mathbf{k}.$$

$$\mathbf{E}\left\{\mathrm{d}X(\nu )\overline{\mathrm{d}Y(\nu )}\right\}=\sum _{\tau =-\infty }^{+\infty }\mathbf{E}\left\{x[t]e^{-\mathbf{j}2\pi \nu \tau }\overline{y[t-\tau ]}\right\}$$

$$\mathbf{E}\left\{\mathrm{d}X(\nu )\mathrm{d}Y(\nu {)}^{\star \mathbf{j}}\right\}=\sum _{\tau =-\infty }^{+\infty }\mathbf{E}\left\{x[t]e^{-\mathbf{j}2\pi \nu \tau }y[t-\tau {]}^{\star \mathbf{j}}\right\}$$

$$\mathcal{S}({\Gamma }_{xx}(\nu ))$$

$$\mathcal{V}({\Gamma }_{xx}(\nu ))$$

$${\gamma }_{x+y,x+y}[\tau ]={\gamma }_{xx}[\tau ]+{\gamma }_{yy}[\tau ]$$

$$S_0(\nu )$$

$$S_1(\nu )$$

$$S_2(\nu )$$

$$S_3(\nu )$$

$${\Gamma }_{xx}(\nu )=S_0(\nu )+\mathbf{i}{S}_3(\nu )+\mathbf{j}{S}_1(\nu )+\mathbf{k}{S}_2(\nu ).$$

$$\begin{array}{rl}\mathbf{E}\left\{|\mathrm{d}X(\nu ){|}^2\right\}& =\mathbf{E}\left\{|\mathrm{d}U(\nu ){|}^2\right\}+\mathbf{E}\left\{|\mathrm{d}V(\nu ){|}^2\right\}\\ & =(P_{uu}(\nu )+P_{vv}(\nu ))\mathrm{d}\nu ,\end{array}$$

Full text

I Introduction

Random bivariate signals are 2D vector timeseries. They appear in a large variety of applications, ranging from oceanography [1, 2], to optics [3], radar [4], geophysics [5] or EEG analysis [6] to name but a few. A bivariate signal is usually decomposed in two orthogonal components $u[t]$ and $v[t]$. Thus a bivariate signal $x[t]$ can be either represented as the vector signal $x[t] =(u[t],\>v[t])^{T}\in\mathbb{R}^{2}$ or the complex-valued signal $x[t]=u[t]+\bm{i}v[t]$.

The statistical analysis of signals with vector-valued samples can be carried out using standard multivariate time series analysis techniques (see e.g [7, chap. 9] or [8, chap. 11]), bivariate signals are no exception. However, in the signal processing community, bivariate signals have often been described using complex-valued models [9, 10, 11, 12, 13]. To account for the full second-order statistical characterization of the complex signal $x[t]=u[t]+\bm{i}v[t]$, a usual approach is to define two quantities: the usual autocovariance function and the complementary-covariance function (the relation function in [9]). This leads to the definition of the augmented vector $(x[t],\overline{x[t]})^{T}\in\mathbb{C}^{2}$ from the signal and its conjugate, and to the related augmented covariance and spectral density matrices [12, chap. 8].

The rotary spectrum analysis [1, 2] is a well-known technique rooted in oceanographic studies. It is based on the decomposition of a complex-valued signal into clockwise and counterclockwise rotating components. This seminal approach has stimulated many theoretical developments [14, 15, 16, 17, 18, 19]. The rotary components method is also related to polarization analysis [14, 12]. The focus on polarization or rotary components usually depends on the field of application: rotary components are more common in oceanography, while optics and radar scientists usually deal with polarization [20, 21, 22, 23, 24].

Existing approaches all rely on the same spectral representation of stationary bivariate processes, which is based on the standard Fourier Transform (FT). For complex signals, it means that negative frequencies must be taken into account, as they provide information about the process. Therefore one needs to consider spectral matrices rather than a scalar spectral density. As a consequence, meaningful physical parameters are not directly “readable” in the state-of-the-art formulations.

We propose a new approach to analyze the spectral content of stationary random bivariate signals. It is based on recent results from [25, 26] and extended to the case of stationary bivariate signals seen as complex-valued signals. This paper provides a well-suited framework for the analysis of stationary bivariate signals which naturally describes the spectral content and the “geometric” (also sometimes called polarization) content of bivariate signals. Thanks to the definition of the dedicated Quaternion Fourier Transform (QFT), it is possible to describe the spectral content of such signals in terms of polarized and unpolarized parts, which both encode meaningful information about the signal.

This paper structure is as follows. Section II reviews the necessary material regarding quaternions and the QFT. Section III gives the central results of this paper: we introduce the scalar quaternion-valued spectral density of a bivariate signal and its subsequent properties. Results are compared with state-of-the art approaches. In particular the differences between second-order circularity, also called properness, and polarization are stressed. Simple explicit examples are finally presented. Section IV deals with nonparametric spectral density estimation, and introduces the polarization periodogram. Our theoretical analysis is supported by numerical experiments in Section V. Section VI gathers concluding remarks.

II Quaternion Fourier Transform

II-A Quaternion algebra

We review the basic material regarding quaternions and refer to more detailed textbooks (e.g. [27]) for a complete overview. Quaternions form a four dimensional noncommutative algebra. Any quaternion $q\in\mathbb{H}$ can be written in its Cartesian form as

[TABLE]

where $a,b,c,d\in\mathbb{R}$ and $\bm{i},\bm{j},\bm{k}$ are roots of $-1$ satisfying

[TABLE]

The canonical elements $\bm{i},\bm{j},\bm{k}$, together with the identity of $\mathbb{H}$ form the quaternion canonical basis given by $\left\{1,\bm{i},\bm{j},\bm{k}\right\}$. We will use the notation $\mathcal{S}(q)=a\in\mathbb{R}$ to define the scalar part of the quaternion $q$, and $\mathcal{V}(q)=q-\mathcal{S}(q)\in\mathrm{span}\left\{\bm{i},\bm{j},\bm{k}\right\}$ to denote its vector part. We can define the real and imaginary parts of a quaternion $q$ as $\mathfrak{R}(q)=a,\>\mathfrak{I}_{\bm{i}}(q)=b,\>\mathfrak{I}_{\bm{j}}(q)=c,\>\mathfrak{I}_{\bm{k}}(q)=d.$ A quaternion is called pure if its real (or scalar) part is equal to zero, that is $a=0$, e.g. $\bm{i},\bm{j},\bm{k}$ are pure quaternions. The quaternion conjugate of $q$ is $\overline{q}=\mathcal{S}(q)-\mathcal{V}(q)$. The modulus of a quaternion $q\in\mathbb{H}$ is defined by $|q|^{2}=q\overline{q}=\overline{q}q=a^{2}+b^{2}+c^{2}+d^{2}$. The inverse of a non-zero quaternion is defined by $q^{-1}=\overline{q}/|q|^{2}$. Importantly quaternion multiplication is noncommutative, that is in general for $p,q\in\mathbb{H}$, one has $pq\neq qp$. Involutions with respect to $\bm{i},\bm{j},\bm{k}$ are defined as $\overline{q}^{\bm{i}}=-\bm{i}q\bm{i},\>\overline{q}^{\bm{j}}=-\bm{j}q\bm{j},\>\overline{q}^{\bm{k}}=-\bm{k}q\bm{k}$. The combination of conjugation and involution with respect to an arbitrary pure quaternion $\bm{\mu}$ is denoted by $q^{\star\bm{\mu}}:=\overline{\left(\overline{q}\right)}^{\bm{\mu}}=\overline{\left(\overline{q}^{\bm{\mu}}\right)}$ and for instance $(a+b\bm{i}+c\bm{j}+d\bm{k})^{\star\bm{j}}=a+b\bm{i}-c\bm{j}+d\bm{k}$. For later use, we also introduce the notation $\left|q\right|_{\bm{j}}^{2}:=qq^{\star\bm{j}}$.

Quaternions encompass complex numbers. One can construct complex subfields of $\mathbb{H}$, e.g… $\mathbb{C}_{\bm{j}}=\mathrm{span}\left\{1, \bm{j}\right\}$ or $\mathbb{C}_{\bm{i}}=\mathrm{span}\left\{1,\bm{i}\right\}$ which are isomorphic to $\mathbb{C}$. Any quaternion can be seen as a pair of complex numbers: let us mention the symplectic decomposition $q=q_{1}+\bm{i}q_{2},\>q_{1},q_{2}\in\mathbb{C}_{\bm{j}}$, where the quaternion $q$ is splitted into two $\mathbb{C}_{\bm{j}}$-valued complex numbers. This form is particularly suited for computations performed later on with the quaternion Fourier transform.

Polar forms of quaternions exist. For an arbitrary pure unit quaternion $\bm{\mu}$ and $\theta\in\mathbb{R}$, we have $\exp(\bm{\mu}\theta)=\cos\theta+\bm{\mu}\sin\theta$. It generalizes the notion of complex exponentials and the following polar form was proposed in [28]:

[TABLE]

with $(\theta,\chi,\varphi)\in[-\pi/2,\pi/2]\times[-\pi/4,\pi/4]\times[-\pi,\pi]$. This form is particularly useful for quaternion embedding of complex signals, see [25], and in the spectral description of stationary signals in Section III-G.

II-B Quaternion Fourier Transform

We review here briefly the Quaternion Fourier Transform (QFT) introduced in [29] and studied in detail recently in [25]. We refer the reader to these articles for proofs and a detailed presentation.

Here we consider only discrete-time (DT) signals: $t$ is a time index such that $x(t\Delta)=x[t]$, where $\Delta$ is the sampling step. We assume $\Delta=1$ in the rest of this paper.

We define the Discrete Time Quaternion Fourier Transform (hereafter denoted QFT) of axis $\bm{j}$ of a signal $x:\mathbb{R}\rightarrow\mathbb{H}$ by

[TABLE]

The inverse QFT is given by

[TABLE]

The above relations are directly obtained by discretizing the continuous-time QFT presented in [25, 29] using similar arguments as with the usual FT.

Write $x[t]=u[t]+\bm{i}v[t]$, $u[t],v[t]\in\mathbb{C}_{\bm{j}}$, its QFT reads

[TABLE]

where $U(\nu),V(\nu)$ are the standard FTs of $u,v$: the QFT is performing two standard FT. This may explain why this QFT shares most properties of the classical FT, see [25].

Another benefit of the QFT is that it provides a well-suited framework for bivariate signals. A bivariate signal can be written as a $\mathbb{C}_{\bm{i}}$-valued signal $x[t]=u[t]+\bm{i}v[t]$, with $u,v$ real signals. The QFT of such a signal exhibits an $\bm{i}$-hermitian symmetry [29]:

[TABLE]

Eq. (7) shows that, when using the QFT with $x[t]\in\mathbb{C}_{\bm{i}}$, negative frequencies carry no additional information to positive frequencies about the signal.

The $\bm{i}$-hermitian symmetry (7) permits the construction of the quaternion embedding of a complex signal, by canceling out negative frequencies of the spectrum. The quaternion embedding of a complex signal is a direct bivariate counterpart of the usual analytic signal and permits to identify both instantaneous phase and polarization (*i.e. *geometric) properties of a complex signal.

III Spectral representation of bivariate stationary random processes

Any bivariate discrete-time random process $x[t]$ can be decomposed as $x[t]=u[t] +\bm{i}v[t]$, where $u[t],v[t]$ are real-valued discrete-time random processes. The process $x[t]$ is said to be second-order stationary if $u[t]$ and $v[t]$ are jointly second-order stationary, that is [7, p. 655]:

[TABLE]

Here $\mathbf{E}\left\{\cdot\right\}$ denotes the mathematical expectation and $R_{uu},R_{vv}$ and $R_{uv}$ denote usual autocovariance sequences (ACVS) and crosscovariance sequences (CCVS) between real-valued sequences. Second-order stationarity ensures that first and second-order moments are finite and do not depend on $t$.

Let us make some further assumptions. First all processes are assumed to be zero-mean, *i.e. *$m=0$. Also, we suppose that any of the ACVS or CCVS are absolutely summable, so that their usual Fourier transform exist. Second-order stationarity is simply referred to as stationarity in the sequel.

III-A Main result

Using the QFT, we derive a spectral representation theorem for any bivariate stationary process $x[t]=u[t]+\bm{i}v[t]$. The existence of the spectral increments $\mathrm{d}X(\nu)$ follows from the existence of the (usual) spectral increments of $u[t]$ and $v[t]$, see for instance [30, p. 36], [7, p. 246] or [31, p. 344]. Namely, we can write $x[t]$ as the (quaternion) Fourier-Stieltjes integral

[TABLE]

where the spectral increments $\mathrm{d}X(\nu)=X(\nu+\mathrm{d}\nu)-X(\nu)$ are quaternion-valued and $X(\nu)$ is an independent additive random measure. The existence roots in two main properties. The QFT restricted to $\mathbb{C}_{\bm{j}}$-valued processes is isomorphic to the standard FT so that the spectral increments $\mathrm{d}U(\nu),\mathrm{d}V(\nu)$ of $u[t],v[t]$ are $\mathbb{C}_{\bm{j}}$-valued. Moreover the QFT is left-quaternion-linear, that is $\forall\lambda\in\mathbb{H}$, the QFT of $\lambda x[t]$ is $\lambda X(\nu)$ so that

[TABLE]

Theorem 1** (Spectral representation of bivariate stationary random processes).**

Let $x[t]=u[t]+\bm{i}v[t]$ be a bivariate stationary process. Suppose that $u[t]$ and $v[t]$ are both harmonizable. Then there exists a quaternion-valued orthogonal process $X(\nu)$ such that

[TABLE]

the integral being defined in the mean-square sense. The process $X(\nu)$ has the following properties:

$\forall\nu,\>\mathbf{E}\left\{\mathrm{d}X(\nu)\right\}=0$, 2. 2.

$\forall\nu,\>\mathbf{E}\left\{|\mathrm{d}X(\nu)|^{2}\right\}+\mathbf{E}\left\{\left|\mathrm{d}X(\nu)\right|_{\bm{j}}^{2}\right\}\bm{j}=\Gamma_{xx}(\nu)\mathrm{d}\nu$, where $\Gamma_{xx}(\nu)$ is the spectral density of $x$, 3. 3.

For any $\nu\neq\nu^{\prime}$, we have

[TABLE]

which shows that the spectral increments $dX(\nu)$ are two times orthogonal.

Appendix A proves this theorem that follows the derivation given by Priestley [7], adapted to the QFT setting. Properties 1), 2) and 3) essentially come from the self and joint properties of the spectral increments $\mathrm{d}U(\nu)$ and $\mathrm{d}V(\nu)$. Property 2) introduces the quaternion-valued spectral density of $x$. This is in fact a power spectral density since (see Appendix A) one has

[TABLE]

where the right-hand side contains all the power information of the process $x[t]$.

Note that four (real) power-related quantities are necessary to describe the second-order properties of a bivariate signal, see Section III-C. This argues in favor of the definition of relevant scalar quaternion-valued quantities such as $\Gamma_{xx}(\nu)$.

Note also that in our definition $\Gamma_{xx}(\nu)$ is a generalized function: we can use Dirac distributions. Thus $\Gamma_{xx}(\nu)$ does not need to be continuous everywhere.

Since $x[t]=u[t]+\bm{i}v[t]\in\mathbb{C}_{\bm{i}}$, the spectral increments additionally satisfy the same $\bm{i}$-hermitian symmetry as the QFT of deterministic $\mathbb{C}_{\bm{i}}$-valued signals, i.e.

[TABLE]

As a result, the spectral density $\Gamma_{xx}(\nu)$ has symmetry

[TABLE]

This result shows again that the study of $\mathbb{C}_{\bm{i}}$-valued (bivariate) signals can be performed using only positive frequencies of its quaternion-valued spectral representation. At each (positive) frequency, a quaternion-valued quantity summarizes both power and polarization properties of $x$. This will be detailed in Sections III-C and III-D.

The spectral increments $\mathrm{d}X(\nu)$ are quaternion-valued random variables (RV). It is usual to describe the full second-order statistical structure of a quaternion RV $q$ by the four covariances $\mathbf{E}\left\{qq^{\star\bm{\mu}}\right\}$, $\bm{\mu}=\bm{i},\bm{j},\bm{k}$. These covariances often obey some symmetries characterized by the notion of properness. Properness levels of quaternion RV have been investigated by several authors [32, 33, 34] and reviewed recently in [35]. The spectral increments of a $\mathbb{C}_{\bm{i}}$-valued process $x[t]$ satisfy the symmetry (16) and thus property 3) of theorem 1 with $\nu^{\prime}=-\nu$ yields

[TABLE]

Eq. (18) shows that the spectral increments $\mathrm{d}X(\nu)$ are $(1,\bm{j})$-proper in the classification of [35], also denoted as $\mathbb{C}_{\bm{j}}$-properness in [33].

III-B Covariances, Wiener-Khintchine theorem

Thanks to theorem 1, we are able to describe the spectral content of random bivariate signals by their spectral density. Usually, this spectral density is introduced by using Wiener-Khintchine theorem once the autocorrelation of the process has been defined. This is not the case here due to the non-commutativity of $\mathbb{H}$: the notion of (auto-)covariance must be carefully defined if one wants to recover a Wiener-Khintchine theorem for quaternion valued processes. To this aim, a natural approach can be to define the autocovariance of $x$ by inverse QFT of $\Gamma_{xx}(\nu)$:

[TABLE]

so that $\gamma_{xx}[\tau]$ is explicitly given by (see Appendix A)

[TABLE]

The autocovariance function $\gamma_{xx}[\tau]$ takes its values in $\mathrm{span}\left\{1,\bm{j},\bm{k}\right\}$. It is not symmetric in $\tau$, as the term $R_{uv}[\tau]$ is not symmetric in general. More generally one can define the cross-spectral density between two bivariate stationary random processes $x[t]=u_{x}[t]+\bm{i}v_{x}[t]$ and $y[t]=u_{y}[t]+\bm{i}v_{y}[t]$, where $u_{k}[t],v_{k}[t]\in\mathbb{R}$, $k=x,y$. Let us denote by $\mathrm{d}X(\nu)$ and $\mathrm{d}Y(\nu)$ their spectral increments. The cross-spectral density between $x$ and $y$ is

[TABLE]

so that their cross-covariance defined as its inverse QFT reads:

[TABLE]

This quaternion-valued cross-covariance encodes the full statistical information about $x$ and $y$. Eq. (20) and (III-B) may sound disappointing at first glance, but there is no simple expression of those equations in terms of usual covariance $\mathbf{E}\left\{x[t]\overline{y[t-\tau]}\right\}$ and complementary covariance $\mathbf{E}\left\{x[t]y[t-\tau]\right\}$. However the following Wiener-Khintchine like theorem directly connects $x[t]$ and $y[t]$ to the cross-spectral density (21).

Theorem 2**.**

Let $x$ and $y$ be two stationary random bivariate signals. Then

[TABLE]

In the special case $y[t]=x[t]$ one has

[TABLE]

Sketch of proof.

We start by developing both sides using expressions of $\mathrm{d}X(\nu)$ et $\mathrm{d}Y(\nu)$ in terms of $\mathrm{d}U_{k}(\nu),\mathrm{d}V_{k}(\nu)$, $k=x,y$ and $x[t],y[t]$ in terms of $u_{k}[t],v_{k}[t]$, $k=x,y$. Then usual rules of quaternion calculus (e.g. $\overline{q_{1}q_{2}}=\overline{q_{2}}\>\overline{q_{1}}$ for $q_{1},q_{2}\in\mathbb{H}$; $\bm{i}q=\overline{q}\bm{i}$ if $q\in\mathbb{C}_{\bm{j}}$) permit to simplify both sides. Standard Wiener-Khintchine theorems for real signals lead to the result. ∎

Let us note finally the following property.

Proposition 1** (Autocorrelation of a sum of independent signals).**

If $x$ and $w$ are independent, $\mathbb{C}_{\bm{i}}$-valued, stationary processes then

[TABLE]

Proof.

By direct calculation. ∎

Proposition 1 is a desirable result, which permits to manipulate quaternion valued autocovariance functions like standard autocovariance functions. Note that this result applies with spectral densities of independent signals $x$ and $y$ as well: the spectral density of $x+y$ is the sum of their spectral densities.

III-C Spectral density and Stokes parameters

The spectral density $\Gamma_{xx}(\nu)$ is directly related to Stokes parameters, which are fundamental quantities used in polarization of electromagnetic waves [36, 37]. Stokes parameters are given by [12, 37]

[TABLE]

where we have introduced the usual spectral densities of $u$ and $v$, $P_{uu}$ and $P_{vv}$, as well as the usual cross-spectral density $P_{uv}$.

Theorem 3**.**

Let $\Gamma_{xx}(\nu)$ be defined by theorem 1. It can be re-expressed like

[TABLE]

where $S_{\alpha}(\nu)$, $\alpha=0,1,2,3$ are the Stokes parameters of $x$.

Proof.

The two terms appearing in theorem 1 can be expressed in terms of the spectral increments $\mathrm{d}U(\nu)$ and $\mathrm{d}V(\nu)$:

[TABLE]

and

[TABLE]

Then (32) follows from (28) – (31). ∎

Eq. (32) has a powerful geometric interpretation. Stokes parameters permit to separate the contribution of polarized and unpolarized components, as discussed next in Section III-E. The scalar part $S_{0}(\nu)$ of $\Gamma_{xx}(\nu)$ is the total power at frequency $\nu$, *i.e. *the sum of the power of the polarized and unpolarized parts. The vector part describes only the polarized part of $x$ at frequency $\nu$.

III-D Poincaré sphere and degree of polarization

Fig. 1 depicts the Poincaré sphere of polarization states [20, p. 125][38, 37]. At frequency $\nu$, the vector part of $\Gamma_{xx}(\nu)$ – normalized by its scalar part $S_{0}(\nu)$ – identifies a point on this Poincaré sphere. The angular coordinates $(2\theta,2\chi)$ are directly related to the mean ellipse properties of the signal, *i.e. *$\theta$ is the mean orientation and $\chi$ is the mean ellipticity.

At each frequency, the radius of the Poincaré sphere is called the degree of polarization $\Phi(\nu)$. Namely,

[TABLE]

where $\mathcal{V}(\cdot)$ and $\mathcal{S}(\cdot)$ denote the vector and scalar part, respectively. It follows from the definition of $\Phi(\nu)$ that $0\leq\Phi(\nu)\leq 1$ for all $\nu$. The degree of polarization $\Phi(\nu)$ quantifies the repartition between polarized and unpolarized components. This motivates the following vocabulary: the process $x$ is said to be

• •

fully polarized at frequency $\nu$ if $\Phi(\nu)=1$,

• •

unpolarized at frequency $\nu$ if $\Phi(\nu)=0$,

• •

partially polarized at frequency $\nu$ if $0<\Phi(\nu)<1$.

The degree of polarization is a quantity of fundamental interest in many fields (see e.g. [39, 40]).It is invariant by change of reference frame, making it a robust parameter of interest.

III-E Decomposition into polarized and unpolarized parts

The Unpolarized/Polarized part decomposition (hereafter termed UP decomposition) has been evoked in Section III-C. Now we formally show that any bivariate stationary signal can be uniquely decomposed, as the sum of two uncorrelated processes: a completely polarized process and an unpolarized process. Let us rewrite the spectral density of theorem 1 as

[TABLE]

where the ${\tt u}$ and ${\tt p}$ superscripts stand for unpolarized and polarized part, respectively. The decomposition (36) is unique. Using Stokes parameters, we get

[TABLE]

which is the quaternion counterpart of the decomposition given in standard optics textbooks, see e.g. [20, p. 127], [37, p. 551]. Eq. (37) highlights how the degree of polarization rules the power repartition between the polarized and unpolarized parts of the spectral density.

Back to the time-domain, there exists two uncorrelated processes $x_{\tt p}$ and $x_{\tt u}$ such that $x[t]=x_{\tt p}[t]+x_{\tt u}[t]$ and

[TABLE]

with $\mathrm{d}X_{\tt p}$ and $\mathrm{d}X_{\tt u}$ the spectral increments corresponding to the polarized and unpolarized part, respectively.

III-F Comparison with previous work

III-F1 Proper and improper signals

The notion of (im)properness of a complex signal has attracted much interest in the signal processing community over the last two decades, see [9, 13, 12] and references therein. To account for the full second-order statistical structure of a stationary complex signal, one has to consider both the usual autocovariance $R_{xx}[\tau]$ and the complementary covariance $\tilde{R}_{xx}[\tau]$ such that:

[TABLE]

Proper signals are characterized by a zero complementary covariance sequence, meaning that a signal $x[t]$ is not correlated with its complex conjugate $\overline{x[t-\tau]}$, for all $\tau$. It follows that

[TABLE]

A direct consequence is that the spectral density (32) of a proper signal $x[t]$ reads

[TABLE]

as conditions (41) are equivalent to $S_{1}(\nu)=S_{2}(\nu)=0$ for all $\nu$. Eq. (42) shows that a proper signal is in general partially circularly polarized. This highlights that polarization and properness of complex random signals are distinct concepts and therefore shall not be confused.

III-F2 Relation to the rotary spectrum approach

To illustrate the relevance of the quaternion-valued spectral density $\Gamma_{xx}(\nu)$ defined in theorem 1, we compare it to the well-known rotary spectrum approach [1, 2, 15]. This later method decomposes a bivariate signal into a sum of clockwise (CW) and counterclockwise phasors (CCW) – the so-called rotary components [15]. The determination of the rotary components relies on the usual spectral density $P_{xx}(\nu)$ and complementary spectral density $\tilde{P}_{xx}(\nu)$ defined by standard FTs of $R_{xx}[\tau]$ and $\tilde{R}_{xx}[\tau]$, respectively [12].

For $\nu>0$, the CW rotary power spectrum is given by $P_{xx}(\nu)$, while $P_{xx}(-\nu)$ gives the CCW rotary power spectrum. The rotary coherence (*i.e. *correlation between CW and CCW components) is controled by $\tilde{P}_{xx}(\nu)$, which is in general complex-valued.

The rotary spectra can be expressed in terms of Stokes parameters like [12, p. 213]

[TABLE]

Since $S_{0}(\nu)$ is even and $S_{3}(\nu)$ is odd, $P_{xx}(\nu)$ shows no particular symmetry. Moreover, we see that $P_{xx}(\nu)$ combines in one real scalar two very different quantities: $S_{0}(\nu)$ is related to the total power, and $S_{3}(\nu)$ gives the (signed) power of the circularly polarized part. This is not surprising since the pair $(P_{xx}(\nu),\tilde{P}_{xx}(\nu))$ was introduced to account for improperness properties of complex-valued signals, not polarization properties. In contrast the use of $\Gamma_{xx}(\nu)$ provides directly the total power information decoupled from the polarization information.

The rotary spectrum and the quaternion-valued approach provide equivalent representations. However, the quaternion-valued spectral density $\Gamma_{xx}$ provides a direct interpretation of physical quantities, the Stokes parameters. These parameters appear naturally in the components of $\Gamma_{xx}$. The use of the QFT to study bivariate signals provides a well-suited framework for a meaningful and rooted in physics “geometric spectral analysis”. This approach demonstrates that the use of higher dimensional algebra in the definition of the FT permits to avoid ad hoc constructions while revealing intrinsic relevant parameters of bivariate signals.

III-G Examples

III-G1 Deterministic signals

In this case the spectral density reads

[TABLE]

where expectations have been dropped out. Moreover one gets

[TABLE]

so that a deterministic signal is always fully polarized at all frequencies for which $X(\nu)\neq 0$. The other polarization parameters $a,\theta$, $\chi$ or $S_{1},S_{2},S_{3}$ can be determined as well, depending on the signal. Let the bivariate monochromatic signal $x$ such that

[TABLE]

with $a,\chi,\theta$ the parameters of the elliptic polarization. The QFT of $x$ reads111For the sake of simplicity, we restrict our analysis to the frequency support $[-1/2,1/2]$, although $X(\nu)$ is periodic. Note also the use of the Euler polar form (3) with a null $\bm{j}$-phase.

[TABLE]

which gives the following spectral density $\Gamma_{xx}(\nu)$:

[TABLE]

where $\Gamma_{xx}^{\nu_{0}}=S_{0}+\bm{i}S_{3}+\bm{j}S_{1}+\bm{k}S_{2}$. The autocovariance sequence is then obtained by inverse QFT:

[TABLE]

The autocovariance is not symmetric. The value of $S_{3}$ controls the odd contribution, whereas the remaining terms are all even. It is interesting to note that the autocovariance function of a monochromatic signal is even if and only if the signal is linearly polarized.

III-G2 Bivariate white noise

Consider the process $w[t]=u[t]+ \bm{i}v[t]$ where $u,v$ are both real, i.i.d. and jointly second-order stationary with properties:

[TABLE]

Since $u,v$ are i.i.d., their autocovariances are

[TABLE]

which yields the autocorrelation of $x$ using (20)

[TABLE]

The spectral density is obtained by QFT:

[TABLE]

This spectral density is constant. It has no $\bm{i}$-component, so that $S_{3}(\nu)=0$ for all $\nu$. As a consequence, an i.i.d. arbitrary second-order stationary bivariate or complex white noise shows no ellipticity. It is always either unpolarized, or linearly polarized (fully or partially). The polarization properties are identical at all frequencies.

The polarization degree defined by (35) is:

[TABLE]

where we see that $x[t]$ is unpolarized at all frequencies iff $\sigma_{u}=\sigma_{v}$ and $\rho_{uv}=0$. When $\Phi\neq 0$, the angle $\theta$ of the linear polarization is given by $\theta=0$ if $\rho_{uv}=0$ and by

[TABLE]

when $\rho_{uv}\neq 0$ and where $\mathrm{atan2}$ denotes the four-quadrant inverse tangent.

The UP decomposition of bivariate white noise gives a simple procedure to simulate bivariate white noise with desired polarization properties. Let $0\leq\Phi\leq 1$ be the desired degree of polarization, $\theta\in[-\pi/2,\pi/2]$ the orientation angle and $S_{0}>0$ the total intensity. Let $w^{\texttt{u}}[t]$ be an unpolarized white noise, $\mathbb{C}_{\bm{i}}$-valued, such that $R_{w^{\texttt{u}}w^{\texttt{u}}}[\tau]=\delta_{\tau,0}$. Let $w^{\texttt{p}}[t]$ be a real-valued white noise sequence of unit variance. Assume further that $w^{\texttt{u}}[t]$ and $w^{\texttt{p}}[t]$ are independent. Then the white noise $w[t]$ constructed as

[TABLE]

has exactly the desired polarization properties.

III-G3 Bivariate monochromatic signal in white noise

Consider the signal $y[t]=x[t]+w[t]$, where $x[t]$ is a bivariate monochromatic signal and $w[t]$ is a bivariate white noise (i.i.d. sequence). Assume moreover that $x[t]$ and $w[t]$ are independent. From Proposition 1 the spectral density of $y$ is directly given by the sum of their spectral densities:

[TABLE]

where $\Gamma_{ww}(\nu)$ is given by (53) and $\Gamma_{xx}^{\nu_{0}}=S_{0,x}+\bm{i}S_{3,x}+\bm{j}S_{1,x}+\bm{k}S_{2,x}$; the $S_{\alpha,x}$ denote the Stokes parameters of $x$. At $\nu=\nu_{0}$ the spectral density writes $\Gamma_{yy}^{\nu_{0}}=S_{0,y}+\bm{i}S_{3,y}+\bm{j}S_{1,y}+\bm{k}S_{2,y}$ where

[TABLE]

First, we see that $S_{1}$ and $S_{2}$ parameters are mixing polarization properties of $x$ and $w$. Since $S_{3}$ is not modified in presence of i.i.d. white noise, only the direction of polarization changes, not the ellipticity. The output degree of polarization (at $\nu=\nu_{0}$) takes a simple form when the noise is unpolarized:

[TABLE]

where $\mathrm{SNR}=S_{0,x}/(\sigma^{2}_{u}+\sigma_{v}^{2})$ is the signal-to-noise ratio (SNR). The degree of polarization decreases with the SNR.

IV Spectral density estimation

We propose two nonparametric spectral density estimation methods, and derive their properties. In particular, we investigate the problem of estimating the degree of polarization. In the remainder of this paper, we consider a bivariate stationary signal $x[t]=u[t]+\bm{i}v[t]$ consisting in $N$ samples such that $t=0,1,\ldots N-1$ and with sampling size $\Delta=1$.

IV-A A naive spectral estimator: the polarization periodogram

The first basic spectral density estimator is the polarization periodogram $\hat{\Gamma}_{xx}^{(p)}(\nu)$. The underlying rationale is very close to the derivation of the usual periodogram. One starts by computing an estimator $\hat{\gamma}_{xx}^{(p)}[\tau]$ of the autocovariance sequence $\gamma_{xx}[\tau]$. It is done by combining usual (biased) estimators of auto- and cross-covariance sequences:

[TABLE]

where $\hat{R}_{uv}^{(p)}[\tau]=\hat{R}_{vu}^{(p)}[-\tau]$ for $\tau=-1,\ldots,-(N-1)$ and $\hat{R}_{uv}^{(p)}[\tau]=0$ for $|\tau|\geq N$. Auto-covariance estimators follow from (60). Then taking the QFT of $\hat{\gamma}_{xx}^{(p)}[\tau]$ given by (20) yields the polarization periodogram, which reads (after simplification)

[TABLE]

Alike the classical periodogram, this estimator is a biased, inconsistent estimator of the spectral density $\Gamma_{xx}(\nu)$. One has

[TABLE]

where $\mathcal{F}_{N}(\nu)=\sin^{2}(\pi N\nu)/[N\sin^{2}(\pi\nu)]$ is the Fejér kernel. It follows that the polarization periodogram is only asymptotically unbiased. Note however that in the case of white noise, as $\Gamma_{ww}(\nu)$ is constant, the polarization periodogram is unbiased for any $N$. This is a bivariate counterpart of a classical result, see e.g. [41, p. 202].

Alike in standard spectral analysis [41], data tapers are to be employed to produce a direct spectral estimator with better bias properties than the naive polarization periodogram.

IV-B Multitapering

The multitaper method is a well established technique [41, 42, 43] which produces a spectral density estimate with reduced variance, while maintaining good bias properties. The basic idea is to compute a series of $K$ direct estimators $\hat{\Gamma}_{xx}^{k}(\nu)$, $k=0,1,\ldots K-1$ that are approximately uncorrelated [41]. The $k$-th spectral estimator reads

[TABLE]

where the $h_{k}$’s are real-valued sequences of size $N$. They are normalized $(\sum_{t=0}^{N-1}h_{k}[t]^{2}=1)$ and orthogonal

[TABLE]

Functions satisfying these conditions together with good leakage properties are for instance the Slepian tapers [44]. Then the multitaper estimate is obtained by averaging:

[TABLE]

In practice the value of $K$ often satisfies $K<10$ to maintain a good spectral resolution [41].

IV-C Estimation of the degree of polarization

IV-C1 Theoretical properties

The estimation of the degree of polarization (35) has attracted interest in the signal processing community [45, 46] in relation to many fields [39, 40]. A naive estimator (at frequencies where the polarization periodogram is nonzero) based on the polarization periodogram would be trivial since:

[TABLE]

which is systematically biased, except for frequencies where $x[t]$ is fully polarized. In a situation where $M$ approximately uncorrelated estimates of the spectral density are available (having multiple realizations of $x$ or using a multitaper estimate (65), in which case $M=K$) one can form a new estimate of the degree of polarization as

[TABLE]

which is a better estimator of $\Phi$ than (66). Medkour and Walden222Their approach is based on spectral matrices rather than the quaternion-valued spectral density introduced here. However this does not change the nature of their results, since definitions of the degree of polarization are identical. [45] studied theoretically this estimator in a Gaussian setting and showed that it is unbiased in the limit $M\rightarrow\infty$.

IV-C2 Numerical study

We propose to numerically study its performances. To avoid spectral blurring effects, we consider the (Gaussian) white noise case since in that case the polarization periodogram is an unbiased estimator of the spectral density:

[TABLE]

The UP decomposition (56) of bivariate white noise allows to generate white noise with prescribed polarization properties. We fix $S_{0}=1$ and generate $M$ independent bivariate white Gaussian noise sequences, leading to $M$ independent periodogram estimates of the spectral density (61). The degree of polarization is then estimated by (67). The bias is estimated by Monte-Carlo simulation.

Fig. 2 depicts the bias in the estimation of the degree of polarization, for $M=1,2,5,10,20,50$ and $500$. Given $M$, the bias increases as the true degree of polarization goes to 0. The bias decreases with larger values of $M$, and becomes negligible for $M\rightarrow\infty$. Note that for typical values of $M<10$ admissible in multitaper estimation, the bias remains significant for $\Phi\simeq 0.6$. Our results agree with those of [45].

IV-D Summary

Spectral density estimation of bivariate signals suffers from two biases. As in the univariate case, the naive polarization periodogram $\hat{\Gamma}^{(p)}_{xx}(\nu)$ is a biased estimator of the spectral density since the signal is of finite length $N$. Well-known multitapering techniques can be adapted to handle this bias and reduce variance.

The second source of bias relates to the estimation of polarization properties. Precisely, a key quantity is the degree of polarization, as it relates the ratio between polarized and unpolarized parts of the signal. Spectra of polarization attributes are more difficult to obtain than simple power spectra. They require the observation of many realizations to reach a good accuracy.

V Numerical validation

We consider the synthetic bivariate signal of Section III-G3 $y[t]=x[t]+w[t]$, where $x$ is a bivariate monochromatic signal defined by (46) and where $w$ is a bivariate white Gaussian noise given in (56). All signals are of length $N=1024$. We consider positive frequencies only, as negative frequencies can be obtained by symmetry (17).

The frequency of the monochromatic signal $x$ is set to $\nu_{0}=128/N=0.125$. The signal $x$ has Stokes parameters at frequency $\nu_{0}$: $S_{0,x}(\nu_{0})=1$, $s_{1,x}(\nu_{0})=S_{1,x}(\nu_{0})/S_{0,x}(\nu_{0})=-0.354$, $s_{2,x}(\nu_{0})=S_{2x}(\nu_{0})/S_{0,x}(\nu_{0})=-0.612$, $s_{3,x}(\nu_{0})=S_{3,x}(\nu_{0})/S_{0,x}(\nu_{0})=0.707$ and $\Phi_{x}(\nu_{0})=1$ since $x$ is deterministic. For $\nu\neq\nu_{0}$, all Stokes parameters $S_{\alpha,x}$, $\alpha=0,...,3$ are zero. Equivalently, $x$ is defined by $a_{x}=1$, $\theta_{x}=-\pi/3$ and $\chi_{x}=\pi/8$ using the polar form (46).

The white noise signal $w$ has constant-frequency Stokes parameters (see Section III-G2) such that $S_{0,w}(\nu)=10/N$, $\Phi_{w}(\nu)=0.2$, $\theta_{w}(\nu)=\pi/8$ or equivalently, $s_{1,w}(\nu)=0.141$ and $s_{2,w}(\nu)=0.141$. Since $w$ is a white noise sequence (and thus has no memory), $s_{3,w}(\nu)=0$ for any $\nu$.

The spectral description of $y[t]=x[t]+w[t]$ was derived explicitely in Section III-G3. Using expressions (57) and (58) we see that the resulting Stokes parameters at $\nu\neq\nu_{0}$ are those of $w$. At frequency $\nu_{0}$ we have

[TABLE]

Normalized Stokes parameters and polarization degree then read for frequency $\nu_{0}$

[TABLE]

Due to the polarization properties of $w$, the polarization properties at $\nu_{0}$ of $y$ are not the same as those of $x$.

To investigate the estimation of the spectral density of $y$, we have generated $M=20$ realizations. For each realization, we compute the polarization periodogram and the multitaper estimate using $K=5$ Slepian tapers. To reduce the bias of the degree of polarization estimate, those estimates are averaged to produce a polarization periodogram estimate and a multitaper estimate. This also reduces bias on normalized Stokes parameters estimates, as one expects the estimates of $s_{1},s_{2},s_{3}$ to be biased as they depend on the value of the degree of polarization.

Fig. 3 shows simulation results, and one realization of the process $y$ is shown. As expected the (averaged) multitaper estimate has less variance than the (averaged) polarization periodogram. Thin lines indicate the theoretical values of the components of the spectral density, showing accurate estimation of the spectral quantities.

VI Conclusion

This paper provides a powerful and relevant framework for an interpretable and efficient spectral analysis of stationary bivariate processes. The richness of the quaternion algebra permits a fruitful interplay between mathematical tools and physical features. Using the QFT, we have introduced the quaternion valued-spectral representation of a bivariate stationary random signal. As a result, the quaternion-valued spectral density is defined, which leads naturally to a spectral analysis in terms of Stokes parameters. It permits a direct interpretation of both power and polarization features of the signal. Simple theoretical examples demonstrate the efficiency of the approach. Nonparametric spectral density estimation has been investigated. The limitations of the proposed tools can be studied using standard techniques of univariate spectrum analysis. Moreover, we have stressed the issue raised by the degree of polarization estimation and polarization attributes. These key quantities are relevant to the analysis of bivariate signals but require more care than standard spectral analysis of univariate signals. Our approach is very generic, and generalizes the standard toolbox of spectral analysis to bivariate stationary signals. It paves the way to new developments in the simulation, estimation and filtering of bivariate signals.

Appendix A Proof of the spectral representation theorem 1

The proof is divided in two parts, for clarity.

Existence

Let $x[t]=u[t]+\bm{i}v[t]$, where $u[t],v[t]$ are real-valued, zero-mean, harmonizable stationary processes. These real processes admit a spectral representation, such that

[TABLE]

since the QFT applied to $\mathbb{C}_{\bm{j}}$-valued signals is equivalent to the usual Fourier transform. By linearity of the QFT, the spectral increments of $x$ are $\mathrm{d}X(\nu)=\mathrm{d}U(\nu)+\bm{i}\mathrm{d}V(\nu)$, so that

[TABLE]

holds for all $t$ in the mean-square sense.

Properties of the spectral increments

The properties of the spectral increments $dX(\nu)$ are a direct consequence of the properties of the spectral increments of $u$ and $v$, respectively. If $x$ is assumed zero-mean stationary,

[TABLE]

Turning to the second-order properties of the spectral increments, let us consider the spectral representation of $u$ and $v$. Second-order stationarity implies that (see [7] for details)

[TABLE]

and autocorrelation functions of $u,v$ read

[TABLE]

The quantity $\mathbf{E}\left\{|\mathrm{d}{U}(\nu)|^{2}\right\}$ is interpreted as the spectral density $P_{uu}$ of $u$ times $\mathrm{d}\nu$. The same result holds for $v$.

To fully characterize the spectral increments of $x$, we also need the covariance between the spectral increments of $u$ and $v$. Since $u$ and $v$ are jointly second-order stationary,

[TABLE]

and the cross-correlation function reads

[TABLE]

As a result we have from (74) and (77):

[TABLE]

When $\nu^{\prime}=\nu$, the properties are summarized by the spectral density $\Gamma_{xx}(\nu)$

[TABLE]

which separates in quaternion algebra the information contained in the two moments of the spectral increments. This theorem holds also for quaternion-valued stationary signals by simply adapting the proof. As a corollary, combining (75), (76) and (78) for $\tau=0$ yields (15).