# Decimative Spectral Estimation with Unconstrained Model Order

Authors: Stavroula-Evita Fotinea; Ioannis Dologlou; Stylianos Bakamidis; Theologos Athanaselis
Journal: Computational and Mathematical Methods in Medicine
 (2012)

---

## Abstract

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.


---

## Body

## 1. 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–8]. The generalized model we use is given by
(1)s(n)=∑i=1p(biej(φ0+φi))e(-di+j2pfi)n=∑i=1pgizin,n=0,…,N-1,
where p is the number of complex damped sinusoids that comprise the measured signal. The objective is to estimate the frequencies fi, damping factors di, amplitudes bi, 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 factorD), 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.Section2.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.


## 2. The DESE_D Method

### 2.1. 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.LetSH be the L×M Hankel signal observation matrix of our deterministic signal of p exponentials s(n),n=0,…,N-1(2)SH=[s̃0s̃1⋯s̃M-1]
with L-D≤M,p<L-D and L+M-1=N. Note that s̃n  are the column vectors of SH, for n=0,1,…,M-1. In particular s̃nT=[s(n),s(n+1),…,s(n+L-1)], with s(n) being the nth sample of the signal s. Let the L-D×M matrices S↓D and S↑D be the D order lower shift (top D rows deleted) and the D order upper shift (bottom D rows deleted) equivalents of SH. The best choice for L and M is discussed in Section 2.2.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,
(3)XS↑D=S↓D
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↓Dpinv(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 becauseSH 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 thatp of the eigenvalues of the matrix
(4)X=S↓Dpinv(S↓D)
are equal to the decimated signal pole estimates ziD,i=1,…,p.
Let us consider the well-known Vandermonde decomposition ofSH:(5)SH=AGBT,
where superscript T denotes transpose and the L×p matrix A, p×p matrix G and M×p matrix B are defined as follows:
(6)A=(αH(z1)⋯αH(zp)),aH(z)=(z0z1⋯zL-1)T,B=(bH(z1)⋯bH(zp)),bH(z)=(z0z1⋯zM-1)T,G=diag⁡(g1,…,gp).
We can then easily write S↓D=A↓DGBT and S↑D=A↑DGBT where A↓D and A↑D are defined from A, similarly to the way S↓D and S↑D are defined from SH.
Hence, (3) can be written as
(7)(XA↑D-A↓D)GBT=0⟺XA↑D=A↓D.
The latter system of linear equations has an infinite number of solutions given by
(8)X=A↓D(AH↑DA↑D)-1AH↑D+ΔH,
where superscript H denotes Hermitian conjugate and ΔHA↑D=0. Let us consider now the p×p diagonal matrix ΦD containing the decimated signal poles ziD,i=1,…,p. It is easy then to see that A↓D=A↑DΦD. Moreover, A↑D=A↓DΦD-1 resulting in ΔHA↓D=0.
Hence, since∥X∥F2=∥X0∥F2+∥ΔH∥F2 (∥·∥F denotes the Frobenius norm), the minimum-norm solution to (7) which we compute is
(9)X0=A↓D(AH↑DA↑D)-1AH↑D.
Furthermore, knowing that A↓D=A↑DΦD it is also valid that AH↑DA↓D=AH↑DA↑DΦD which leads to (AH↑DA↑D)-1AH↑DA↓D=ΦD.
By settingW=A↓D and Y=(AH↑DA↑D)-1AH↑D it is straightforward to see that X0=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 X0(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 forD=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 matrixSH 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 ofS↑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↑De has rank p. Then X is computed from XS↑De≈S↓D which gives rise to an overdetermined system of equations with the following solution X=S↓Dpinv(S↑De).
Note that since matrixS↑De has rank p, X is also of rank p (minimum of ranks of S↑De 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 fi (frequency estimates) and di (damping factor estimates) and a computation in a total least squares sense of estimates gi then takes place. In this way, complex-valued linear parameter estimates of gi are calculated, from which amplitude bi and phases φ0+φi estimates are determined as the magnitudes and angles of gi, respectively.


### 2.2. DESE_D Algorithmic Presentation

LetSH be the L×M Hankel signal observation matrix of our deterministic signal of p exponentials s(n), n=0,…,N-1, with L-D≤M,p<L-D and L+M-1=N, where D denotes the decimation factor.The proposed algorithm, involves the following five steps.Step 1 (DESE_D).
We compute theL×M matrix SH of (2) from the N data points s(n) of (1).Step 2 (DESE_D).
We compute theS↓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 SH.
The best results are obtained when we use the(L-D)×M matrices S↓D and S↑D as square as possible [11, 19–21].Step 3 (DESE_D).
We compute the enhanced versionS↑De of S↑D in the following way. We employ the SVD of S↑D, S↑D=U↑DΣ↑DV↑DH and we truncate to order p by retaining only the largest p singular values.Step 4 (DESE_D).
We compute matrixX=S↓Dpinv(S↑De).
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 zi replaced by the estimates and s(n) given by the signal data points.
MatrixX 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.


### 2.3. 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 upper shift (bottom row deleted) of the original Hankel SH of (2) with L-1≤M,p<L-1 and L+M-1=N.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 versionS↓1e of S↓1 in the following way. We employ the SVD of S↓1, S↓1=U↓1Σ↓1V↓1Hand we truncate to order p by retaining only the largest p singular values.
In the same way, we compute the enhanced versionS↑1e of S↑1.Step 4 (a).
We compute matrixX=(S↓1e)pinv(S↑1e).
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 matrixS↑1), the nondecimative method (D=1, for DESE_D) is identical to a method proposed in [24], the MATPEN method.
For MATPEN(=DESE_1), Step4(a) is replaced by X=S↓1pinv(S↑1e).
The eigenvaluesλ̂i of X give the signal pole estimates, which in turn give the estimates for the damping factors and frequencies of (1).


### 2.4. 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 sizen=(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 orderp and proceeds with estimation of the frequency and damping factor in a total least squares sense.The method’s algorithmic presentation for decimation factorD involves the following five steps.Step 1 (CONDE_D).
We compute theL×M(L=M=N/(2D)) Hankel matrix Ck that corresponds to the kth decimated signal, ck(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 matrixC by concatenating Ck,k=1,…,D as shown below:
(10)C=[C1:C2:⋯:CD].
We then compute a global covariance matrix R=CCH.Step 3 (CONDE_D).
We compute the eigen analysis ofR=UΛUH  to deduce U, which in turn is truncated to order p, thus, yielding Up.Step 4 (CONDE_D).
We compute the solutionQ of U↑pQ=U↓p, in a total least squares sense, where U↓p(U↑p) are derived from Up 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 zi replaced by the estimates and s(n) given by the signal data points.


## 2.1. 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.LetSH be the L×M Hankel signal observation matrix of our deterministic signal of p exponentials s(n),n=0,…,N-1(2)SH=[s̃0s̃1⋯s̃M-1]
with L-D≤M,p<L-D and L+M-1=N. Note that s̃n  are the column vectors of SH, for n=0,1,…,M-1. In particular s̃nT=[s(n),s(n+1),…,s(n+L-1)], with s(n) being the nth sample of the signal s. Let the L-D×M matrices S↓D and S↑D be the D order lower shift (top D rows deleted) and the D order upper shift (bottom D rows deleted) equivalents of SH. The best choice for L and M is discussed in Section 2.2.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,
(3)XS↑D=S↓D
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↓Dpinv(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 becauseSH 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 thatp of the eigenvalues of the matrix
(4)X=S↓Dpinv(S↓D)
are equal to the decimated signal pole estimates ziD,i=1,…,p.
Let us consider the well-known Vandermonde decomposition ofSH:(5)SH=AGBT,
where superscript T denotes transpose and the L×p matrix A, p×p matrix G and M×p matrix B are defined as follows:
(6)A=(αH(z1)⋯αH(zp)),aH(z)=(z0z1⋯zL-1)T,B=(bH(z1)⋯bH(zp)),bH(z)=(z0z1⋯zM-1)T,G=diag⁡(g1,…,gp).
We can then easily write S↓D=A↓DGBT and S↑D=A↑DGBT where A↓D and A↑D are defined from A, similarly to the way S↓D and S↑D are defined from SH.
Hence, (3) can be written as
(7)(XA↑D-A↓D)GBT=0⟺XA↑D=A↓D.
The latter system of linear equations has an infinite number of solutions given by
(8)X=A↓D(AH↑DA↑D)-1AH↑D+ΔH,
where superscript H denotes Hermitian conjugate and ΔHA↑D=0. Let us consider now the p×p diagonal matrix ΦD containing the decimated signal poles ziD,i=1,…,p. It is easy then to see that A↓D=A↑DΦD. Moreover, A↑D=A↓DΦD-1 resulting in ΔHA↓D=0.
Hence, since∥X∥F2=∥X0∥F2+∥ΔH∥F2 (∥·∥F denotes the Frobenius norm), the minimum-norm solution to (7) which we compute is
(9)X0=A↓D(AH↑DA↑D)-1AH↑D.
Furthermore, knowing that A↓D=A↑DΦD it is also valid that AH↑DA↓D=AH↑DA↑DΦD which leads to (AH↑DA↑D)-1AH↑DA↓D=ΦD.
By settingW=A↓D and Y=(AH↑DA↑D)-1AH↑D it is straightforward to see that X0=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 X0(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 forD=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 matrixSH 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 ofS↑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↑De has rank p. Then X is computed from XS↑De≈S↓D which gives rise to an overdetermined system of equations with the following solution X=S↓Dpinv(S↑De).
Note that since matrixS↑De has rank p, X is also of rank p (minimum of ranks of S↑De 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 fi (frequency estimates) and di (damping factor estimates) and a computation in a total least squares sense of estimates gi then takes place. In this way, complex-valued linear parameter estimates of gi are calculated, from which amplitude bi and phases φ0+φi estimates are determined as the magnitudes and angles of gi, respectively.


## 2.2. DESE_D Algorithmic Presentation

LetSH be the L×M Hankel signal observation matrix of our deterministic signal of p exponentials s(n), n=0,…,N-1, with L-D≤M,p<L-D and L+M-1=N, where D denotes the decimation factor.The proposed algorithm, involves the following five steps.Step 1 (DESE_D).
We compute theL×M matrix SH of (2) from the N data points s(n) of (1).Step 2 (DESE_D).
We compute theS↓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 SH.
The best results are obtained when we use the(L-D)×M matrices S↓D and S↑D as square as possible [11, 19–21].Step 3 (DESE_D).
We compute the enhanced versionS↑De of S↑D in the following way. We employ the SVD of S↑D, S↑D=U↑DΣ↑DV↑DH and we truncate to order p by retaining only the largest p singular values.Step 4 (DESE_D).
We compute matrixX=S↓Dpinv(S↑De).
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 zi replaced by the estimates and s(n) given by the signal data points.
MatrixX 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.


## 2.3. 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 upper shift (bottom row deleted) of the original Hankel SH of (2) with L-1≤M,p<L-1 and L+M-1=N.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 versionS↓1e of S↓1 in the following way. We employ the SVD of S↓1, S↓1=U↓1Σ↓1V↓1Hand we truncate to order p by retaining only the largest p singular values.
In the same way, we compute the enhanced versionS↑1e of S↑1.Step 4 (a).
We compute matrixX=(S↓1e)pinv(S↑1e).
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 matrixS↑1), the nondecimative method (D=1, for DESE_D) is identical to a method proposed in [24], the MATPEN method.
For MATPEN(=DESE_1), Step4(a) is replaced by X=S↓1pinv(S↑1e).
The eigenvaluesλ̂i of X give the signal pole estimates, which in turn give the estimates for the damping factors and frequencies of (1).


## 2.4. 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 sizen=(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 orderp and proceeds with estimation of the frequency and damping factor in a total least squares sense.The method’s algorithmic presentation for decimation factorD involves the following five steps.Step 1 (CONDE_D).
We compute theL×M(L=M=N/(2D)) Hankel matrix Ck that corresponds to the kth decimated signal, ck(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 matrixC by concatenating Ck,k=1,…,D as shown below:
(10)C=[C1:C2:⋯:CD].
We then compute a global covariance matrix R=CCH.Step 3 (CONDE_D).
We compute the eigen analysis ofR=UΛUH  to deduce U, which in turn is truncated to order p, thus, yielding Up.Step 4 (CONDE_D).
We compute the solutionQ of U↑pQ=U↓p, in a total least squares sense, where U↓p(U↑p) are derived from Up 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 zi replaced by the estimates and s(n) given by the signal data points.


## 3. 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 factorD, 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.


### 3.1. 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 di and the frequencies fi 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 theL×M Hankel matrix SH of (2) from the N data points s(n) of (1):
(11)SH=UL×LΣL×MVM×MH,
where L≤M. The best results are obtained when we use L=M(+1)=N/2.Step 2 (TLS-ESPRIT).
We truncateU,Σ,V to order p and compute Sp=UpΣpVpH where Up,Σp,VPare the first p columns of UL×L,ΣL×M,VM×MH.Step 3 (TLS-ESPRIT).
We compute the solutionQ of U↑pQ=U↓p, in a total least squares sense, where U↓p(U↑p) are derived from Up 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 zi 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.


### 3.2. 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 factorD, the decimative version of TLS-ESPRIT can be easily derived. More precisely, appropriate formation of matrix Up—similar to that of matrix SH 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 theL×M Hankel matrix SH of (2) from the N data points s(n) of (1):
(12)SH=UL×LΣL×MVM×MH,
where L≤M. The best results are obtained when we use L=M(+1)=N/2.Step 2 (TLS-ESPRIT).
We truncateU,Σ,V to order p and compute: Sp=UpΣpVpH where Up,Σp,Vp are the first p columns of UL×L,ΣL×M,VM×MH.Step 3 (TLS-ESPRIT_D).
We compute the solutionQ of U↑DQ=U↓D, in a total least squares sense, where U↓D(U↑D) are derived from UP 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 zi 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.


### 3.3. The LS-ESPRIT Algorithm

If instead of computing in a total least squares sense the solutionQ of U↑pQ=U↓p, one employs the least squares solution, one uses LS-ESPRIT [12].In this case, Step3_TLS-ESPRIT of the Section 3.1 above is replaced by the following.Step 3 (LS-ESPRIT).
We compute the solutionQ of U↑pQ=U↓p, in a least squares sense, where U↓p(U↑p) are derived from Up 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).


### 3.4. The LS-ESPRIT_D Algorithm

One can easily derive LS-ESPRIT_D, as it was done for TLS-ESPRIT_D, by appropriate formation of matrixUp for decimation factor D and by solving in a least squares sense for matrix Q.In this case, Step3_TLS-ESPRIT_D of the Section 3.2 above is replaced by the following.Step 3 (LS-ESPRIT_D).
We compute the solutionQ of U↑DQ=U↓D, in a least squares sense, where U↓D(U↑D) are derived from Up 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).


### 3.5. 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 thepinv 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 (Step 3_DESE_D). Moreover, only p of the eigenvalues of matrix X involved in the EVD, are nonzero, which can considerably reduce the complexity of the large EVD in DESE_D if a fast algorithm is used. On the other hand, the TLS-ESPRIT_D/LS-ESPRIT_D decimative versions require one large SVD and one small EVD. Consequently the difference of the complexity of DESE_D versus TLS-ESPRIT_D/LS-ESPRIT_D cannot be considered as a drawback for DESE_D. Regarding the complexity of CONDE_D, the fact that it uses lower dimension Hankel matrices (as decimation increases) improves its computational characteristics with respect to the other methods.


## 3.1. 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 di and the frequencies fi 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 theL×M Hankel matrix SH of (2) from the N data points s(n) of (1):
(11)SH=UL×LΣL×MVM×MH,
where L≤M. The best results are obtained when we use L=M(+1)=N/2.Step 2 (TLS-ESPRIT).
We truncateU,Σ,V to order p and compute Sp=UpΣpVpH where Up,Σp,VPare the first p columns of UL×L,ΣL×M,VM×MH.Step 3 (TLS-ESPRIT).
We compute the solutionQ of U↑pQ=U↓p, in a total least squares sense, where U↓p(U↑p) are derived from Up 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 zi 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.


## 3.2. 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 factorD, the decimative version of TLS-ESPRIT can be easily derived. More precisely, appropriate formation of matrix Up—similar to that of matrix SH 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 theL×M Hankel matrix SH of (2) from the N data points s(n) of (1):
(12)SH=UL×LΣL×MVM×MH,
where L≤M. The best results are obtained when we use L=M(+1)=N/2.Step 2 (TLS-ESPRIT).
We truncateU,Σ,V to order p and compute: Sp=UpΣpVpH where Up,Σp,Vp are the first p columns of UL×L,ΣL×M,VM×MH.Step 3 (TLS-ESPRIT_D).
We compute the solutionQ of U↑DQ=U↓D, in a total least squares sense, where U↓D(U↑D) are derived from UP 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 zi 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.


## 3.3. The LS-ESPRIT Algorithm

If instead of computing in a total least squares sense the solutionQ of U↑pQ=U↓p, one employs the least squares solution, one uses LS-ESPRIT [12].In this case, Step3_TLS-ESPRIT of the Section 3.1 above is replaced by the following.Step 3 (LS-ESPRIT).
We compute the solutionQ of U↑pQ=U↓p, in a least squares sense, where U↓p(U↑p) are derived from Up 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).


## 3.4. The LS-ESPRIT_D Algorithm

One can easily derive LS-ESPRIT_D, as it was done for TLS-ESPRIT_D, by appropriate formation of matrixUp for decimation factor D and by solving in a least squares sense for matrix Q.In this case, Step3_TLS-ESPRIT_D of the Section 3.2 above is replaced by the following.Step 3 (LS-ESPRIT_D).
We compute the solutionQ of U↑DQ=U↓D, in a least squares sense, where U↓D(U↑D) are derived from Up 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).


## 3.5. 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 thepinv 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 (Step 3_DESE_D). Moreover, only p of the eigenvalues of matrix X involved in the EVD, are nonzero, which can considerably reduce the complexity of the large EVD in DESE_D if a fast algorithm is used. On the other hand, the TLS-ESPRIT_D/LS-ESPRIT_D decimative versions require one large SVD and one small EVD. Consequently the difference of the complexity of DESE_D versus TLS-ESPRIT_D/LS-ESPRIT_D cannot be considered as a drawback for DESE_D. Regarding the complexity of CONDE_D, the fact that it uses lower dimension Hankel matrices (as decimation increases) improves its computational characteristics with respect to the other methods.


## 4. 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 typical31P 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 Table1, while the sampling frequency is considered 1.Table 1
Exact parameter values of the standard two peak reference signal withφi=0.
Peakifi (Hz)di (rad/s)biψi(a)10.20.011020.220.0210(a)ψi=φ0180π expresses the phase in degrees.The second signal is a representative example simulating a typical31P NMR signal of perfused rat liver. This 31P 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.Table 2
Exact parameter values of the five peak simulated31P NMR signal, modelled by (1) with φi=0.
Peakifi (Hz)di (rad/s)biψi(a)1−13792086.1152−6852569.9153−2711976.01543531172.815547880817.015(a)ψi=φ0180π expresses the phase in degrees.Moreover, the third signal is also a representative example simulating a typical31P 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.Table 3
Exact parameter values of the eleven peak simulated31P NMR signal, modelled by (1) with φi=0.
Peakifi (Hz)di (rad/s)biψi(a)1−8650751352−70501501353−5450751354152501501355168501501356292501501357308501501358360251501359440285.7140013510490256013511530200500135(a)ψi=φ0180π expresses the phase in degrees.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 (excludingfailures) for different noise levels. A failure occurs when not all peaksare resolved within specified intervals lying symmetrically around the exact frequencies and when the estimated damping factors are nonpositive.For the two peak reference signal (Figure1), the half-widths of the intervals are, respectively, 0.0094 Hz and 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.Figure 1
Fast Fourier transform (magnitude) of the two peak signal.For the five-peak31P NMR signal (Figure 2), the half-widths 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.Figure 2
Fast Fourier transform (magnitude) of the five peak simulated31P NMR signal.For the eleven peak31P NMR signal (Figure 3), the half-widths 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.Figure 3
Fast Fourier transform (magnitude) of the eleven peak simulated31P NMR signal.Comparative results between all methods are presented below for different noise standard deviations.For all methods, we have usedN=128 and M=N/2=64 except for CONDE_D, for which M=32(21) for D=2(3), respectively.In Figure4 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.Number of Bad runs in 500 realisations forσv∈(0,2.6] for the two-peak reference signal for methods DESE2, LS-ESPRIT2, TLS-ESPRIT2 and CONDE2.
(a)(b)The same quantity, namely, number of bad runs in 500 realizations, is depicted in Figure5 for the eleven peak simulated 31P NMR signal, for methods DESE2 & 1, CONDE2 & 1, LS-ESPRIT2 & 1, and TLS-ESPRIT2.Number of Bad runs in 500 realisations forσv∈(0,2.6] for the eleven peak simulated 31P NMR signal for methods DESE2, DESE1, CONDE2, CONDE1, LS-ESPRIT2, LS-ESPRIT1, and TLS-ESPRIT2.
(a)(b)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 Figure6 the number of bad runs for 3000 trials are presented for the five peak simulated 31P 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).Number of bad runs in 3000 realisations forσv∈(0,2.6]  for the five peak simulated 31P NMR signal for methods DESE2, DESE3, LS-ESPRIT2, CONDE2 and CONDE3.
(a)(b)Number of bad runs in 3000 realisations forσv∈(0,2.6] for the five peak simulated 31P NMR signal for methods DESE2, DESE1, TLS-ESPRIT2 and CONDE1 (=TLS-ESPRIT).
(a)(b)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 Figure8 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 31P 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.Figure 8
Frequency root mean square errors in 3000 realisations for Peaks 1 and 4 of the five peak simulated31P NMR signal, for methods DESE2 & TLS-ESPRIT and for σv∈(0,2.6].The same quantities (also for 3000 trials) for methods DESE2 versus LS-ESPRIT2 are presented in Figure9, for DESE2 versus CONDE2 in Figure 10, and for methods DESE3 versus CONDE3 in Figure 11.Figure 9
Frequency root mean square errors in 3000 realisations for Peaks 1 and 4 of the five peak simulated31P NMR signal, for methods DESE2 & LS-ESPRIT2 and for σv∈(0,2.6].Figure 10
Frequency root mean square errors in 3000 realisations for Peaks 1 and 4 of the five peak simulated31P NMR signal, for methods DESE2 & CONDE2 and for σv∈(0,2.6].Figure 11
Frequency root mean square errors in 3000 realisations for Peaks 1 and 4 of the five peak simulated31P NMR signal, for methods DESE3 & CONDE3 and for σv∈(0,2.6].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. Figures8 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 mean-squared 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 highS/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.


## 5. Conclusion

In this paper DESE_D, a new state-space decimative method, for spectral estimation was presented. It makes use of decimation by any factorD 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-the-art 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 P31 NMR signals were presented and it was shown that DESE_D performs better than the other methods, especially for low signal-to-noise ratio.



---

*Source: 917695-2012-02-22.xml*