PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Proc Appl Math Mech. Author manuscript; available in PMC 2010 August 16.
Published in final edited form as:
Proc Appl Math Mech. 2009 February 26; 8(1): 10761–10762.
PMCID: PMC2922022
NIHMSID: NIHMS175361

Dependence of accuracy of ESPRIT estimates on signal eigenvalues: the case of a noisy sum of two real exponentials

Abstract

This paper is devoted to estimation of parameters for a noisy sum of two real exponential functions. Singular Spectrum Analysis is used to extract the signal subspace and then the ESPRIT method exploiting signal subspace features is applied to obtain estimates of the desired exponential rates. Dependence of estimation quality on signal eigenvalues is investigated. The special design to test this relation is elaborated.

1 Introduction

In this paper we address the task of estimation of signal parameters having noisy measurements when the signal is the sum of two real exponential functions. For this purpose we apply the method ESPRIT [1] originally aimed at estimation of parameters of complex-valued exponentials (cisoids).

The ESPRIT method exploits the subspace properties of the signal and this is what connects it with Singular Spectrum Analysis (SSA) [2]. SSA is the method of time series analysis aimed at decomposition of the time series into the sum of components, e.g. signal and noise. As an essential step, SSA includes the construction of a subspace which approximates the subspace generated in some way by the signal. ESPRIT uses a basis of this subspace for estimating the rates of the exponential components. After that, the coefficients before the exponentials can be evaluated using the conventional least squares method.

We introduce the signal subspace as follows. For the given series S = (s0, …, sN−1), we select the window length L, where 2 ≤ L ≤ (N + 1)/2. Then S is embedded into RL constructing Sj = (sj−1,…,sj+L−2)T [set membership] RL, j = 1,…, K, K = NL + 1. The linear span of {Sj} is called L-trajectory space of the series S. The L-rank (or simply rank) of the series is defined to be the dimension r of the trajectory space, naturally rL. The rank of the series S is equal to (i) the rank of the trajectory matrix S with the column vectors {Sj}, and (ii) the number of non-zero eigenvalues of the matrix SST, i.e. the number of components of Singular Value Decomposition (SVD) of S. We denote the non-zero eigenvalues as λ1(S) ≥ … ≥ λr(S), and μj(S) = λj(S)/K.

If the signal S has rank r < L, then the trajectory subspace of S is called the signal subspace and determines the shape of the signal up to coefficients before summands. When the series contains both the signal and noise, the signal subspace can be found only approximately (e.g. by means of SSA) that leads to approximate estimation of parameters of the signal.

Algorithms of SSA and ESPRIT are briefly described in Section 2. The design of numerical tests aimed to reveal the relation between properties of ESPRIT estimates and the signal eigenvalues λj(S) is constructed in Section 3. The results on mean square deviations of estimates are presented and discussed in Section 4.

2 Algorithms

Let us consider a noisy signal F= S + [mathematical script N] = (f0, …., fN−1) and the window length selected is L. The first step of both SSA and ESPRIT is the construction of the trajectory matrix F, having the columns Fj = (fj−1, …, fj+L−2)T, j = 1, …, K, K = NL + 1. The second step is the SVD of the matrix F producing, in particular, the eigenvectors Uj of the matrix FFT ordered by decreasing the corresponding eigenvalues λj = λj(F). The next step of SSA is the grouping (summation) of SVD components. For example, to extract the signal S from the given series F, one can select a group of SVD components with indices J, which produces an approximation S of the signal after the diagonal averaging, on condition that the signal is approximately separated from the rest [2]. As a rule, J = {1, …, r}.

Under conditions of approximate separability, the space spanned by {Uj}j[set membership]J approximates the signal subspace. In general, the better is the reconstruction of the signal subspace, the more precise are the estimates of signal parameters.

Let us describe ESPRIT for the signal S = [mathematical script A] + B consisting of the sum of two exponential components [mathematical script A] and B with elements an=caexp(αn)=caranandbn=cbexp(βn)=cbrbn. The aim is to estimate the parameters ra, rb and the coefficients ca, cb. The rank of S equals two and its subspace is approximated by span {U1, U2}. We denote by UΔ the matrix U without its last row, by U[nabla] the matrix U without its first row, and by UΔ+ the Moore-Penrose pseudoinverse of the matrix UΔ. Then the eigenvalues of the matrix UΔ+U of size 2 × 2 are the estimates of the parameters ra and rb of the signal, [r with circumflex]a and [r with circumflex]b [1]. Note that in the noiseless case F = S it holds that [r with circumflex]a = ra and [r with circumflex]b = rb. Having the estimates [r with circumflex]a and [r with circumflex]b, the coefficients ca and cb can be estimated, e.g. by the linear least-squares method.

3 Design of tests

Let us consider the signal S = [mathematical script A] + B, define signal eigenvalues λmin = λ1(S), λmax = λ2(S) and denote their divided by K versions as μmin and μmax. Suppose that we observe F = S + [mathematical script N], where [mathematical script N] is the white gaussian noise with variance σ2. Let S and N be the trajectory matrices of S and [mathematical script N], correspondingly, and ‖ · ‖ is the Frobenius norm which is defined to be the square root of the sum of squares of the entries of the matrix. Hence ‖S2 = λmin + λmax, N2=j=1Lλj(𝒩)LK σ2 and the signal-to-noise ratio ‖S2/‖N2 for the given L is specified by (μmin + μmax)/(Lσ2), where the numerator is actually the weighted sum of squares of the elements of the signal. Below we refer to the Frobenius norm of the trajectory matrix as L-norm of the series. Note that the average of noise eigenvalues λj([mathematical script N]), j = 1, …, L, is approximately equal to Kσ2. To obtain the proper approximation of the two-dimensional signal subspace, we need λmin to be larger than the maximal of noise eigenvalues. Hence, the relation between μmin and σ2 is the key point here.

Obviously, the eigenvalues depend on L. We choose L = N/2 in order to increase μmin and to improve also the separability between signal and noise following the recommendations of SSA theory [2]. Let us fix N = 99 and L = 50. We consider eight test sets of parameters producing exponential components with equal L-norms and the fixed signal L-norms: μ1([mathematical script A]) = μ1(B), μmin + μmax = 2. These test sets differ in the value of μmin, see Table 1 with the test sets ordered by increasing μmin. Note that we consider as the first exponential the one with the largest absolute value of rate: |α| ≥ |β|, i.e. the first exponential component changes not slower than the second one.

Table 1
Eight test sets of the parameters used in simulation and mean square deviations (MSD) for the parameters estimates. For better reading, the significantly large MSD values are highlighted in bold

4 Numerical results

The statistical simulation was performed to estimate the parameters of the model. Two levels of noise were considered, σ2 = 1.e–2 and σ2 = 1.e–4. The signal-to-noise ratio is equal to 4 and 400, correspondingly. The two obtained ESPRIT estimates [r with circumflex]a and [r with circumflex]b were ordered by value to have the same order as ra and rb. Agreement of the results with the signal model was not checked, that is, if the eigenvalue was complex-valued or negative, we formally took its real part.

In our experiments, the noise variance σ2 is relatively low for the signal to be reconstructed with good accuracy. However, based on the previous section, the accuracy of the parameter estimates is assumed to depend not only on σ2 but on μmin also.

The simulation involving 1000 replicates confirmed this guess. Though the signal estimates are precise (not shown in Table 1), the parameters estimates are close to the original values only when μmin is several times larger than σ2, see error values in Table 1. Moreover, it turns out that the estimates of ra and rb are different in the sense that the value of μmin is crucial only for the former one, which corresponds to the faster changing exponential component.

Acknowledgements

This work was supported by NSF/NIGMS BioMath 1-R01-GM072022 grant and GAP award RUB1-1643-ST-06 from the CRDF.

References

1. Roy R, Kailath T. IEEE Trans. Acoust. 1989;37:984–995.
2. Golyandina N, Nekrutkin V, Zhigljavsky A. Analysis of Time Series Structure: SSA and Related Techniques. Chapman&Hall/CRC; 2001.