Compressive Estimation of a Stochastic Process with Unknown Autocorrelation Function

TL;DR

This paper introduces a blind spectral estimation method for predicting stationary Gaussian processes without prior knowledge of the autocorrelation function, achieving performance comparable to traditional MMSE predictors.

Contribution

It proposes a novel spectral estimation approach using atomic-norm minimization for blind prediction of Gaussian processes with unknown autocorrelation.

Findings

The blind predictor performs comparably to the MMSE predictor with known PSD.

Atomic-norm minimization effectively quantizes the spectral content from observed samples.

Simulation results validate the proposed method's effectiveness.

Abstract

In this paper, we study the prediction of a circularly symmetric zero-mean stationary Gaussian process from a window of observations consisting of finitely many samples. This is a prevalent problem in a wide range of applications in communication theory and signal processing. Due to stationarity, when the autocorrelation function or equivalently the power spectral density (PSD) of the process is available, the Minimum Mean Squared Error (MMSE) predictor is readily obtained. In particular, it is given by a linear operator that depends on autocorrelation of the process as well as the noise power in the observed samples. The prediction becomes, however, quite challenging when the PSD of the process is unknown. In this paper, we propose a blind predictor that does not require the a priori knowledge of the PSD of the process and compare its performance with that of an MMSE predictor that has a full knowledge of the PSD. To design such a blind predictor, we use the random spectral representation of a stationary Gaussian process. We apply the well-known atomic-norm minimization technique to the observed samples to obtain a discrete quantization of the underlying random spectrum, which we use to predict the process. Our simulation results show that this estimator has a good performance comparable with that of the MMSE estimator.

Equations

$$y_n=h_n+z_n,n\in [N]$$

$$\hat{\mathbf{g}}=\text{boldmath}{\Sigma }_{\mathbf{g}\mathbf{y}}\text{boldmath}{\Sigma }_{\mathbf{y}\mathbf{y}}^{-1}\mathbf{y}=\text{boldmath}{\Sigma }_{\mathbf{g}\mathbf{h}}\text{boldmath}{\Sigma }_{\mathbf{y}\mathbf{y}}^{-1}\mathbf{y},$$

$$r_h(m)={\int }_{-1/2}^{1/2}{\mathrm{e}}^{\mathrm{j}2\pi m\xi }F(\mathrm{d}\xi )$$

$$r_h(m)={\int }_{-1/2}^{1/2}{\mathrm{e}}^{\mathrm{j}2\pi m\xi }{S}_h(\xi )\mathrm{d}\xi .$$

$$\mathbb{E}[(H(u)-H(v))(H(u^{\prime })-H(v^{\prime }){)}^{\ast }]=0.$$

$$\mathbb{E}[|H(u)-H(v){|}^2]={\mu }_F(v,u]=F(u)-F(v),$$

$$h_n={\int }_{-1/2}^{1/2}{\mathrm{e}}^{\mathrm{j}2\pi \xi n}H(\mathrm{d}\xi )$$

$$h_n^d={\int }_{-1/2}^{1/2}{\mathrm{e}}^{\mathrm{j}2\pi \xi n}{H}^d(\mathrm{d}\xi )=\sum _{i=1}^k{X}_i{\mathrm{e}}^{\mathrm{j}2\pi {\xi }_{i}n},$$

$$h_n^c={\int }_{-1/2}^{1/2}{\mathrm{e}}^{\mathrm{j}2\pi \xi n}{H}^c(\mathrm{d}\xi ),$$

$$\displaystyle r_{h}(m)=\left\{\begin{array}[]{ll}1&m=0,\\ 0&\text{otherwise,}\end{array}\right.$$

$$h_n^c=\int {\mathrm{e}}^{\mathrm{j}2\pi \xi n}{H}^c(d\xi )\approx \int {\mathrm{e}}^{\mathrm{j}2\pi \xi n}{\hat{H}}^c(d\xi ),n\in [N],$$

$${\mu }^{\ast }=$$

$${\mu }^{\ast }=\sum _{i=1}^l{c}_i^{\ast }{\delta }_{{\xi }_i^{\ast }},$$

$${\mathbf{x}}^{\ast }=$$

$$\|{\bf x}\|_{\mathcal{D}}=\inf\Big{\{}\underset{i}{\sum}|c_{i}|:\exists\,\xi_{i}\text{ s.t. }{\bf x}=\underset{i}{\sum}c_{i}{\bf a}(\xi_{i})\Big{\}}$$

$$\{{\mathbf{x}}^{\ast },{\mathbf{u}}^{\ast },{\lambda }^{\ast }\}=$$

$$\text{s.t.}\left[\begin{array}{cc}\text{Toep}(\mathbf{u})& \mathbf{x}\\ {\mathbf{x}}^{\mathsf{H}}& \lambda \end{array}\right]⪰0,\Vert \mathbf{x}-\mathbf{y}{\Vert }_2\le ϵ,$$

$$\{{\mathbf{q}}^{\ast },{\mathbf{H}}^{\ast }\}$$

$$\text{s.t.}\left[\begin{array}{cc}\mathbf{H}& -\mathbf{q}\\ -{\mathbf{q}}^{\mathsf{H}}& 1\end{array}\right]⪰0,\mathbf{H}={\mathbf{H}}^{\mathsf{H}}$$

$$\displaystyle~{}~{}\sum_{k=1}^{N-j}H_{k,k+j}=\left\{\begin{array}[]{ll}1,j=0,\\ 0,j=1,2,\ldots,N-1,\end{array}\right.$$

$$\begin{array}{rl}& Q({\xi }_i^{\ast })=\text{sign}(c_i^{\ast }),i=1,2,\dots ,l\\ & |Q(\xi )|<1,\xi \ne {\xi }_i^{\ast }\forall i\end{array}$$

$${\Delta }_{\xi }=\underset{i\ne j}{\min }|{\xi }_i^{\ast }-{\xi }_j^{\ast }|\ge \frac{1}{\lfloor (N-1)/4\rfloor }.$$

$${\bf A}=\Big{[}{\bf a}(\xi_{1}^{\ast}),{\bf a}(\xi_{2}^{\ast}),\ldots,{\bf a}(\xi_{l}^{\ast})\Big{]}.$$

$${\mathbf{c}}^{\ast }=$$

$$[\hat{\mathbf{g}}{]}_n=\sum _{i=1}^l{c}_i^{\ast }{\mathrm{e}}^{\mathrm{j}2\pi (N+n){\xi }_i^{\ast }}.$$

$${\text{MMSE}}_{\text{genie}}(T)\le \text{MMSE}(T)$$

$${\hat{h}}_{N+T}$$

$${\text{MMSE}}_{\text{genie}}(T)$$

$$=\mathbb{E}[|h_{N+T}^c-{\hat{h}}_{N+T}^c{|}^2]+\mathbb{E}[|h_{N+T}^d-{\hat{h}}_{N+T}^d{|}^2]$$

$$\mathbb{E}[|h_{N+T}^c-{\hat{h}}_{N+T}^c{|}^2]$$

$$=P_c-\mathbb{E}[|{\hat{h}}_{N+T}^c{|}^2]$$

$$\mathbb{E}[|{\hat{h}}_{N+T}^c{|}^2]=\text{boldmath}{\Sigma }_{h_{N+T}^c{\mathbf{h}}^c}\text{boldmath}{\Sigma }_{{\mathbf{h}}^c{\mathbf{h}}^c}^{-1}\text{boldmath}{\Sigma }_{h_{N+T}^c{\mathbf{h}}^c}^{\mathsf{H}}$$

$$\text{boldmath}{\Sigma }_{h_{N+T}^c{\mathbf{h}}^c}(l)$$

Full text

I Introduction

Let ${\mathbb{h}}:=\{h_{n}:n\in{\mathbb{Z}}_{+}\}$ be a zero-mean circularly symmetric stationary Gaussian process. Due to stationarity, the probability law of the process is fully characterized by its autocorrelation function $r_{h}(m)={\mathbb{E}}[h_{n}h^{*}_{n-m}]$. Let $N\in{\mathbb{Z}}_{+}$ and let

[TABLE]

be a window of observations of the process consisting of $N$ samples, where $[N]=\{0,1,\dots,N-1\}$ and where $\{z_{n}:n\in[N]\}$ denotes the observation noise consisting of $N$ i.i.d. zero-mean circularly symmetric Gaussian variables with a variance $\sigma^{2}$. In this paper, we are interested in predicting the process ${\mathbb{h}}$ over the window $\Omega=\{N,N+1,\dots,2N-1\}$ consisting of $N$ future samples of ${\mathbb{h}}$ using the observed samples $y_{n},n\in[N]$. Denoting ${\bf y}=(y_{0},\dots,y_{N-1})^{\sf T}$, ${\bf h}=(h_{0},\dots,h_{N-1})^{\sf T}$, ${\bf z}=(z_{0},\dots,z_{N-1})^{\sf T}$, and ${\bf g}=(h_{N},\dots,h_{2N-1})$, the optimal Minimum Mean Squared Error (MMSE) predictor of ${\bf g}$ from the observations ${\bf y}$ is given by the linear operator

[TABLE]

where \hbox{\boldmath\Sigma}_{{\bf g}{\bf y}}:={\mathbb{E}}[{\bf g}{\bf y}^{{\sf H}}]=\hbox{\boldmath\Sigma}_{{\bf g}{\bf h}} and \hbox{\boldmath\Sigma}_{{\bf y}{\bf y}}={\mathbb{E}}[{\bf y}{\bf y}^{{\sf H}}]=\hbox{\boldmath\Sigma}_{{\bf h}{\bf h}}+\sigma^{2}{\bf I}_{N} denote the cross-correlation matrix of ${\bf g}$ and ${\bf y}$ and the autocorrelation matrix of ${\bf y}$ respectively, and where we used the independence of ${\bf h}$ and the observation noise ${\bf z}$. Notice that the linear MMSE predictor only depends on \hbox{\boldmath\Sigma}_{{\bf g}{\bf h}} and \hbox{\boldmath\Sigma}_{{\bf h}{\bf h}} whose components can be obtained from the autocorrelation function $r_{h}(m)$ of the process ${\mathbb{h}}$. Thus, when the autocorrelation function or equivalently the Power Spectral Density (PSD) of the process is known a priori, the predictor in (2) can be directly computed.

The prediction becomes quite challenging, however, if the autocorrelation function is unknown, which is the case in many applications in communication theory and signal processing. For instance, this problem arises in a wireless communication scenario where one needs to obtain the channel state information, modeled by a Gaussian process, to schedule a set of users, which in turn requires predicting a fading channel in time. The multipath fading is commonly modeled as a superposition of sinusoids with random frequencies, amplitudes and phases [1]. In [2], the authors exploit such a discrete model for the fading channel to estimate the corresponding parameters via ESPRIT, as exploited for Direction of Arrival (DOA) estimation algorithm in [3]. The resulting estimate is then used to extrapolate the process across time. Such a prediction with unknown autocorrelation also arises in the problem of estimating the downlink (DL) channel in a Frequency Division Duplex (FDD) system from the uplink (UL) observations. In such a system, the UL and DL transmissions are performed in different frequency bands, and the Base Station (BS) needs to send pilot signals to the users in the DL and receive their feedback via UL to acquire the channel state information. To reduce the resulting feedback overhead, in [4] the authors proposed an estimator for the delay profile of the channel using the observations over UL, which they exploited to extrapolate the channel state for DL. As in [2], they assume a discrete model (for delay profile), thus, overlooking the possibility of continuous components, which occur in practical scenarios.

In this paper, we generalize the idea proposed by [2, 4] and design a blind predictor that does not require the knowledge of autocorrelation or PSD of the process. To do so, we use the random spectral representation of a stationary Gaussian process whose underlying structure depends on the PSD. Using this representation, we obtain a decomposition of the process ${\mathbb{h}}$ into a discrete and a continuous part denoted by ${\mathbb{h}}^{d}$ and ${\mathbb{h}}^{c}$, which resembles the well-known Wold’s decomposition for stationary processes [5]. We also show that in contrast to ${\mathbb{h}}^{d}$, which can be well predicted from its past samples, prediction of ${\mathbb{h}}^{c}$ typically results in a large error. Motivated by the theory developed in [6, 7, 8, 9], we use an atomic-norm minimization approach to quantize the random spectrum of process ${\mathbb{h}}$, and use the resulting quantization to perform a parametric extrapolation of the process. Compared to the previous works, where atomic-norm minimization is exploited for denoising a sparse signal, in our case the underlying spectrum is not sparse per se and the atomic-norm minimization is used merely for quantization rather than denoising.

The rest of the paper is organized as follows. In sections II and III, we introduce the necessary tools from probability and stochastic processes and state the extrapolation problem to be solved. Section IV describes our spectrum quantization algorithm based on the current results on atomic-norm minimization. Finally, we provide simulation results to evaluate the performance of our proposed algorithm in section V.

II Preliminaries

Earlier we introduced ${\mathbb{h}}$ as a stationary discrete-time circularly symmetric Gaussian process with a zero mean and an autocorrelation function $r_{h}(m)={\mathbb{E}}[h_{n}h^{*}_{n-m}]$. Without any loss of generality, we assume that ${\mathbb{h}}$ has a normalized power, i.e., $r_{h}(0)={\mathbb{E}}[|h_{n}|^{2}]=1$. From the spectral theory of stationary processes [10] and, in particular, the Wiener-Khinchin theorem, we have that

[TABLE]

for some right-continuous and non-decreasing function $F:I\to[0,1]$, with $I:=[-\frac{1}{2},\frac{1}{2}]$, $F(-\frac{1}{2})=0$ and $F(\frac{1}{2})=1$, which is known as the power spectral distribution of ${\mathbb{h}}$ and assigns the positive measure $\mu_{F}(a,b]=F(b)-F(a)$ to any half-open interval $(a,b]$ for $a<b$. When $F$ is dominated by the Lebesque measure (length) over $I$, i.e, $\mu_{F}(A)=0$ for any measurable subset $A$ of $I$ of Lebesgue measure zero, then the derivative $\frac{\mathrm{d}F(\xi)}{\mathrm{d}\xi}$ exists almost everywhere [11], and it is easier to work with the power spectral density (PSD) (also known as the spectral density function) defined as $S_{h}(\xi)=\frac{\mathrm{d}F(\xi)}{\mathrm{d}\xi}$, for which the Wiener-Khinchin theorem takes the more familiar form given by

[TABLE]

The spectral representations in (3) and (4) relate two deterministic functions, i.e., the autocorrelation function and the power spectral density function associated to the Gaussian process ${\mathbb{h}}$. In this paper, we are additionally interested in a random spectral representation of ${\mathbb{h}}$. Such a representation is itself a circularly symmetric Gaussian process ${\mathbb{H}}:=\{H(\xi):\xi\in I\}$ parametrized with $I$. It has independent increments, i.e., for any $v,u,v^{\prime},u^{\prime}\in I$ with $v<u<v^{\prime}<u^{\prime}$,

[TABLE]

Moreover, for such a $v,u$, it satisfies

[TABLE]

which implies that the variance of the increment $H(u)-H(v)$ in ${\mathbb{H}}$ is given by the power of the process ${\mathbb{h}}$ in the spectral interval $(v,u]$. Also, ${\mathbb{h}}$ and ${\mathbb{H}}$ are related via the stochastic integral

[TABLE]

which resembles the Fourier transform relation between the deterministic functions $r_{h}$ and $F$ as in (3) (please refer to [10] for a more rigorous definition of (7)).

The random spectral representation in (7) has an important implication, that is, the process ${\mathbb{h}}$ can be represented as the superposition of discrete-time complex exponentials of the type $\psi_{n}(\xi)=\mathrm{e}^{\mathrm{j}2\pi\xi n}$ for frequencies $\xi\in I$ with circularly symmetric Gaussian coefficients given by the increments $H(\mathrm{d}\xi)$ of the spectral process around every frequency $\xi$. In particular, according to (6), if $F(\xi)$ is continuous at $\xi$ then the coefficient $H(\mathrm{d}\xi)$ is a random variable with infinitesimal variance. In contrast, if $F(\xi)$ has a jump at $\xi$ ($F(\xi^{+})\neq F(\xi^{-})$), then the corresponding complex exponential $\psi_{n}(\xi)=\mathrm{e}^{\mathrm{j}2\pi\xi n}$ has a coefficient of non-infinitesimal variance given by the height of the jump of the spectral distribution at $\xi$, i.e. $F(\xi^{+})-F(\xi^{-})$. A convenient way to represent the jumps in $F(\xi)$ is via using Dirac deltas in the PSD $S_{h}(\xi)$, where the variance of $H(\mathrm{d}\xi)$ at these jump points is given by the coefficient of the delta. Notice that, from Lebesgue decomposition, any Cumulative Distribution Function (CDF) over the real-line and in particular the interval $I$ can be written as a convex combination of a discrete part, a continuous part, and a singular (fractal) part as $F(\xi)=\alpha_{d}F_{d}(\xi)+\alpha_{c}F_{c}(\xi)+\alpha_{s}F_{s}(\xi)$ where $\alpha_{d},\alpha_{c},\alpha_{s}\geq 0$ with $\alpha_{d}+\alpha_{c}+\alpha_{s}=1$. For simplicity, we will focus on the case where $\alpha_{s}=0$ and $F(\xi)=\alpha F_{d}(\xi)+(1-\alpha)F_{c}(\xi)$ for some $\alpha\in[0,1]$. Note that $F_{c}(\xi)$ is a continuous CDF without any jumps, and $F_{d}(\xi)$ is a discrete CDF, i.e., it is right-continuous and piecewise constant with jumps in $I$.

Using the conventions above, we can decompose the stationary process ${\mathbb{h}}$ into a discrete part and a continuous part. The discrete part is equal to the sum of the complex exponentials $\psi_{n}(\xi)$ at points $\{\xi_{i}\}_{i=1}^{k}$ corresponding to the jump points of $F_{d}(\xi)$ given by

[TABLE]

where the variables $\{X_{i}\}_{i=1}^{k}$ are independent zero-mean circularly symmetric Gaussian variables of variance $\alpha(F_{d}(\xi_{i}^{+})-F_{d}(\xi_{i}^{-}))$, and where $H^{d}(\xi)=\sum_{i:\xi_{i}\leq\xi}X_{i}$, for $\xi\in I$, is a spectral Gaussian process with jumps of size $X_{i}$ at location $\xi_{i}$ generated by $\alpha F_{d}(\xi)$ according to (6). Similarly, the continuous part is given by

[TABLE]

where $H^{c}(\xi)$ corresponds to a random spectral process obtained from ${\mathbb{H}}$ after removing the discrete part associated with $\{X_{i},\xi_{i}\}_{i=1}^{k}$ and is generated by the continuous part $(1-\alpha)F_{c}(\xi)$ according to (6). We will use these decompositions in the next section to formulate the prediction problem of the process ${\mathbb{h}}$ more clearly.

III Statement of the Problem

We begin this section with an example. Let ${\mathbb{h}}$ be a stationary process with a purely discrete spectral distribution function $F(\xi)$, that is, $F(\xi)=F_{d}(\xi)$. From the observations we made in section II, we can write $h_{n}=\sum_{i=1}^{k}X_{i}\mathrm{e}^{\mathrm{j}2\pi n\xi_{i}}$, where $\{X_{i}\}_{i=1}^{k}$ are independent zero-mean circularly symmetric Gaussian random variables with variances $\{\sigma_{i}^{2}\}_{i=1}^{k}$, where $\{\xi_{i}\}_{i=1}^{k}$ denote the jump locations in $F_{d}(\xi)$, and where $\sigma_{i}^{2}=F_{d}(\xi_{i}^{+})-F_{d}(\xi_{i}^{-})$. Let $N\in{\mathbb{Z}}_{+}$ be the length of the observation window and let us define the vector-valued function ${\bf a}:I\to{\mathbb{C}}^{N}$ by ${\bf a}(\xi)=(1,\mathrm{e}^{\mathrm{j}2\pi\xi},\dots,\mathrm{e}^{\mathrm{j}2\pi(N-1)\xi})^{\sf T}$. Also, let ${\bf h}=(h_{0},\dots,h_{N-1})^{\sf T}$ be the $N$-dim vector consisting of the first $N$ components of ${\mathbb{h}}$ as before. We can write ${\bf h}=\sum_{i=1}^{k}X_{i}{\bf a}(\xi_{i})$, where it is seen that, for the discrete spectrum $F_{d}(\xi)$, the observation vector ${\bf h}$ resembles the signal received in a uniform linear array with $N$ elements from $k$ Gaussian sources with amplitudes and DoAs $\{X_{i},\xi_{i}\}_{i=1}^{k}$. Thus, provided that $k\ll N$ and the sources are sufficiently separable in $\xi$, the random spectral process can be well estimated from ${\bf h}$ even in presence of noise. The resulting estimate $\widehat{H}(\xi)$ can be used to obtain an estimate of any other sample of the process ${\mathbb{h}}$ as $\widehat{h}_{n}=\int\mathrm{e}^{\mathrm{j}2\pi\xi n}\widehat{H}(d\xi)$, especially those samples ${\bf g}=(h_{N},\dots,h_{2N-1})^{\sf T}$ belonging to our desired prediction window $\Omega=\{N,N+1,\dots,2N-1\}$. Thus, it seems that under some regularity conditions on $F_{d}(\xi)$, the prediction over $\Omega$ is feasible to do.

Now assume that $F(\xi)=F_{c}(\xi)$ consists of only a continuous part. For simplicity of illustration, let us take $F_{c}(\xi)=\xi+\frac{1}{2}$. We can check that $F_{c}$ induces a uniform measure over $I$. In particular, for such a uniform measure

[TABLE]

which implies that ${\mathbb{h}}$ is a white noise consisting of i.i.d. Gaussian samples. As a result, the samples ${\bf g}$ inside the window $\Omega$ can not be estimated from the observed samples. This turns out to be true for continuous distributions $F_{c}(\xi)$ other than the uniform one.

In practice, $F(\xi)$ consists of both continuous and discrete parts. However, in terms of prediction over $\Omega$, we can only hope to predict the discrete part ${\bf g}^{d}$ of ${\bf g}={\bf g}^{d}+{\bf g}^{c}$ over the window $\Omega$ from the discrete part ${\bf h}^{d}$ of the observation ${\bf h}={\bf h}^{d}+{\bf h}^{c}$. In fact, we can apply the estimation technique we mentioned for the discrete part ${\bf h}^{d}$ to estimate $\widehat{H}_{d}(\xi)$, which we can exploit to predict ${\bf g}^{d}$. This results in a prediction error comparable to that of the MMSE predictor, especially, for a long-term prediction. However, an obstacle to do this is that we have access only to a mixture of the discrete part ${\bf h}^{d}$ and the continuous part ${\bf h}^{c}$ through ${\bf h}={\bf h}^{d}+{\bf h}^{c}$ rather than ${\bf h}^{d}$ itself. In particular we have the following result.

Theorem 1

Let ${\bf h}=(h_{0},\dots,h_{N-1})^{\sf T}$ be a vector consisting of $N$ samples of a zero-mean circularly symmetric stationary Gaussian process ${\mathbb{h}}$ with the power spectral distribution $F(\xi)=P_{c}F_{c}(\xi)+(1-P_{c})F_{d}(\xi)$ with $0\leq P_{c}\leq 1$. Let $\text{MMSE}(T)$ denote the minimum mean squared error for estimating sample $h_{N+T}$ of process ${\mathbb{h}}$ given the observation vector ${\bf h}$. Then $\mathop{\rm lim\ inf}_{T\rightarrow\infty}\text{MMSE}(T)\geq P_{c}$. $\square$

Proof:

The proof is provided in Appendix I.     \sqcap$$\sqcup

In the next section, we propose an algorithm that applies the atomic-norm minimization to the full observation ${\bf h}$ to obtain a quantized estimate $\widehat{H}(\xi)$ of the random spectral process $H(\xi)=H^{d}(\xi)+H^{c}(\xi)$. The main idea is that under suitable regularity conditions on the discrete part $F_{d}(\xi)$, we can decompose the resulting estimate as $\widehat{H}(\xi)=\widehat{H}^{d}(\xi)+\widehat{H}^{c}(\xi)$, where $\widehat{H}^{d}(\xi)$ provides a quite precise estimate of the true process $H^{d}(\xi)$, which we can exploit to predict the discrete part ${\bf g}^{d}$. In contrast, $\widehat{H}^{c}(\xi)$ is only a discrete approximation of $H^{c}(\xi)$. In particular, although $\widehat{H}^{c}(\xi)$ provides a good approximation of ${\bf h}^{c}$ over the observation window, i.e.,

[TABLE]

from our earlier explanation, we expect that the long-term prediction of ${\bf g}^{c}$ obtained via $\widehat{H}^{c}(\xi)$ be uncorrelated with ${\bf g}^{c}$. Overall, using the proposed algorithm we hope to predict ${\bf g}$ from observations ${\bf h}$ with an error that grows by twice the variance of the continuous part ${\bf g}^{c}$, thus, twice the best prediction error achieved by the MMSE predictor.

IV Proposed algorithm

As motivated in section III, our goal is to quantize the spectrum of the process ${\mathbb{h}}$. To do so, we find among all the complex measures that fit our observation vector ${\bf h}$, the one with lowest total-variation norm. More precisely, we consider

[TABLE]

where $\mathcal{F}_{N}$ is a linear map returning the first $N$ frequency coefficients of the measure $\mu(\mathrm{d}\xi)$, i.e., $[\mathcal{F}_{N}\mu]_{n}=\int_{-1/2}^{1/2}\mathrm{e}^{\mathrm{j}2\pi n\xi}\mu(\mathrm{d}\xi),~{}n\in[N]$, and $\epsilon$ is an estimate of the $\ell_{2}$-norm of the additive Gaussian noise vector. Optimization (14) returns a discrete measure of the form

[TABLE]

where $\delta_{\xi_{i}}^{\ast}$ denotes a delta measure at $\xi_{i}^{\ast}\in I$. Intuitively speaking, in the noiseless case, minimizing the total-variation norm returns a consistent measure, i.e., we expect that $\xi_{i}^{\ast}$s in (15) belong to the support of the PSD. Considering the complex measure $\mu^{\ast}$ in (15), we have that $\mathcal{F}_{N}\mu^{\ast}=\sum_{i}c_{i}^{\ast}{\bf a}(\xi_{i}^{\ast})$ and total-variation norm of this measure is simply given by $\sum_{i}|c_{i}^{\ast}|$. This implies that (14) can be equivalently written as

[TABLE]

where $\|{\bf x}\|_{\mathcal{D}}$ denotes the atomic norm of ${\bf x}$ over the dictionary \mathcal{D}=\Big{\{}{\bf a}(\xi)\in\mathbb{C}^{N},\xi\in I\Big{\}} defined by

[TABLE]

and where the optimal solution ${\bf x}^{\ast}$ of (16) is given by ${\bf x}^{\ast}=\mathcal{F}_{N}\mu^{\ast}$ in terms of the optimal solution $\mu^{\ast}$ of (14). From the solution of (16) one can determine the frequencies $\xi_{i}^{\ast}$ and their corresponding coefficients $c_{i}^{\ast}$, which immediately specifies the measure $\mu^{\ast}$, showing that (16) and (14) are equivalent.

In particular, solving (16) serves two purposes: first, it reduces the effect of white additive noise and second, it forces the estimated vector to have a sparse representation in $\mathcal{D}$, which implies a sparse quantization of the PSD. Fortunately, compared to (14) which requires an infinite-dimensional optimization over the space of complex measures, (16) has an equivalent SDP form as

[TABLE]

where $\text{Toep}({\bf u})$ denotes a Toeplitz hermitian matrix with ${\bf u}$ as its first column. Once this optimization problem is solved, one can specify the active atoms in $\mathcal{D}$, or equivalently the discrete measure $\mu^{\ast}$, by solving the dual problem to (16), which can also be represented as the following SDP [8]:

[TABLE]

Using the optimal solution ${\bf q}^{\ast}$ of the dual problem, we can construct the dual polynomial $Q(\xi)={\langle{\bf q}^{\ast},{\bf a}(\xi)\rangle}=\sum_{n=0}^{N-1}q_{n}^{\ast}\mathrm{e}^{-\mathrm{j}2\pi n\xi}$ which satisfies the following [8]

[TABLE]

where $\{c_{i}^{\ast}\}_{i=1}^{l}$ are the coefficients of the optimal complex measure $\mu^{\ast}$ in (15) and where $\text{sign}(c_{i}^{\ast})=\frac{c_{i}^{\ast}}{|c_{i}^{\ast}|}$. Therefore, the support $\{\xi_{i}^{\ast}\}_{i=1}^{l}$ of $\mu^{\ast}$ can be obtained from solving the equation $|Q(\xi)|=1$. If there is no noise and ${\bf h}^{c}=\bf 0$, solving $|Q(\xi)|=1$ results in accurate spectrum localization provided that $k\leq\frac{N}{2}$ and $\xi_{i}^{\ast}$s are separated enough such that

[TABLE]

However, this optimization applied to a mixture of discrete and continuous components, generates additional spectral elements to compensate for the the continuous component as well as noise. Fig. 1 illustrates this fact more clearly. In this figure, we have plotted the corresponding dual polynomial in an experiment, where we observed a process consisting of a mixture of discrete and continuous components with $N=64$ and a PSD function which consists of a discrete part located at $\xi=-0.4$ and $\xi=-0.2$ and a continuous part uniformly distributed over $[0.2,0.3]$. Fig. 1 illustrates how $\{\xi_{i}^{\ast}\}_{i=1}^{l}$ are identified via the dual polynomial, hence resulting in spectrum quantization. To obtain the coefficients corresponding to the set $\{\xi_{i}^{\ast}\}_{i=1}^{l}$, we form the matrix of active atoms as

[TABLE]

Knowing the active elements, the atomic-norm minimization in (16) can be written as an $\ell_{1}$-norm minimization

[TABLE]

where $\|{\bf c}\|_{1}\coloneqq\sum_{i=1}^{l}|c_{i}|$. This optimization yields the corresponding coefficients ${\bf c}^{\ast}=(c_{1}^{\ast},\ldots,c_{l}^{\ast})^{\sf T}$. Finally, we use the estimated parameters $\{c_{i}^{\ast}\}_{i=1}^{l}$ and $\{\xi_{i}^{\ast}\}_{i=1}^{l}$ to estimate ${\bf g}$ as

[TABLE]

In the next section, we will show that the proposed method has a good performance in different scenarios.

V Simulation Results

In this section, we provide simulation results to compare the performance of our proposed extrapolation method to that of the MMSE extrapolator. We focus on cases in which the process consists of a discrete and a continuous part. For simulations, we consider a noisy observation vector ${\bf y}$ of size $N=64$ with SNR = $20$ dB. Using this observation, we predict vector ${\bf g}$ consisting of the next $N=64$ samples of the process. Our performance metric in these simulations is the normalized mean squared error defined as $E=\frac{1}{N}{\mathbb{E}}\{\|\hat{{\bf g}}-{\bf g}\|^{2}\}$.

V-A Effect of the continuous component

As before, we assume that $\int F(\mathrm{d}\xi)=1$ and denote the fraction of power in the continuous part by $P_{c}\in[0,1]$. The PSD function is chosen to have a discrete part located at $\xi=-0.4$ and $\xi=-0.2$ and a continuous part uniformly distributed over $[0.05,0.15]$, as plotted in Fig. 2. We increase $P_{c}$ from 0 to 0.5 and for each $P_{c}$ we generate ${\bf y}$, estimate ${\bf g}$ via both the MMSE estimator and our proposed estimator, and calculate the normalized estimation error $E$. We find an estimate of $E$ by averaging it over 1000 independent realizations of the process. Fig. 3 illustrates $E$ as a function of $P_{c}$. It is seen that, the MMSE estimator has an error which is converging to $P_{c}$. The error of our proposed estimator is approximately twice the MMSE and, in fact, is quite low when $P_{c}$ is small.

V-B Effect of increasing the number of jumps in $F_{d}(\xi)$

To see how the number of jumps in the spectral distribution affects the performance of our estimator, we consider a PSD with a continuous part as before, i.e., uniformly distributed over $[0.05,0.15]$ and with a fixed power $P_{c}=0.3$. We add to this PSD $k$ Dirac deltas with equal amplitudes $\frac{1-P_{c}}{k}=\frac{0.7}{k}$ and random frequencies $\{\xi_{i}\}_{i=1}^{k}\subset[0,1]$ with minimum separation larger than $\frac{1}{N}$. We perform the experiment for $k=1,\ldots,5$, and for each $k$ we repeat it for 1000 trials. Fig. 4 illustrates the result of this experiment. As the number of discrete elements grows, the error of the proposed method increases with a slight slope. This is because the power is distributed over a greater number of spikes and the noise effect appears relatively stronger. However, the error is still less than twice the MMSE.

VI Conclusion

Using the framework of total-variation norm and atomic-norm minimization, we proposed a method to quantize the random spectrum of a stationary stochastic process. This quantization is then exploited to predict the process. We investigated the empirical performance of our proposed algorithm via numerical simulations. We illustrated that the prediction error is relatively low, and is roughly proportional to the MMSE up to a factor of 2 when there exists a continuous component in the PSD with non-negligible power.

VII Appendices

VII-A Proof of Theorem 1

The process ${\bf h}$ can be decomposed into independent discrete and continuous components, ${\bf h}^{d}$ and ${\bf h}^{c}$ with corresponding power spectral densities $(1-P_{c})F_{d}(\xi)$ and $P_{c}F_{c}(\xi)$, respectively. Suppose there exists a genie-aided estimator which has access to the observation vectors ${\bf h}_{c}$ and ${\bf h}_{d}$ separately, as opposed to the conventional MMSE considered so far, which has access only to their sum, that is ${\bf h}={\bf h}_{c}+{\bf h}_{d}$. We denote by $\text{MMSE}(T)$ the minimum mean squared estimation error for sample $h_{N+T}$. We mention the mean squared error of the genie-aided estimator by $\text{MMSE}_{\text{genie}}(T)$. According to data processing inequality

[TABLE]

for all $T\geq 0$. Furthermore, the genie-aided estimator can be written as

[TABLE]

and since the continuous and discrete parts are independent, MMSE of the continuous part is only dependent on the continuous part and similarly MMSE of the discrete part is only dependent on the discrete part. As a result, the genie-aided estimator has the error

[TABLE]

In addition, it can be easily shown that

[TABLE]

where

[TABLE]

Note that \hbox{\boldmath\Sigma}_{h_{N+T}^{c}{\bf h}^{c}} is an $N\times 1$ vector whose $l\textsuperscript{th}$ component $l=0,\ldots,N-1$ is given by

[TABLE]

where $S_{c}=\frac{dF}{d\xi}$ denotes the Radon-Nikodym derivative of the absolutely continuous measure $F_{c}$ with respect to the Lebesgue measure [11]. Since $\int|S_{c}(\xi)|d\xi=\int S_{c}(\xi)d\xi=P_{c}<\infty$, from Rimann-Lebesgue lemma [12], it results that

[TABLE]

This implies that \lim_{T\rightarrow\infty}\hbox{\boldmath\Sigma}_{h_{N+T}^{c}{\bf h}^{c}}=\bf 0 and as a result $\lim_{T\rightarrow\infty}{\mathbb{E}}[|\hat{h}_{N+T}^{c}|^{2}]=0$. From (28) it results that $\lim_{T\rightarrow\infty}{\mathbb{E}}[|h_{N+T}^{c}-\hat{h}_{N+T}^{c}|^{2}]=P_{c}-\lim_{T\rightarrow\infty}{\mathbb{E}}[|\hat{h}_{N+T}^{c}|^{2}]=P_{c}$, which implies

[TABLE]

Plugging this inequality in (25) we obtain

[TABLE]

This completes the proof. $\blacksquare$