Decimative Spectral Estimation with Unconstrained Model Order

This paper presents a new state-space method for spectral estimation that performs decimation by any factor, it makes use of the full set of data and brings further apart the poles under consideration, while imposing almost no constraints to the size of the Hankel matrix (model order), as decimation increases. It is compared against two previously proposed techniques for spectral estimation (along with derived decimative versions), that lie among the most promising methods in the field of spectroscopy, where accuracy of parameter estimation is of utmost importance. Moreover, it is compared against a state-of-the-art purely decimative method proposed in literature. Experiments performed on simulated NMR signals prove the new method to be more robust, especially for low signal-to-noise ratio.


Introduction
Various applications in the field of digital signal processing, including speech processing [1] as well as spectroscopy, that is, quantification of NMR signals, are employing complex damped sinusoidal models in order to represent a signal as a sum of exponentially damped complex-valued sinusoids [2][3][4][5][6][7][8]. The generalized model we use is given by where p is the number of complex damped sinusoids that comprise the measured signal. The objective is to estimate the frequencies f i , damping factors d i , amplitudes b i , and phases ϕ 0 + ϕ i , i = 1, . . . , p. ϕ 0 is the zero order phase, whereas ϕ i represents extra degrees of freedom.
The new method proposed here is called DESE D (DEcimative Spectral Estimation by factor D), which can perform decimation by any factor and exploits the full data set, whereas it is not obliged to reduce the size of the Hankel matrix as D increases, allowing the use of size N/2 approximately. The advantage of DESE D relies on the fact that it can benefit from the higher pole resolution obtained by decimation [9], while at the same time is not bound to use smaller sizes of Hankel matrices, as other decimative approaches are. The new method makes use of Singular Value Decomposition. The DESE D is a generalization of the DESE2 method proposed in [10], which performs decimation by factor 2.
The new method has been tested and compared to TLS-ESPRIT and LS-ESPRIT, that lie among the most promising methods for parameter estimation solving the same overdetermined system of equations in a total least squares and least squares sense, respectively. Moreover, their decimative versions are being presented and compared with DESE D for the same decimation factors. In addition, the new method has been tested against a purely decimative method existing in the literature. In the sections that follow the proposed DESE D method as well as the methods against which it is tested are presented and the superior performance of DESE D is shown through Monte-Carlo-based experiments. This can be explained from signal processing theory where it is proved that decimation increases spectral resolution as it brings the frequencies of the sinusoids further apart. The MonteCarlo technique is used to alleviate the random effects of the noisy poles. It is also expected the improvement due to decimation to be more important at low signal-to-noise ratio SNR.
Note that the subspace estimation method TLS-ESPRIT [11] that we are using does not act on the covariance matrix but on the corresponding data matrix, the latter presenting some numerical advantages. Similarly, we are using in this paper the LS-ESPRIT method [12], which is also acting on the data matrix in a Least Squares sense (instead of Total Least Squares). The same method in NMR literature is known as HSVD (Hankel Singular Value Decomposition) [13]. We use the names TLS-ESPRIT and LS-ESPRIT as more suitable for the signal processing the literature. Also note that the conventional decimation method we are referring to as CONDE D [9,14,15] can be also seen as a decimated version of the ESPRIT. However, the TLS-ESPRIT and LS-ESPRIT decimate in a different way than CONDE D, according to the principles introduced for the DESE D method. Section 2.1 introduces DESE and presents a derivation for the decimation factor D case. Section 2.2 contains the algorithmic presentation for DESE D, while special cases are discussed in Section 2.3. The relation between DESE D and a previously proposed conventional decimative spectral estimation method (CONDE D) along with its algorithmic presentation is described in Section 2.4. In Section 3 aspects regarding the decimative versions of two existing spectral estimation methods the TLS-ESPRIT and LS-ESPRIT methods are presented. More specifically, the TLS-ESPRIT algorithm is briefly presented in Section 3.1, while its decimative version called TLS-ESPRIT D is shown in Section 3.2. Similarly, Section 3.3 presents the LS-ESPRIT algorithm and Section 3.4 its decimative version, called LS-ESPRIT D. In Section 3.5 computational considerations for DESE D versus the other methods are discussed. Experimentation and results are found in Section 4, while concluding remarks follow in Section 5.

Derivation.
In [10] we have presented a derivation for the DESE 2 (decimation factor 2) method. A different derivation by employing the well-known Vandermonde decomposition as well as generalization to the decimation factor D case is presented in this paper.
Let S H be the L × M Hankel signal observation matrix of our deterministic signal of p exponentials s(n), n = 0, . . . , N − 1

Theorem 1.
Assuming that the signal s in (1) has unique poles (multiplicity one) and it is noise free, there is an (L − D) order matrix X, such that, and all the signal's decimated poles (i.e., the poles of the signal multiplied by the factor D), are equal to the nonzero eigenvalues of X. A solution of (3) is given by X = S↓ D pinv(S↑ D ) and contains the decimated poles of the signal where pinv denotes the Moore-Penrose pseudoinverse.
Proof. The first claim of the theorem is true because S H is constructed of p complex damped sinusoids and therefore any row of S↓ D can be written as a linear combination of the rows of S↑ D .
We want then to prove that p of the eigenvalues of the matrix are equal to the decimated signal pole estimates z D i , i = 1, . . . , p.
Let us consider the well-known Vandermonde decomposition of S H : where superscript T denotes transpose and the L × p matrix A, p × p matrix G and M × p matrix B are defined as follows: We can then easily write S↓ D = A↓ D GB T and S↑ D = A↑ D GB T where A↓ D and A↑ D are defined from A, similarly to the way S↓ D and S↑ D are defined from S H . Hence, (3) can be written as The latter system of linear equations has an infinite number of solutions given by where superscript H denotes Hermitian conjugate and Δ H A↑ D = 0. Let us consider now the p × p diagonal matrix Φ D containing the decimated signal poles z D i , i = 1, . . . , p. It is easy then to see that denotes the Frobenius norm), the minimum-norm solution to (7) which we compute is Computational and Mathematical Methods in Medicine 3 Furthermore, knowing that is straightforward to see that X 0 = WY and Φ D = YW.
According to the theorem in [16] the nonzero eigenvalues of WY are equal to the nonzero eigenvalues of YW. Consequently, the p nonzero eigenvalues of X 0 (rank(A↑ D ) = p) are equal to the eigenvalues of Φ D , that is, the decimated poles.
The novelty of the above theorem relies on the generalization of the simpler problem for D = 1 as described in [17] to any factor D. In general, decimation may introduce aliasing effects that one should take into account in the algorithm. This is easy to do when prior knowledge (exact or approximate) for the frequencies of the complex damped sinusoids is available. In this case, one can undo the effects of aliasing that might occur when high decimation is used by employing filtering techniques (as described in [18]) to extract the useful sinusoids prior to estimating their parameters. Alternatively, when the frequencies are clustered together one faces the so-called "high resolution" scenario described in [9].
The Proposed Method in the Presence of Noise. In case of real life signals (i.e., signals impaired with additive noise) the peaks are embedded in noise and the rank of matrix S H is full. The equality (3) does not hold any longer because the signal does not obey linear models. If the number of complex peaks to estimate is p, the matrix S↑ D can be enhanced by reducing its rank to p.
To do so we employ the SVD of S↑ D and we retain the p largest singular values based on the assumption that the noise energy is lower than the energy of the p sinusoids of the signal s. The resulting matrix S↑ D e has rank p. Then X is computed from XS↑ D e ≈ S↓ D which gives rise to an overdetermined system of equations with the following solution X = S↓ D pinv(S↑ D e).
Note that since matrix S↑ D e has rank p, X is also of rank p (minimum of ranks of S↑ D e and S↓ D ) and this guarantees that only p of the eigenvalues of X is nonzero and corresponds to the decimated signal pole estimates. This yields the desired decimated estimates of frequencies and damping factors from the angles and magnitudes, respectively, of the eigenvalues of X. These decimated estimates are converted to their nondecimated equivalents f i (frequency estimates) and d i (damping factor estimates) and a computation in a total least squares sense of estimates g i then takes place. In this way, complex-valued linear parameter estimates of g i are calculated, from which amplitude b i and phases ϕ 0 + ϕ i estimates are determined as the magnitudes and angles of g i , respectively. The proposed algorithm, involves the following five steps.

DESE D Algorithmic Presentation. Let
Step 1 (DESE D). We compute the L × M matrix S H of (2) from the N data points s(n) of (1).
Step 2 (DESE D). We compute the S↓ D and S↑ D as the D order lower shift (top D rows deleted) and the D order upper shift (bottom D rows deleted) equivalents of S H .
The best results are obtained when we use the (L−D)×M matrices S↓ D and S↑ D as square as possible [11,[19][20][21].
Step 3 (DESE D). We compute the enhanced version S↑ D e of S↑ D in the following way. We employ the SVD of S↑ D , S↑ D = U↑ D Σ↑ D V ↑ H D and we truncate to order p by retaining only the largest p singular values.
The eigenvalues λ i of X give the decimated signal pole estimates, which in turn give the estimates for the damping factors and frequencies of (1).
Step 5 (DESE D). The last step is to compute the phases and the amplitudes. This is done by finding a least squares solution to (1), with z i replaced by the estimates and s(n) given by the signal data points.
Matrix X of Step 4 in the above described version of DESE D, is computed in a least squares sense. We could however, compute matrix X in a total least squares sense using the Theorem 3.10 presented in [17]. We can, hence, obtain the DESE D TLS method, presented in [22], where the obtained results suggest that DESE D and DESE D TLS perform rather similarly for small noise standard deviations, whereas DESE D seems slightly more robust for large noise standard deviations than its total squares counterpart. This is the reason why only DESE D was included in the experimentation reported here.

DESE D Special Cases.
The above presented method can also serve as a state-space method for spectral estimation, if seen and implemented with no decimation whatsoever (D = 1). In this case, matrices S↓ 1 and S↑ 1 are, respectively, the first-order lower shift (top row deleted) and first-order A variation of such a nondecimative method, called CSE, was proposed in [23]. In this case both matrices S↓ 1 and S↑ 1 (of Step 2) were enhanced (truncated to order p) with the use of SVD. Thus, Steps 3 and 4 presented above are replaced by the following step.
Step 3 (a). We compute the enhanced version S↓ 1 e of S↓ 1 in the following way. We employ the SVD of S↓ 1 , S↓ 1 = U↓ 1 Σ↓ 1 V ↓ H 1 and we truncate to order p by retaining only the largest p singular values.
In the same way, we compute the enhanced version S↑ 1 e of S↑ 1 .

Computational and Mathematical Methods in Medicine
The eigenvalues λ i of X give the signal pole estimates, which in turn give the estimates for the damping factors and frequencies of (1). Note that the CSE method was also proven to be more robust in terms of bad runs when compared to TLS-ESPRIT; however, compared to DESE 1, similar results were obtained while the complexity was increased due to the second enhancement.
When only one enhancement is performed (to matrix S↑ 1 ), the nondecimative method (D = 1, for DESE D) is identical to a method proposed in [24], the MATPEN method.
The eigenvalues λ i of X give the signal pole estimates, which in turn give the estimates for the damping factors and frequencies of (1).

DESE D versus Other Decimative Methods.
The drawbacks of conventional decimative methods are related to the size of the data set and to the overdetermined model order that can be used. Already proposed decimative methods (e.g., [9]), even though they make use of the full data set available, are obliged to reduce the maximum possible matrix size as D increases. Hence, they relate the size of the Hankel matrix n with D, according to n = N/(2D).
This implies that the efficiency of the overdetermined model is reduced. On the contrary, DESE D does not present this drawback and allows the use of matrix size n = (N + 1)/2 − D/2, that change very slowly with respect to decimation factor D.
The DESE D has been tested against a decimation method proposed in [9,14,25], which we call below CONDE D (CONventional DEcimative method for decimation factor D).
The method makes use of the auto-and cross-covariance matrices of the input signal, and decimated sequences of the input signal. Then, averaged covariance matrices are used for parameter estimation of the complex damped sinusoids. Next, it employs Singular Value Decomposition of the resulting matrix to truncate to order p and proceeds with estimation of the frequency and damping factor in a total least squares sense.
The method's algorithmic presentation for decimation factor D involves the following five steps.
Step 1 (CONDE D). We compute the L × M (L = M = N/(2D)) Hankel matrix C k that corresponds to the kth decimated signal, c k (n) = s(k : D : N), where D is the decimation factor, from the N data points s(n) of (1).
Step 2 (CONDE D). We compute a global matrix C by concatenating C k , k = 1, . . . , D as shown below: We then compute a global covariance matrix R = CC H .
Step 3 (CONDE D). We compute the eigen analysis of R = UΛU H to deduce U, which in turn is truncated to order p, thus, yielding U p .
Step 4 (CONDE D). We compute the solution Q of U↑ p Q = U↓ p , in a total least squares sense, where U↓ p (U↑ p ) are derived from U p by deleting its top (bottom) row. The eigenvalues λ i of Q give the decimated signal pole estimates, which in turn give the estimates for the damping factors and frequencies of (1).
Step 5 (CONDE D). The last step is to compute the phases and the amplitudes. This is done by finding a least squares solution to (1), with z i replaced by the estimates and s(n) given by the signal data points.

Decimative Versions of TLS-ESPRIT and LS-ESPRIT
The new concept of using all data samples available while practically imposing no constraints between the size of HANKEL matrix and decimation factor D, included in the DESE D method, can be implemented in other state-space spectral estimation methods, thus, deriving a new family of methodologies.
The subsections that follow present the TLS-ESPRIT and LS-ESPRIT methods along with their decimative versions TLS-ESPRIT D and LS-ESPRIT D, respectively.

The TLS-ESPRIT Algorithm.
The TLS-ESPRIT method, reported in [11], consists of using the Hankel matrix, performing an SVD decomposition and reducing the size of matrices to order p. Damping factors d i and the frequencies f i are estimated in a total least squares sense. Phases and amplitudes are estimated using the least squares method.
Step 1 (TLS-ESPRIT). We compute the SVD of the L × M Hankel matrix S H of (2) from the N data points s(n) of (1): where L ≤ M. The best results are obtained when we use L = M(+1) = N/2.
Step 2 (TLS-ESPRIT). We truncate U, Σ, V to order p and compute Step 3 (TLS-ESPRIT). We compute the solution Q of U↑ p Q = U↓ p , in a total least squares sense, where U↓ p (U↑ p ) are derived from U p by deleting its top (bottom) row. The eigenvalues λ i of Q give the signal pole estimates, which in turn give the estimates for the damping factors and frequencies of (1).
Step 4 (TLS-ESPRIT). The last step is to compute the phases and the amplitudes. This is done by finding a least squares solution to (1), with z i replaced by the estimates and s(n) given by the signal data points. It is worth noting that TLS-ESPRIT and CONDE 1 (no decimation whatsoever) are identical.
Computational and Mathematical Methods in Medicine 5

The TLS-ESPRIT D Algorithm. By using the notion introduced by DESE D that implies minor reduction of the
Hankel matrix size with respect to decimation factor D, the decimative version of TLS-ESPRIT can be easily derived. More precisely, appropriate formation of matrix U p -similar to that of matrix S H in the proposed DESE D-, creates the TLS-ESPRIT D decimative version for decimation factor D.
Its algorithmic presentation is shown below. Please note that Steps 1and 2 are identical in the two approaches.
Step 1 (TLS-ESPRIT). We compute the SVD of the L × M Hankel matrix S H of (2) from the N data points s(n) of (1): where L ≤ M. The best results are obtained when we use L = M(+1) = N/2.
Step 2 (TLS-ESPRIT). We truncate U, Σ, V to order p and compute: Step 3 (TLS-ESPRIT D). We compute the solution Q of U↑ D Q = U↓ D , in a total least squares sense, where U↓ D (U↑ D ) are derived from U P by deleting its top D (bottom D) rows, respectively. The eigenvalues λ i of Q give the decimated signal pole estimates, which in turn give the estimates for the damping factors and frequencies of (1).
Step 4 (TLS-ESPRIT D). The last step is to compute the phases and the amplitudes. This is done by finding a least square solution to (1), with z i replaced by the estimates and s(n) given by the signal data points. Note that in [26], a different decimated version of TLS-ESPRIT (=HTLS in NMR literature) is presented, which is not treated in this paper.

The LS-ESPRIT Algorithm.
If instead of computing in a total least squares sense the solution Q of U↑ p Q = U↓ p , one employs the least squares solution, one uses LS-ESPRIT [12]. In this case, Step 3 TLS-ESPRIT of the Section 3.1 above is replaced by the following.
Step 3 (LS-ESPRIT). We compute the solution Q of U↑ p Q = U↓ p , in a least squares sense, where U↓ p (U↑ p ) are derived from U p by deleting its top (bottom) row. Hence, Q = pinv(U↑ P )U↓ P .
The eigenvalues λ i of Q give the signal pole estimates, which in turn give the estimates for the damping factors and frequencies of (1).

The LS-ESPRIT D Algorithm.
One can easily derive LS-ESPRIT D, as it was done for TLS-ESPRIT D, by appropriate formation of matrix U p for decimation factor D and by solving in a least squares sense for matrix Q.
In this case, Step 3 TLS-ESPRIT D of the Section 3.2 above is replaced by the following.
Step 3 (LS-ESPRIT D). We compute the solution Q of U↑ D Q = U↓ D , in a least squares sense, where U↓ D (U↑ D ) are derived from U p by deleting its top D (bottom D) rows, respectively. Hence, Q = pinv(U↑ D )U↓ D . The eigenvalues λ i of Q give the decimated signal pole estimates, which in turn give the estimates for the damping factors and frequencies of (1).

Computational Considerations for DESE D.
Regarding the computational complexity, DESE D involves one large SVD (singular value decomposition) and one large EVD (eigenvalue decomposition). Note that the pinv operation in Step 4 of the DESE D algorithm is of no computational load, since the SVD of the matrix involved in the pseudoinverse operation is already available from the previous step

Experimental Results
All methods, namely, DESE D, CONDE D, LS-ESPRIT D, and TLS-ESPRIT D have been tested via simulations on a typical two peak reference signal, and two typical 31 P NMR signals, in order to evaluate both robustness and accuracy of parameter estimation in the problem defined by (1). All the experiments have been conducted using the Matlab software.
The first signal is a two-peak signal often used in the literature (reference signal), the exact parameter values of which are presented in Table 1, while the sampling frequency is considered 1.
The second signal is a representative example simulating a typical 31 P NMR signal of perfused rat liver. This 31 P NMR signal comprises a fifth-order model function given in Table 2 by which N data points uniformly sampled at 10 KHz are exactly modeled.
Moreover, the third signal is also a representative example simulating a typical 31 P NMR signal which comprises an eleventh-order model function given in Table 3 by which N data points uniformly sampled at 3 KHz are exactly modeled.
The data points of all signals are perturbed by white Gaussian noise whose real and imaginary components have standard deviation σ v .
Root mean-squared errors of the estimates of all signal parameters are computed using either 500 noise realizations or 3000 noise realizations (excluding failures) for different noise levels. A failure occurs when not all peaks are resolved within specified intervals lying symmetrically around the exact frequencies and when the estimated damping factors are nonpositive.
0.0106 Hz and were deduced from the Cramer-Rao lower bounds of the two peaks at the noise standard deviation where these intervals are touching each other. These values are only used to determine when a failure (bad run) occurs. They depend on the signal parameters and the noise energy and show how far one can go before the two peaks cannot be resolved.
For the five-peak 31 P NMR signal (Figure 2), the halfwidths of the intervals are, respectively, 82, 82, 82, 43, and 82 Hz, the values being derived from the Cramer-Rao lower bounds of the closest peaks 4 and 5 at the noise standard deviation where these intervals are touching each other. The number of complex damped sinusoids to be estimated is set to 5. The Cramer-Rao lower bounds are derived from the exact parameter values and σ v .
For the eleven peak 31 P NMR signal (Figure 3), the halfwidths of the intervals are 8.6, 7.3, 8.6, 3.2, 3.2, 3.4, 3.6, 7.4, 5.5, 2.3, and, 7.7 Hz, for peak 1 to 11, respectively. These values are derived from the Cramer-Rao lower bounds of the closest peaks 4 and 5 at the noise standard deviation where these intervals are touching each other.
Comparative results between all methods are presented below for different noise standard deviations. In Figure 4 failure rates (bad runs) in 500 realizations are depicted as a function of noise standard deviation for the two peak reference signal. In this graphical representation, results are presented for methods DESE2, LS-ESPRIT2, TLS-ESPRIT2, and CONDE2.
Note that DESE2 and DESE1 have fewer bad runs than the other methods under consideration. This is even more evident as the noise increases.
Moreover, we have conducted further experiments to improve the statistical behavior of RMSE figures when the number of bad runs tends to be very large. We decided to increase considerably the number of Monte Carlo trials to   insure that the number of good runs that are taken into account to deduce the RMSE is big enough. In Figure 6 the number of bad runs for 3000 trials are presented for the five peak simulated 31 P NMR signal, for σ v ∈ (0, 2.6], for methods DESE2, DESE3, LS-ESPRIT2, CONDE2, and CONDE3. The same quantity is depicted in Figure 7 for methods DESE2, DESE1, TLS-ESPRIT2, and CONDE1 (=TLS-ESPRIT).
Note that here again DESE2 and DESE3 present fewer bad runs than the other methods, which becomes more evident as the noise increases. In general, DESE3 performs better in terms of robustness than DESE2 as expected, because it brings the peaks even more further apart due to the decimation factor 3 instead of 2, while the matrix size remains practically unchanged.
In Figure 8 root mean square errors of the frequency estimates are presented for DESE2 and TLS-ESPRIT and for σ v ∈ (0, 2.6], for peaks 1 and 4 of the five peak simulated 31 P NMR signal, for 3000 trials. Peak 4 of this signal is considered the most difficult to estimate since it is relatively close to peak 5. The same quantities (also for 3000 trials) for methods DESE2 versus LS-ESPRIT2 are presented in Figure 9, for DESE2 versus CONDE2 in Figure 10, and for methods DESE3 versus CONDE3 in Figure 11.
These graphs show the same trend, according to which, the DESE D technique outperforms the other methods for both peaks 1 and 4, especially for low SNR. In particular, the difference in performance is more evident in the case of  peak 4, which is more difficult to estimate due to its short distance from peak 5. Figures 8 to 11 show that DESE D has always better or at least similar performance compared to the other techniques. More detailed results involving noise standard deviation, number of bad runs, and root meansquared errors of frequency, damping factor, amplitude, and phase estimates for all signals are presented in tabular forms in [27].
It is worth noting that the increased number of trials improves the statistical behavior of the RMSE variance.
The above results suggest in all cases that the DESE D approach performs similarly to the other methods for high S/N ratio. However, for low S/N ratio, despite the similarity of the root mean-squared errors of the estimated parameters, the DESE D technique performs better due to its lower failure rate.
Note that the calculation of the root mean-squared error does not take into account failures and it is normal that methods with small number of bad runs will present larger error than those with big number of bad runs.

Conclusion
In this paper DESE D, a new state-space decimative method, for spectral estimation was presented. It makes use of    decimation by any factor D and SVD, to estimate frequencies, damping factors, amplitudes, and phases of complex damped sinusoids. DESE D makes use of the full data set available and, unlike conventional decimation methods, it imposes no constraints to the size of the Hankel matrix, as decimation increases. It was tested in spectroscopy, one of the most demanding applications of digital signal processing in terms of accuracy. DESE D was compared to a state-of-theart decimative method, along with other state-of-the-art nondecimative ones together with their derived decimative counterparts. Examples on a two-peak reference signal as well as on two typical 31 P NMR signals were presented and it was shown that DESE D performs better than the other methods, especially for low signal-to-noise ratio.