Compressive Covariance Sensing-Based Power Spectrum Estimation of Real-Valued Signals Subject to Sub- Nyquist Sampling

In this work, an estimate of the power spectrum of a real-valued wide-sense stationary autoregressive signal is computed from subNyquist or compressed measurements in additive white Gaussian noise. The problem is formulated using the concepts of compressive covariance sensing and Blackman-Tukey nonparametric spectrum estimation. Only the second-order statistics of the original signal, rather than the signal itself, need to be recovered from the compressed signal. This is achieved by solving the resulting overdetermined system of equations by application of least squares, thereby circumventing the need for applying the complicated l1-minimization otherwise required for the reconstruction of the original signal. Moreover, the signal need not be spectrally sparse. A study of the performance of the power spectral estimator is conducted taking into account the properties of the different bases of the covariance subspace needed for compressive covariance sensing, as well as different linear sparse rulers by which compression is achieved. A method is proposed to benefit from the possible computational efficiency resulting from the use of the Fourier basis of the covariance subspace without considerably affecting the spectrum estimation performance.


Introduction
The efficiency of signal acquisition systems has been greatly improved by the introduction of the concept of compressive sensing (CS) [1]. CS is a technique that enables simultaneous acquisition and compression via sub-Nyquist signal sampling to reduce the costs of data sensing, communications, and storage. The original signal can be recovered by ℓ 1 -minimization provided that the signal is spectrally sparse. Reconstruction of second-order statistics from the compressed measurements is possible without the need for the original signal reconstruction and without the condition of spectral sparsity. Such a technique is referred to as compressive covariance sensing (CCS) [2]. Structure information present in the statistical domain is captured during compression such as the common positive semidefinite Hermitian Toeplitz (HT) structure of the signal covariance matrix. This information enables the reconstruction of wide-sense stationary (WSS) signal statistics by CCS, and the technique has been extended to nonstationary signals as well using online CCS [3]. Compression is assumed to be performed using a linear sparse ruler (LSR) [2,3] that enables subsequent recovery of the Toeplitz covariance matrix by sampling the delays at least once, and the development of the CCS problem leads to an overdetermined system of equations that can be solved by least squares to obtain the covariance values of the original signal.
Many signal processing methods and techniques only require knowledge of second-order statistics of a signal rather than the signal itself. Examples are power spectrum estimation [4], frequency estimation [5,6], and direction-ofarrival estimation in array signal processing [7]. The power spectrum shows the distribution of signal power as a function of frequency. Power spectrum estimation has many practical applications such as in speech analysis where the estimated spectrum helps determine the formants of the vocal tract, in medical diagnosis and electroencephalography where the estimated spectrum helps in the study of sleep disorders for instance, and also in Doppler radar, seismology, and geophysics. In this work, power spectrum estimation is achieved based on CCS to compute or recover the autocorrelation function (simply referred to as covariance) of a zero-mean discretetime signal from compressed measurements. Then, the fast Fourier transform (FFT) of the recovered covariance values is taken to find the power spectrum since the latter constitutes a Fourier transform pair with the signal autocorrelation. Recently, CCS has found application in wideband spectrum sensing for cognitive radio (CR) [4,[8][9][10]. In line with this CR direction, sensing of OFDM signals from sub-Nyquist samples has also been targeted [11]. In CR, the wideband spectrum with subbands is scanned in search of unused bands for transmission exploitation possibilities. Therefore, sensing the power spectrum of the entire frequency band is required. The present work shares many facets of this approach but is more generally directed towards power spectrum estimation of discrete-time signals in which the power spectral density (PSD) estimator is of the type of the standard periodogram and its modified variants [12]. This CCS-related problem can be formulated using covariance matching which, in essence, consists of the estimation of a linearly parameterized covariance matrix from compressed measurements [13,14]; that is, the covariance matrix is assumed to be a linear combination of a set of HT matrices. The subspace spanned by these constituent matrices is referred to as the covariance subspace, and it captures the structural information needed for covariance reconstruction from the available compressed signal.
Without loss of generality, the signals considered in this paper are assumed to be real-valued. Two different bases for the covariance subspace are employed and tested through simulation to obtain the covariance values of the original signal, and subsequently, Blackman-Tukey spectral estimation is applied to find the power spectrum. Performance is measured in terms of normalized mean square spectrum estimation error. The advantages entailed by using any of the two bases are demonstrated regarding accuracy and computational complexity. A method is introduced that allows a reduction in computational complexity using the Fourier basis without considerably degrading the estimation accuracy. Further, to achieve the highest possible compression, a minimal sparse ruler is used for compression [2], and it is shown that two different minimal sparse rulers affect the CCS process differently albeit with the same compression ratio. A range of several compression ratios is also tested.
Recently, various covariance matrix estimation methods from compressed measurements have been introduced. In [15], an algorithm is proposed to recover covariance matrices of high-dimensional signals from random compressive projections based on the projection gradient method. Similarly, compressive projection principal component analysis (PCA) algorithm has been proposed in [16], and the works in [17,18] analyzed the effect of introduced bias on covariance matrices by the projection matrix. The abovementioned compression covariance matrix estimation methods were tested in applications such as hyperspectral imaging [15]. Extending these methods to the context of power spectrum estimation remains a topic for future research. However, it is expected that the present least squares CCS method is more compatible with power spectrum estimation, especially when using the Fourier basis of the covariance subspace which allows the deduction of signal spectral components midway through the computations, thereby reducing the computational burden. Another related work is the implementation of cooperative compressive power spectral estimation [19] exploiting channel diversity to account for received user signal fading in wireless sensor communications. Similar distributed processing methods for power spectral estimation can be found in [20]. The work in [21] proposes a sub-Nyquist wideband spectrum sensing method based on CCS for rapid wideband spectrum sensing that is suitable for fast spectral detection for CR. However, accurate power spectral density estimation is not the primary requirement in this course-grained sensing method.
To the best of the author's knowledge, the most relevant work in the literature to the present work is the power spectrum estimation method presented in [4] in which an approach to the problem is similarly developed based on sub-Nyquist sampling with no requirement of sparsity constraints. The solution in [4] is found by solving a least squares problem and using an LSR for sampling to achieve two objectives. These are spectrum estimation and spectrum detection, where the latter is common in CR applications. The present work deals only with spectrum estimation and has the following contributions. First, the Blackman-Tukey approach of widowing the covariance is used to achieve more accurate spectrum estimation while implementing least squares CCS. Second, and most importantly, a CCS-related method referred to as the improved alpha-based method is devised whose performance in terms of mean squared estimation error approaches that of the original CCS method and has the advantage of reduced computational complexity.
In order to present the experimental methods used in power spectrum estimation based on CCS, the theory of least squares CCS is first explained.

Least Squares Compressed Covariance Sensing (LS-CCS)
In this work, the CCS problem is formulated to recover the covariance of a zero-mean real random signal vector x ∈ R N from the vector y ∈ R M which is a compressed version of x according to where Φ is the M × N compression or sampling matrix. Typically, M < <N and there may be several available realizations of y. The compression ratio is given by M/N. We consider sparse compression matrices Φ with at most one nonzero entry, usually of unit value, in each row or column. The covariance matrix of x, containing the second-order statistics, is given by It is assumed to be a linear combination of generally HT matrices Σ s that are elements of the set I = fΣ 0 , Σ 1 , ⋯, Σ S−1 g ⊂ R N×N and can be written as 2 Modelling and Simulation in Engineering where α s are real scalars.
Since the elements of the vector x are samples from a WSS process, the positive semidefinite HT structure of Σ allows this compression. The subspace I is called the covariance subspace and is a linearly independent set of HT matrices. The elements of I constitute a basis for covariance such that the decomposition of Equation (3) is unique and a knowledge of the scalars α s leads to the determination of Σ. The inequality S ≤ 2N − 1 must be met for I to be linearly independent because Σ s are HT and α s are real.
The covariance matrix of the compressed vector y is Substituting Equations (1) and (2) in Equation (4) where This indicates that Σ is a linear combination of the Hermitian (but not necessarily Toeplitz) matrices in the set If compression preserves all relevant information, then I is linearly independent and α s can be found by applying Equation (6) then (5). Thereafter, Σ can be found from Equation (3). However, if compression is strong, the linear independence of I may be lost and the covariance of x cannot be fully recovered. This method of finding α s is called covariance matching [13,14], and it leads to an overdetermined system of equations (represented by Equation (5)) that can be solved by least squares. Since the matrices Σ are M × M, this system of equations has M 2 equations and only S unknowns, where usually S ≤ 2N − 1 so that it is typically an overdetermined system of equations. The solution is mathematically proved in detail in [2] and by way of a clarifying example in [22]. A flowchart showing the algorithmic steps of the LS-CCS method in sequential order is demonstrated in Figure 1.
The least squares solution for the estimate of Σ, denoted by Σ∧ LS , as arrived at in [2], is Equation (7) mathematically summarizes Figure 1, where the operator vec −1 denotes the restacking of a vector into a square matrix, ⊗ is the Kronecker product, and the dagger † denotes the pseudoinverse. The vector b σ y is the estimate of σ y which is the vectorization (the reverse process of vec −1 ) of the matrix Σ, and the matrix S has the vectors σ s as its columns. The vectors σ s are the vectorizations of the matrices Σ s .
Note that the solution of Equation (7) requires the estimation of the covariance of the compressed signal whose measurements y are available. The effect of different estimation methods for the covariance of the compressed signal Σ on the whole CCS process leading to the recovery of Σ is investigated in [22]. (1) can be explained in terms of a linear sparse ruler. This is a classical problem in number theory.

Linear Sparse Ruler (LSR). The function of the compression matrix in Equation
Compress vector x to obtain vector y using compression matrix Φ : Find the estimate LS of the original covariance matrix Σ: Step 1 Step 2 Step 3 Step 4 Step 5 Step 3 Modelling and Simulation in Engineering the condition that for all lag values τ between 0 and N − 1, there is at least one pair of elements in K with the difference between the two elements equal to τ. It is called a minimal LSR if no other length-ðN − 1Þ LSR has a smaller cardinality M [2]. The LSR can be thought of as a ruler with missing marks, and the remaining M marks allow all integer distances between 0 and N − 1 to be measured. As defined before, the compression ratio is M/N. To compare with the operation of the compression matrix Φ, each of the M rows of Φ contains a "one" at one of the corresponding M marks of the LSR and zeros elsewhere.

The Toeplitz Covariance Subspace.
A Toeplitz matrix is constant along its diagonals. The Toeplitz subspace is the set of all N × N HT matrices, and it is a subspace of C N×N over the real scalar field. Two bases will be considered for this subspace, denoted by the sets I 1 and I 2 [13]. These are given below.
where T k is the HT matrix with ones on the diagonals +k and -k and zeros elsewhere andT k is the HT matrix with the imaginary unit j on the diagonal +k, the imaginary negative unit -j on the diagonal -k, and zeros elsewhere. This basis is called the standard basis. The other important basis for this subspace is the Fourier basis: where Clearly, the number of elements in each of I 1 and I 2 is 2N − 1.
The following section presents the experimental methods for nonparametric power spectrum estimation based on CCS.

Experimental Methods
Nonparametric power spectrum estimation does not require a priori information of the signal model so that the spectral estimation is governed entirely by the observed data. The PSD is the discrete-time Fourier transform (DTFT) of the covariance of the zero-mean signal. The covariance sequence of the N-element vector x is the first column of the Toeplitz covariance matrix and is given by To obtain an estimateĉ x of the above covariance sequence of the original signal x from compressed measurements y by applying LS-CCS, we need an estimate of the covarianceĉ y of the compressed signal y as in Equation (7). In [22], three methods of covariance estimation have been tested resulting in the unbiased covariance estimation method giving the best results. However, for spectrum estimation applications, we use the biased covariance estimation given byĉ where T is the number of terms in the sum, τ is the lag, and The estimate of Equation (12), though biased, is preferred over other estimates because it is nonnegative definite in contrast to the unbiased covariance estimate [23,24]. The original signal covariance estimates are based on the compressed signal covariance estimates through the CCS process. The FFT of the unbiased covariance estimates would exhibit negative values for some frequencies whereas a PSD should be positive for all frequencies. Thus, the unbiased covariance estimate is not related to any possible spectral estimate [24]. In any case, the covariance estimates of small lags are more reliable than those of large lags because the number of summed terms is smaller for the latter covariance values making them less accurate. The Blackman-Tukey method of nonparametric spectrum estimation windows the The estimated PSD is found from the DTFT of the estimated covariance: where PŜD BT ðωÞ is the Blackman-Tukey spectral estimate. In the present work, both the covariance values and the PSD are real and even functions of their respective arguments. Naturally, the window used such as the Bartlett window is also an even function of the lag. In the simulations, the FFT is taken rather than the DTFT, and zero padding of the covariance values is applied to get a better plot when the FFT values are interpolated. Many realizations of the random signal vector x are used to compute the power spectrum and an ensemble average is taken [12].
The above procedure of estimatingĉ x from y by CCS and then applying a Fourier transform is termed compressive power spectrum estimation [13]. A flowchart for the Blackman-Tukey version of this procedure is shown in Figure 2. Constraining the covariance matrix of the original signal to be, in general, HT and positive semidefinite means that the covariance subspace is the Toeplitz subspace. One can employ any of the two bases for this subspace that are explained in Section 2. However, the Fourier basis allows the real scalars α s to represent the power spectrum at frequencies 2πs/ð2N − 1Þ in a straightforward manner. This Fourier transformation result becomes clear upon inspecting Equation (3) with Σ s given by Equation (10). This method will be called the alpha-based method of power spectrum estimation. In this case, there is no need to apply Equation (3) (Step 6 of Figure 1) to find the CCS-estimated covariance matrix from which to find the PSD, which would be the case if the standard basis was used. Thus, the alpha-based method would save N 2 ð2N − 1Þ multiplications. This computational advantage, however, comes at the expense of lowered performance. The reason is that there is no chance of decreasing the effect of erroneous high-lag covariance values if α s are taken directly to represent the spectral coefficients.  Figure 1 S-1 s=0 Figure 3: Flowchart of the alpha-based and improved alpha-based methods of power spectrum estimation using the Fourier basis for the Toeplitz subspace. The dashed boxes concern the proposed improved alpha-based method. In this work, a means to circumvent the lowered performance of the alpha-based method of spectrum estimation is proposed in order to benefit from the ensuing reduced computational complexity. The performance may be improved if the compressed signal estimated covariance is windowed by a Bartlett or triangular window before applying least squares. Figure 3 is a flowchart of the alpha-based and improved alpha-based power spectrum estimation methods. The efficacy of the proposed improved alpha-based method is demonstrated in the following section on results.

Simulation Results and Discussion
Simulations are carried out in MATLAB. A length-21 realvalued signal vector x is taken from a first-order autoregressive random process with parameter -0.8. The autocorrelation function of the process would then be given LS-CCS is then applied to recover the original covariance and then compute the estimated power spectrum by applying a FFT operation. The biased method is used for estimating the covariance of the compressed signal, and the estimated original covariance values are windowed by a Bartlett window before computing the spectrum. The number of elements in both the x-related covariance and power spectrum vectors is 41 or (2N − 1). Both the standard and the Fourier bases of the Toeplitz covariance subspace are tested and found to yield identical results when the spectrum in both cases is computed from the estimated or recovered original covariance. The performance criterion for the spectrum estimator implemented is the normalized mean square spectrum error given by where PSD BT and PŜD BT are the true and estimated Blackman-Tukey length-41 power spectrum vectors, respectively. All simulation results (Figures 4-9) were averaged over 20 realizations of the original AR random vector x, and for   each realization, the results were ensemble-averaged over 1000 runs to account for the noisy compressed observations. Figure 4 shows the NMSE vs. SNR when spectrum estimation is performed by LS-CCS with the three used LSRs. Although K 1 and K 2 are both minimal LSRs, K 1 clearly performs better. Due to the higher compression ratio and hence less compression of K 3 , this LSR gives the least NMSE. When the Fourier basis is used, the computed α s can be considered as the spectral components, which is the alphabased method as discussed in Section 3. However, this results in suboptimal performance when compared with using α s to recover the covariance and then windowing and taking the FFT. If α s are taken directly to represent the spectral coefficients, there would be no chance of windowing the covariance values to improve the spectrum estimation. Figure 5 clarifies this issue of lowered performance of the alphabased method compared to the original LS-CCS method. It is a plot of NMSE vs. SNR for both methods of spectrum estimation using the Fourier basis and for two different compression ratios. It is noticed that when compression is less, that is, a higher compression ratio, the difference between the two approaches to spectrum estimation decreases because, then, α s become a better representation of the spectral coefficients.
As discussed in Section 3 on experimental methods, the proposed improved alpha-based method is tested and found to give comparable results to the LS-CCS method in addition to its being more computationally efficient. The improved performance stems from windowing the compressed signal estimated covariance by using a Bartlett window as the flowchart of Figure 3 demonstrates. The results of this attempt are shown in Figure 6. It is a plot of NMSE vs. SNR using the Fourier basis and the alpha-based methods with a minimal LSR. The performance is significantly improved when the compressed signal estimated covariance is windowed. However, the method of recovering the original signal covariance is still better. The performance of the improved alpha-based method is expected to increase with higher compression ratios. Figure 7 is a plot of the NMSE vs. compression ratio M/N for the three methods of Blackman-Tukey LS-CCS, the alpha-based method, and the improved alpha-based method all using the Fourier basis and for two noise scenarios: SNR = 30 dB and SNR = 5 dB. The five compression ratios considered in the figure are 8/21, 9/21, 13/21, 15/21, and 19/21. The K 1 LSR with ratio 8/21 is used, and then, the larger ratio cases are implemented by randomly adding marks to the LSR. This is the procedure followed in [4] to implement the larger compression ratios. It is clear from the figure that the improved alpha-based method performance compares very well with the original CCS method and guarantees a significant improvement over the simple alpha-based method. The slightly reduced performance of the new method compared to the original CCS method is counterbalanced by the gained advantage of lower computational complexity. As such, the devised method compares favorably with [4] while gaining the advantage of reduced computational burden. Figure 8 shows the Blackman-Tukey estimated power spectrum from compressed measurements by LS-CCS using the Fourier basis and by original signal covariance recovery. The compression ratio is 11/21 employing the K 3 LSR and the SNR is 30 dB. The smooth spectrum is a result of zero padding by 60 zeros and then interpolating. Figure 9 shows the true and recovered original signal covariance values under these conditions.

Conclusion
Power spectrum estimation of an AR random signal is achieved from compressed measurements using LS-CCS to recover the signal covariance and compute the Blackman-Tukey spectral estimate. Different compression minimal LSRs with the same compression ratio do not necessarily perform similarly. Standard and Fourier bases of the covariance  Modelling and Simulation in Engineering subspace yield identical performance. When the Fourier base is used, spectrum estimation from the recovered covariance values gave better (less) mean square spectrum estimation error than the alpha-based method at the expense of increased computational complexity. However, the present work proposes an improvement of the latter method by windowing the compressed-signal covariance estimate needed in the LS-CCS process. In terms of mean square spectrum estimation error, the proposed improved alpha-based method has close performance to the method of estimating the spectrum from recovered covariance by LS-CCS and with the additional advantage of reduced computational complexity.

Data Availability
The datasets used in this study are computer-generated. The computer simulation MATLAB code is available from the author on reasonable request.

Conflicts of Interest
The author declares no conflicts of interest in the publication of this paper.