Enhanced resolution based on minimum variance estimation and exponential data modeling

Signal Processing 33 (1993) 333-355 Elsevier

333

Enhanced resolution based on minimum variance estimation and exponential data modeling Sabine Van Huffel* ESAT Laboratory, Department of Electrical Engineering, Katholieke Universiteit Leuven, Kardinaal Mercierlaan 94, 3001 Heverlee, Belgium Received 7 May 1992 Revised 3 December 1992

Abstract. In various signal processing applications, it is desired to appropriately modify a given data set so that the modified data set possesses prescribed properties. The modification of the given data set serves as a preprocessing step of 'cleaning up' the data before estimating the values of the signal parameters. In this paper, evaluation and improvement of a signal enhancement algorithm, originally proposed by Tufts, Kumaresan and Kirsteins and recently generalized by Cadzow, are presented. In essence, the newly proposed algorithm first arranges the data in a very rectangular (instead of a square) Hankel structured matrix in order to make the corresponding signal-only data matrix orthogonal to the noise, then computes a minimum variance (instead of a least squares) estimate of the signalonly data matrix and finally restores the Hankel structure of the computed minimum variance estimate. An extensive set of simulations is given demonstrating a significant improvement in resolution performance over Cadzow's method at a comparable parameter accuracy. Moreover, arranging the data in a very rectangular matrix reduces drastically the required computation time. In particular, the newly proposed signal enhancement algorithm can be successfully applied to the quantitative time-domain analysis of Nuclear Magnetic Resonance (NMR) data. Zusammenfassnng. In verschiedenen Signalverarbeitungs-Anwendungen ist es erwiinscht, einen gegebenen Datensatz passend zu modifizieren, so dab der ge~inderte Datensatz vorgeschriebene Eigenschaften besitzt. Die Datenmodifikation dient als Vorverarbeitungsschritt zur 'Datenbereinigung' vor der Sch~itzung der Werte von Signalparametern. In diesem Beitrag werden die Auswertung und Verbesserung eines Signalaufbereitungs-Verfahrens vorgesteUt, das urspriinglich von Tufts, Kumaresan und Kirsteins vorgeschlagen und ktirzlich von Cadzow verallgemeinert wurde. Im wesentlichen ordnet der neu vorgeschlagene Algorithmus zuerst die Daten in einer stark rechteckftrmigen ( start quadratischen ) Hankelmatrix an, um die entsprechende reine Signaldaten-Matrix orthogonal zum Rauschen zu machen; dann berechnet er eine Schiitzung der Signaldaten-Matrix orthogonal zum Rauschen zu machen; dann berechnet er eine Sch~itzung der Signaldaten-Matrix mit miminaler Varianz (statt mit kleinstem Quadrat), und schliel31ich stellt er die Hankelstruktur der berechneten Minimum-Varianz-Schiitzung wieder her. Ein umfangreicher Satz von Simulationen wird angegeben. Er demonstriert eine merkliche Verbesserung beziiglich der Aufltsung gegentther Cadzows Methode bei vergleichbarer Parameter-Genauigkeit. Dariiberhinans verringert das Anordnen der Daten in einer stark rechteckigen Matrix die effordediche Rechenzeit drastisch. Insbesondere l~13tsich der ueu vorgeschlagene Signalverbesserungs-Algorithmus erfolgreich auf die quantitative Zeitbereichs-Analyse von Kerspinresonanz(NMR-) Daten anwenden. R t s u m t . Dans de nombreuses applications de traitement du signal, on dtsire modifier de mani~re approprite un ensemble de donntes de telle sorte que l'ensemble de donntes modifite poss~de des proprittts prescrites. La modification de l'ensemble de donntes sert d'ttape de prt-traitement an 'nettoyage' des donntes avant estimation des valeurs des param~tres du signal. Nous prtsentons dans cet article une 6valuation et une amtlioration d'un algorithme de rehaussement du signal origineilement propos6 par Tufts, Kumaresan et Kirsteins et rtcemment gtntralis6 par Cadzow. En essence, le nouvel algorithme propos6 arrange tout d'abord les donntes dans une matrice structurte de Hankel tr~s rectangulaire (au lieu d'Stre carrte) de faqon ~ rendre la matrice de donntes 'signal seulement' orthogonale au bruit, puis calcule une estimte au minimum de variance (an lieu d'&re anx moindres carrts) de la matrice de donntes 'signal seulement' et finalement restaure la structure de Hankel de l'estimte au minimum de variance. Nous donnons un ensemble 6tendu

Correspondenc to: Dr. Sabine van Huffel, ESAT Laboratory, Department of Electrical Engineering, Katholieke Universiteit Leuven, Kardinaal Mercierlaan 94, 3001 Heverlee, Belgium. Email: [email protected] * Research Associate of the Belgian NFWO (National Fund for Scientific Research).

0165-1684/93/$06.00 © 1993 Elsevier Science Publishers B.V. All rights reserved

334

S. Van Huffel / Enhanced resolution and exponentials

de simulations qui montre une amtlioration significativedes performances de rtsolution par rapport ~tla mtthode de Cadzow pour une prtcision des param~trescomparable.De plus, l'arrangement des donntes dans une matrice tr~srectangulairertduit de mani~redrastique le temps de calcul requis. En particulier, l'algorithme de rehaussement de signal nouvellementpropos6 peut 8tre appliqu6 avec succ~s ~t l'analyse quantitative temporelle de donntes de rtsonance magnttique nucltalre (NMR). Keywords. Signal enhancement; singular value decomposition; total least squares; minimum variance; Monte-Carlo simulation.

1. Introduction

x=.~+Xw,

In 1982, Tufts et al. [ 16] presented a method, based on the singular value decomposition, for estimating the signal component of an observed data vector. Recently, Cadzow [2] has generalized this method and applied to the enhancement of sinusoidal and exponential signals. In this paper, we show how to improve the performance of these matrix-based methods. Little prior information is required and the proposed method can be applied as a preprocessing step to a large class of subspace-based signal estimation methods, e.g. in sinusoidal and exponential data modeling and in the automatic removal of annoying low-level interfering components for enhancement of speech signals and NMR data. In essence the data are arranged in a matrix which is then modified in order to possess prescribed properties, i.e. a specific rank determined by the model order and a Hankel structure. The resulting modified data set can be considered as a 'cleaned up' version of the original data whereby corrupting noise, measurement distortion or theoretical mismatch present in the original data set are removed. A closely related method, based on a different approach, has very recently been presented by Tufts and Shah [17]. In that paper, matrix-perturbation results of Vaccaro, Li, Tufts and Kirsteins [ 13, 14, 10] are used to improve the existing signal enhancement algorithms provided some prior knowledge or estimates of the power spectrum density or correlation structure of the perturbing noise are available.

where .~ contains the exact signal component and Xw represents the noise. We are interested in estimating the signal component )( from the observed data vector X. Let us assume that the exact signal satisfies a model of order K. For example, in exponential data modeling this assumption implies that the N data points $n of satisfy the following time-domain model function:

2. Algorithms Consider a vector X = Ix0 . . . . . xN- ~] T of N observations which can be represented as SignalProcessing

(1)

K

$n = Y~ CkZ'k" k~l K

= ~_, (akeJg.)e(--dk+j2~rfk)tn, k=l

n=0 ..... N-I.

(2)

where j = ~ / - 1 and t. is the time lapse between the effective time origin and sample x.. The objective is to estimate the frequenciesfk, damping factors dk, amplitudes ak and phases ~bk,k = 1..... K. The newly proposed signal enhancement algorithm can be best understood by describing a general signal enhancement algorithm, called P-estim (stands for 'signal parameter estimation with preprocessed data' ). that includes Cadzow' s method and the newly proposed signal enhancement algorithm as its special cases.

2.1. The general signal enhancement algorithm Pestim Data preprocessing This preprocessing procedure consists of the following three steps.

Step 1. Form a Hankel data matrix HL × M, L >/M > K and N = L + M - 1, from the observed data vector X such that

335

S. Van Huffel/ Enhanced resolution and exponentials

tively by starting from the cleaned-up

data vector 2

obtained in the previous iteration step.

XL,-_l

x,

Parameter estimation

.‘.

Finally, the signal parameters are extracted from the Srep 2. Modify the Hankel matrix H as follows superscript H denotes the conjugated

(the

(SVD)

[81: -

(4)

ULXL~LLXMV:xM,

whereZ=diag(a,, . . ..~~).a~&... aa,. ( b) Correct the singular values by applying a correction functionf,,, cffolr’s purpose will be discussed in Section 2.2) and truncate to rank K: HK = ~,fco,A~l)

(5)

VP,

where

fJ=[U*;

v,

K

L-K

V=[V,;

I,

data modeling

algorithm. For the 4K signal

fkrdk, ak and +k in (2) can be estimated

through linear prediction

(see e.g. [ 11, 3,

techniques

5, 201). A computationally more efficient and more precise alternative is Kung et al.‘s method [ 121, known as HSVD in NMR [ 4,5], which circumvents polynomial rooting and root selection by representing the signal in a state space model setting. An improved variant, based on total least squares (TLS) and called HTLS, is presented in [ 201. In summary, these algorithms are as follows. Step 1. Compute the singular value decomposition (SVD) of the f. X fi Hankel structured data matrix fi, i.e.

v, 1

K

signal estimation

example, in exponential parameters

H LxM

data vector 8 by means of an appro-

priate subspace-based

transpose) :

(a) Compute the singular value decomposition

final cleaned-up

M-K

and

K

(7)

M-K

Step 3. Compute each element 9,, of the cleaned-up data vector 2, which estimates the signal component 2, from the rank-reduced

data matrix HK by arithmetic

aging alongs its antidiagonals,

n=O,

. . . . N-l,

aver-

i.e.

(6)

where hlj denotes the (fj) the element of HKI cu=max(l,n-L+2) andp=min(M,n+l).Thisis equivalent to making a Hankel matrix approximation to HK prior to extracting the estimated signal parameters. If desired, these three steps can be repeated itera-

and truncate A to a rank K matrix fiK:

OK,eK, pKare the first K columns of 0, 2, p respectively. In order to obtain the best parameter accuracy fi (and fiK) is as square as possible, i.e. t- fi (although the minimum is very flat). Step 2. Compute the solution & of the (incompatible) set * ?.

(8)

&Q = UK,

where aK (respectively & are derived from t?K (or I!?~$.~‘* if a balanced splitting is used) by omitting its first (respectively last) row. Solving (8) with least squares (LS) results in Kung et al.‘s algorithm [ 121, called HSVD here. Computing Vol. 33, No. 3,

September 1993

S. Van Huffel / Enhanced resolution and exponentials

336

instead the TLS solution ~ of (8) results in an improved variant [20], called HTLS here. Once Q is estimated, its eigenvalues A~ give the signal pole estimates, i.e. Ak =Z~ =e-d~+J2~/k,

k = 1. . . . . K,

yielding the desired frequency and damping factor estimated dk andfk. Step 3. Finally, compute the LS solution (~= [ ~ . . . . . dK] T of the set

obtained by fitting the N model equations (2) to the data points .~ (where the zk are replaced by the computed estimates ~k). ~ yields the complex-valued linear parameter estimates ~k----dkej'k which contain the amplitude and phase estimated t~kand ~k. Note that this signal enhancement algorithm (including all special cases listed below) can be straightforwardly adapted to differently structured da m matrices, e.g. a Toeplitz-Hankel matrix as occurring in sinusoidal modeling [3], or a Toeplitz-Toeplitz matrix as occurring in transfer function modeling [ 18] (also called recursive modeling in [3] ). 2.2. Special cases

The algorithm P-estim reduced to Cadzow's method [ 2, 3 ], originally proposed by Tufts et al. in a noniterative form [ 16], when and

L=£,

M=M,

(9)

i.e. the singular values are not corrected and the Hankel matrix used in the preprocessing procedure has the same dimensions as the Hankel matrix used in the chosen signal parameter estimation method. For example, if HSVD or HTLS are used a square Hankel matrix H is considered in the preprocessing procedure. We call the method a least squares (LS) estimation method (following the terminology of [7] ), if f¢o~(,~)=,~ i.e.

but

L>>M>K

L4=/~, M4=/14),

since we seek to approximate H in (3) by a matrix of SignalProcessing

L>> M > K and f¢or~('~) = (_y2_ Lo.21) 1/2 with tr 2 >Lcr 2,

ANxKC~ XNx 1,

fco=(,~)=,~

rank K in least squares sense. Of course, if L =/~ and M = i¢/the LS estimation method coincides with Cadzow's method defined by (9). The condition L > > M ensures that the signal-only subspace Range(/~) given H = / ) + W will approximately be orthogonal to the noise subspace Range(W) in the sense that ~I~w--0 so that both spaces can be separated. Consider now the case where

(10)

where o'~ is an estimate of the variance of the (possibly complex) noise added to the signal. This is a wellknown correction function [ 12]. As before, the condition L>>M ensures that the signal subspace will approximately be orthogonal to the noise subspace but, additionally, the correction function applied to the singular values removes the noise present in the data provided the noise is white. This method is called here the singular value adaptation (SVA) method. Next, we consider the minimum variance (MV) estimation method which requires L>>M>K with

and fco~r('~)-- ("~2--L°'v21)~ -1

2 trr2 > Lo%.

(11)

This method computes the MV estimate of the signal-only data matrix /~ given H = I = I + W provided W H W = lzlm with/z an unknown scalar (IM is the identity matrix of order M) and the signal-only data are 'orthogonal' to the noise in the sense that/~HW= 0 [7]. In practice, these conditions are never satisfied exactly, but it can be shown (see Section 3) that these conditions are more and more satisfied with increasing overdetermination provided the noise is white and has bounded fourth moments. That is the reason why L >> M is required. As for the SVA method, or2 is the estimated noise variance and can be known exactly or computed either from the last data points provided they are pure noise (tail of the signal) or from the noise singular values contained in ~2 in (5), i.e. o"v-2( ( M - K ) L ) - IE/M=x+1tr2- Note that try changes if we iterate the method. Hence, in the higher iterations o-vis no longer known but can only be computed as explained above.

S. Van Huffel / Enhanced resolution and exponentials

2.3. Computational considerations and extensions Based on the number of operations involved (see e.g. [ 8, p. 239] ), it is easy to see that the computationally most intensive part of the algorithm, namely the SVD (4) and the computation of (5) of the truncated matrix HK (second step of the preprocessing procedure), requires significantly less work whenever L >> M. Indeed, define by one flop the amount of work associated with one floating point operation (one addition, multiplication, division . . . . ). If only third order terms are considered, the number of flops involved in preprocessing a real L × M rank K matrix H of the form (3) is given by

6LM2+2OM3+2MKL

ifL>~5M/3,

14LM2+8M3+2MKL

ifL<5M/3,

22M3+2M2K

if L = M .

Denote by Lc (Me) the number of rows (columns) as used in Cadzow's method, and by Lmv (Mmv) the number of rows (columns) as used in the MV estimation method. We must have N = Me + L~ - 1 = Mmv + Lmv- 1, i.e. the number of data points is the same. If/~ in (7) is approximately square (/~ = ~ ) then also L~-~Mc. This means that the MV estimation (LS estimation and SVA) method requires roughly only 3LmvM2v + lOM3v +MmvKLmv 100% 7LJ14~ + 4M 3 + M~KLc

(12)

of the number of flops performed by Cadzow's method in the preprocessing procedure. Since only the highest order terms are considered the expressions above only hold for large L, M. For complex matrices the same factor can be used (depending on the algorithm used for the complex SVD this ratio can alter a little bit). Although flop counting is a crude approach to the measuring of algorithm efficiency [8] it gives us nevertheless a good idea of the possible gain in computational efficiency. For example, assume that N = 25 and K = 2 (as in [ 16, Example 2] and Example 2 of this paper). Choose Me = Lc = 13, Lmv = 22 and Mmv ~- 4. Then (12) implies that the MV estimation method requires 8% of the computation time needed for Cadzow's method. In

337

practical computations using Matlab [ 15 ], we obtained a ratio of 12.8% for real matrices and a ratio of 14.9% for complex matrices. Of course, since N, M, L are small (12) is too crude. If N = 128 and K = 5, as in Example 1 given further on, then the work involved in preprocessing a 119 × 10 data matrix H by the MV estimation method should, according to (12), only be 2% of the work involved in processing a 65 × 64 matrix as done in Cadzow's method. With Maflab, ratios of 3.6% and 4.4% were obtained for real and complex matrices, respectively. Finally, the condition that the noise should be white (e.g. in applying an SVA or MV estimation method) is not a restriction. If the noise is not white then W RW = / . d is no longer valid, so that the singular value correction function is no longer correct. However, if ,Yw= W " W is known or can be estimated up to an unknown scalar, then we can factor ~ w = R S R (or W=QR, QHQ=I if w is directly available) where R ~ ~M×M is a nonsingular (false) square root. In this case, we consider the matrix ffI= H R - I = h + ~ I = IYlR - I + W R -1.

Observe that W " W is again a multiple of the identity matrix I so that we can again apply the same procedure, i.e. compute the SVD of/-t, truncate to rank K and correct the singular values ( according to (10) or ( 11 ) ) in order to obtain/1K. The only modification is that/dK should be postmultiplied with R, i.e. IlK = fflKR, before restoring the Hankel structure. Another approach that takes account of the correlation structure of the noise has been presented in [ 17]. The fact that WHW or an estimate should be known up to a scalar is not always a trivial requirement. In some cases, this matrix can be estimated during the absence of the signal. For instance, when modeling damped exponential signals as is the case in NMR spectroscopy [4, 20], W can be computed from the tail of the signal where the exact signal has died out. Another remark concerns the numerical computation of the SVD of H R - I . An explicit formation of the inverse R - ~followed by the explicit calculation of the product HR - 1 may result in dramatic loss of accuracy in the data. This can be avoided by using the quotient SVD (called the generalized SVD in [8] ), which Vol. 33, No. 3, September 1993

338

S. Van Huffel / Enhanced resolution and exponentials

delivers the required factorizations without forming quotients and products.

PROPERTY 1. The data matrix Z has a Hankel structu re. PROPERTY 2. The data matrix Z has rank K.

3. Signal enhancement properties 3.1. Composite property mapping Signal enhancement algorithms make use of signal properties that identify the underlying useful information in the data set. Most algorithms require that these properties correspond to closed and convex sets or be describable by linear projection operators. Unfortunately, such restrictions are often too severe for signal properties of interest as exemplified by the SVD characterization of data matrices. A recent signal enhancement algorithm, developed by Cadzow [2], removes these restrictions and only requires that the signal properties be closed. Our signal enhancement algorithm Pestim, outlined in Section 2.1, is essentially a specific variant of Cadzow's method. In essence, the original problem is decomposed into P simpler subproblems related to each individual signal property. The basic requirement of this simplified approach entails the ability to generate the following property solution sets: F p ( Z ) = a r g inf d ( Z - Y ) ,

p = l ..... P,

(13)

YESp

for an arbitrary signal (matrix) Z. d(. ) is a non-negative real-valued metric and Sp specifies the set Sp = {Z ~ _ ~ ] Z has property p }..2" is the universal signal set and is composed of all signals that can arise in a given application, e . g . . ~ = ~ L × M is the set of all complex L × M matrices. Fp must be interpreted as a mapping operator which maps the signal Z into the set of optimum signal approximations. This mapping need not be a point-to-point mapping but is a point-to-set mapping due to the non-uniqueness of solution possibility. In the exponential data modeling problem, the data set {Xo..... xN- 1} is assumed to be exactly representable as a Kth order exponential signal. According to [ 3, Theorem 2] this is the case if and only if the associated Hankel structured data matrix (3) has exactly rank K. Hence, the following properties can be used. SignalProcessing

The matrix operators F1 and F2 in (13), which correspond to Properties 1 and 2, are specified by the third and second step of the preprocessing procedure of the algorithm P-estim. The rank K operator F2 = F ~K) is specified by (5) whereby the M - K smallest singular values in the SVD representation of the data matrix H are set to zero and the largest K singular values are possibly corrected. The special signal enhancement algorithms discussed here differ in the choice of this rank K operator possessing different properties depending on the method used (see Section 3.2). From [2, 3] it follows that the mapping F ~r) is closed. The Hankel property mapping FI = Fh is also closed and readily implemented as outlined in (6). As shown in [2], this mapping Fh(Z) seeks for a Hankel structured matrix that lies closest to Z and is given by the minimizer ;~h of the following optimization problem: min

Z~ ~ St

IIZ- zh IIF,

where S1 = {Zh ~ ~L ×MIZh has a Hankel structure }. The preprocessing procedure of the algorithm Pestim which computes a rank K Hankel approximation of the original data matrix (3) is thus generated by the iterative procedure

H(i) =FhFtK)(Hfi-I)),

with H~o) = H and i>~ 1.

One first gets the rank K approximation of the data matrix H. The corresponding matrix F ~m (H) is generally non-Hankel in structure. To recover the required Hankel structure, we next apply mapping Fh tO matrix F
S. Van Huffel / Enhanced resolution and exponentials

kel structured data matrix of rank K, usually in a rapid fashion. Furthermore, the resulting enhanced data matrix has data elements that generally provide a better representation of the underlying signal components than did the original data. The enhancement process has therefore effectively rejected noise which contaminates the original data. 3.2. The SVD and signal subspace estimation

From ( 1) it follows that we can write the observed data matrix H in (3) as (14)

H=tYl+W,

where B and W represent the exact and pure noise Hankel structured matrices derived from ~ and Xw in ( 1), respectively. A detailed derivation of the algebraic and geometric conditions that allow to derive the exact model from the SVD of a data matrix is given in [7]. Using a multivariate version of the classical Pythagorian lemma for triangles it is shown that the original exact signal can not be recovered consistently when the exact data are perturbed by additive noise. This asymmetry in the consistency of the estimates of the model and the signal is often overlooked. Assume that

339

consistently, while the original exact column space of H, represented by 01, cannot be recovered from the SVD of H, even not asymptotically. The success of SVD-based subspace methods critically depends on these three assumptions, which in practice however are never satisfied exactly (except maybe for condition C3). A nice feature of the SVD however is its robustness with respect to (mild) violations of these conditions. If for instance UB"wIIis small (where I1"II is any unitarily invariant norm), the SVD of H will still deliver good approximations to Range(~71) and Range(17'2). The smaller IIBHWII gets, the better will be the approximations. Despite the fact that/.~ cannot be consistently estimated from the SVD of H, it is possible to find the minimum variance (MV) estimate of B from the SVD of H if all the assumptions C 1, C2, C3 are satisfied and W H W = ~r~l;

(16)

this implies that W itself is an orthogonal matrix and every column of W has norm ~rw. For ease of understanding we reformulate here the MV estimation problem as given in [7]. M V PROBLEM. Given the L X M matrix H as in (14) with rank(B) = K < M . Find the matrix T that mini-

mizes

l-)=0.~l?u= [0, 021 [ ~' 0 ]K ['~/t V2]n K

0 K

0

K

min

TE~'~MxM

IIHT-IZllI2F.

(17)

(15) is the SVD of the exact rank K matrix H. Under the following three conditions: C 1: the exact data should be orthogonal to the noise in the sense that HHW = 0; C2: the matrices }'~ and I72must be orthogonal to each other in the inner product WHW, i.e l)l W HWV2= 0. This is for instance the case when W H W = !*I,/z a scalar; C3: the smallest signal singular value o'r of ,~ must be larger than the largest noise singular value ¢rK+ 1 of ~2 in the SVD (4) of H. Otherwise we cannot separate Range(I) l) from Range( 17'2); the rank K of B and the subspaces generated by ~'1 and t72 (the row and null space of H) can be estimated

Setting to zero the derivatives of the object function trace[THHHHT+/~HB-- 2TnHHB] with respect to the elements of T results in T = (HHH) - 1HHH,

(18)

so that the MV estimate of B is given by H T = H ( H n H ) - 1HH/q.

(19)

Geometrically, the MV estimate of B given H is the orthogonal projection of B onto the column space of H. Observe that rank(HT) = rank(~) = K. Despite the fact that B is not known it is possible to compute the MV estimate (19) from the SVD of H. Indeed, provided C 1, C2, C3 and (16) are satisfied and using (14), (15), the SVD (4) of H is given by VoL 33, No. 3, September 1993

340

S. Van H u f f e l / E n h a n c e d resolution a n d e x p o n e n t i a l s

H=IYI + W=O,Y~, (ZH + WV, (,'~ + WV2 (,"H 2 = [ ( 0, $, + wg,) (~,~ +,rwI,,)

0 1 if,

,,w,,O 0

=[U l

U2][~01

--1/2

OrwlM_KA ;2][VI

~ --I wv2,rw ]

v2]"

(20)

V2] H.

Substituting (20) in (19) yields the desired MV estimate:

HT= UU"t:I= U, U~ (], $, (Znl = Ul($12+ O'wlx)2

min IIH-H~II 2,

-1/2

where $2 = {H r ~ ~L × u IHx has rank K}. The solution is a well-known result [ 8 ], originally due to Eckart and Young, and follows from the SVD (4) of H as

Vn~.

(21)

From (20), (21) it can be seen that we do not obtain a consistent estimate of the column space o f / ~ since U, ~ UI, i.e. the estimate is always biased. The singular values of the MV estimate (21 ) are given by .i,~(~:~

+

2 2 -2 O.wlK)--l/2=lK(IK+Orw$l ) 1/2

(22)

= (,~l ,~ 11),,~1 = C,~1,

where C are the cosines of the canonical angles between the exact column space Range(01) of [ / a n d the column space Range(U1) of the noisy data matrix H generated by its first K left singular vectors. These canonical angles characterize the bias and depend on the signal-to-noise ratio defined by &x/O'w, as shown in (22). The larger this ratio the smaller will be the

Table 1 Exact parameter values of Example 1; a.u. means arbitrary units and $~= ~bll80/'tr expresses the phase in degrees Peak k 1 2 3 4 5

A (Hz) - 1379 -685 - 271 353 478

Sign~Proc¢~ing

dk

(23)

HK ~ $2

= u, $:,(~,~ + ~ i,,) -'/2v 1" = U, (~,~ - O'~IK)~,'(-'

canonical angles. In the noise-free case, for trw = 0, all angles are zero. It is remarkable that the squared diagonal matrix (22) of singular values of the MV estimate (21) is precisely the optimal weighting matrix used in weighted subspace fitting methods to minimize the estimation error variance on the directions-of-arrival [ 21 ]. The MV estimate (21 ), as used in the MV estimation method, is of course not the only possible estimate of the exact signal-only matrix/). Another estimate, as used in Cadzow's method and the LS estimation method, is obtained when we approximate H by a matrix of rank K in least squares sense:

(Hz)

208 256 197 117 808

at (a.u.)

~bk(deg.)

6.1 9.9 6.0 2.8 17.0

15 15 15 15 15

H , , = U , ~ 1V H.

(24)

Observe that the left and right singular vectors of the LS estimate are the same as those of the MV estimate (21), but the singular values are different. In summary, the mapping operator r,~ mv, (x) as used in the MV estimation method, is defined by (17) and F(~)(H) is given by (21) under the conditions C1, C2, C3 and (16). The mapping operator F ~x), as used in Cadzow's and the LS estimation method, is defined by (23) and F ~ ) ( H ) is given by (24). Similarly, it can be shown that the mapping operator, F ~vK],as used in the SVA method, is defined by min IIH- W--HxII2F,

(25)

HKES2

where W is the noise matrix. Provided C1, C2, C3 and (16) are satisfied, F ~vr2(H) follows directly from the Eckart-Young theorem [ 8 ] applied to H - W and is given by nsva = U l ( , ~ l2 - - O__2 r ) \ 1/21rrH wlK F 1•

(26)

The conditions C1, C2, C3 and (16) are never exactly satisfied, but it is shown in [7] that asymptotically, as the number of rows L ~ ~, these conditions are satisfied, i.e. (plim means probability limit)

S. Van Huffel / Enhanced resolution and exponentials

4OO

341

Ori~r~i(-) / N~isy(:) / Recm~m~i~i(-.) i ~ l by PO-HTLS

(a)

~i iii!7 ~i iiiiil 1000

~

O

-500

-1500

-1000

Frequency (I-17.)

ori~.~(-) / No~:) / R,~o.~.,~,,~(-.) . ~

40O

by CADZ.OW-HTLS

(b) 3.~

7 1500

1000

500

0

-500

--1000

"1500

Frequency (Hz)

Original(-) / Noisy(:) / Reconstructed(-.) ilpal by MV-HTLS 400

350

(c)

,!!~i 15o

il

i

i ~iir

1 i

100

iiiiii iJ

,,

Frequency (Hz)

Fig. l. The reconstructed signal (dashed line) of Example l ( N = 128) as estimated by (a) P0-HTLS, (b) C A - H T L S and (c) M V - H T L S . The full line represents the exact signal, the dotted line represents the noisy signal (o',, = 2). Vol. 33, No. 3, September" 1993

S. Van Huffel / Enhanced resolution and exponentials

342

Failure rate for signal SIMU1 (using 5 kinds of P-HTLS) 90 80 70

CA-HTLS: SV-HTLS: MV-HTLS: LS-HTLS: P0-HTLS:

+ * x o .

......,

---. -. :

. ....-

,4- . . . . .

60

/... •" / •

40

.-

30

..-

s

..."

,

/ / 0

J /s/,

,0"" .,,'""

/

50

g~

.- .-'" /

"7

,. .,o"

"*" . . . . . . .

."

j

.,'

,'f'-".,

, ,."

,:'f"

a~

..

20 ..,

+.,,,

,,-

..,,..-.

•

,.4~ 7. . . . . .

10 .,.,..

0

1.2

1.4

.,,

1.6

i

i

i

i

I

i

1.8

2

2.2

2.4

2.6

2.8

3

Noise standard deviation Fig. 2. Percentageof times that P0--HTLS,CA-HTLS, LS-HTLS, SV-HTLS and MV-HTLS fail to resolve all peaks of Example 1 versus the noise standard deviation or,,.N= 128.

C l ' : plim 1 ~ H w = 0; p~OO

p

C2': plim I W n W = o.2i ' which implies that p ---~ ao

p

perfectly well conditioned [ 8], it follows from the perturbation theory for singular values that with this additive noise the noise threshold will become more and more pronounced for increased L. This allows to estimate the noise variance o .2 from the SVD ( 4 ) - ( 2 0 ) of H with o.w =Lo.2 as follows:

plim I l~,HWHW~2 = 0 ; p ---) o¢

p

2=

cry C3':pplimAi(;HnH)=Ai(~/4H/4)

+ O'2,

1

M

L(M-K)

~-+ o.2,

(27)

i=K+I

and the exact singular values ~; of ~ are estimated by with A;(M) the ith eigenvalue of matrix M; 6 - , = f ~ Z - L o -2, provided the elements of the noise vector Xw, which generate the noise matrix W, have zero mean, bounded

fourth moments and are independently and identically distributed with possibly unknown variance 0"2. This means that as L ~ oo we gradually approach the geometric conditions where / t a w = 0 and W H W = htl, implying that the SVD of H g i v e n in (4) also gradually approaches the SVD (20). Since singular values are Signal Processing

i = 1 . . . . . K.

(28)

Together, (27) and (28) allow to obtain an estimate of the signal-to-noise ratio 6"x/o.w and of the canonical angles (22) that characterize the bias of the signal subspace estimate U1. The lower this ratio the larger the canonical angles and hence the larger will be the bias of the signal subspace estimate. Summarizing, it is shown here that it is impossible

S. Van Huffel / Enhanced resolution and exponentials

343

Frequency of Peak No.4 for signal SIMU1 (using 5 kinds of P-HTLS) 10 2

(a)

CA-HTLS: SV-HTLS: MV-HTLS: LS-HTLS: P0-HTLS: CR-bound:

+ * x o .

---. -. :

101 m o~

I0 o

-4

-~

;

,

,

~

~

,

6

8

SNR ( d B ) A m p l i t u d e of Peak No.4 for signal SIMU1 (using 5 kinds of P-HTLS)

101

(b) CA-HTLS: SV-HTLS: MV-HTLS: LS-HTLS: P0-HTLS: CR-bound:

10 0

+ * x o .

---. -. :

~--":.~.i ......

.... a _ - ~ - ~ < ~

,o.~

-4

, -2

~

~

, 4

, 6

8

SNR (riB) Fig. 3. (Relative) R M S E values, obtained for estimating (a) f4 = 353 Hz and (b) a4 = 2.8 of Example 1 by P0-HTLS, C A - H T L S , LS-HTLS, SV-HTLS and M V - H T L S , versus the peak SNR. The solid line is the C R bound.

f r o m the S V D o f H to o b t a i n e d a c o n s i s t e n t e s t i m a t e o f

c o l u m n s b e l o n g to the c o l u m n s p a c e o f / ~ a n d o n l y an

the original s i g n a l . ~ w h i c h are the c o l u m n s of/-t. T h e s e

LS or M V e s t i m a t e can b e obtained.

Vol. 33. No. 3. September1993

344

S. Van Huffel / Enhanced resolution and exponentials

Singular Values of SIMU1 (noise std=0.1, K=5, N=L+M-l=128) 140

i

**f ,,,*'*"'*

120

*,******** "***,**********%,°

sigma4

"'*'*,,,

*** ***

100

** **%

80

g ~

ao sigma5

"7 ~

•

40

*******************e~**********

**********'*'*

•

•

"**'**,,.,..**

•

19 • ***** **

,

***% **

,

20 , ,

0

0

*

.

,..

*

sigma6 20

40

60

**

80

1O0

, ,

120

140

Number of Columns : M Fig. 4. The 4th, 5th and 6th singular values of the ( 129 -M) ×M Hankel matrix built up with 128 data points, modeled by Example 1 and perturbed by white Gaussiannoise with 03= 0.1.

3.3. Signal-to-noise ratio and the size of the Hankel matrix In practice, the number of data points N is limited because the signals die out rapidly (in case of damped exponentials) or because of lack of measurement time. Suppose that we have N data points xn available. Those are put into an L × M Hankel matrix H with N = M + L - 1 where L/> M > K. We can ask the question which choice of (L, M) is best. The geometricalstatistical analysis presented in Section 3.2 suggests choosing L as large as possible. If L increases, the singular values of the pure noise Hankel matrix W decrease weakly until their existence ends by lack of dimensions. The remaining ones come close to each other, making the noise spectrum approximately isotropic. However, as M decreases (i.e. L increases since N = L + M - 1 is constant) the singular values o f ~ also decrease, the smallest one quite strongly ifM = K, since Signal Processing

//contains less elements. This implies that the signalto-noise ratio 6"K/O'K÷~ decreases because the gap between the smallest signal singular value 6"K and the largest noise singular value trK+ ~ of the L X M Hankel matrix H considered in (3) narrows enlarging the bias (characterized by the canonical angles in (22)) of the signal subspace estimate o f / t . The (quick) decrease of the smallest signal singular value of H may have a deteriorating effect on the accuracy of the estimated signal subspace spanned by the columns of U~ and subsequently affect the accuracy of the signal parameter estimates. This implies that M cannot be too small. On the other hand, M may not be large since then the signal-only data cannot be 'orthogonal' to the noise implying that the singular value correction function is no longer valid. For that reason, we often choose M = 2K. Rationale for this choice of M is also provided in the following section.

S. Van Huffel / Enhanced resolution and exponentials

345

Table 2 B i a s + s t a n d a r d d e v i a t i o n v a l u e s o f t h e p a r a m e t e r e s t i m a t e s o f e a c h p e a k o f E x a m p l e 1; N = 128 a n d o I = 1.8 k

Method

1

CA-HTLS

- 0.079 + 5.166

1.616 + 3 3 . 3 2

0.070 + 0.639

0.128 + 5.884

SV-HTLS

0.987 + 4.929

- 0.891 + 31.13

- 0 . 2 7 0 + 0.573

- 1.020 + 5.611

-0.512

2

4

5

dk ( H z )

ak ( a . u . )

~k ( d e g . )

MV-HTLS

1.371 -I-4.961

0 . 1 8 7 -i- 31.35

-t-0.564

- 1.395 + 5.645

LS-HTLS

0.331 + 4 . 9 6 3

- 1.558 + 30.96

- 0.011 + 0 . 5 7 8

- 0.262 + 5.682

P0-HTLS

- 0.425 + 5.095

1.099 __. 3 2 . 0 0

0,025 _ 0.591

0.398 + 5.758

CA-HTLS

0.272 + 4.057

- 0 . 1 0 4 -I- 2 4 . 2 6

0.040 5:0.671

- 0.344 + 4.019 - 0.420 + 3.960

SV-HTLS

0.462 + 4.026

- 0 . 0 4 8 + 24.17

- 0.171 ___0 . 6 5 8

MV-HTLS

0.506 + 4.046

- 0 . 1 0 4 + 24.25

- 0.383 + 0.652

- 0.411 + 3.975

LS-HTLS

0.459 + 4.073

- 0.211 + 2 4 . 0 2

0.050 + 0.662

- 0 . 4 3 8 -t- 3.955 - 0 . 0 6 3 -I- 4 . 0 0 7

P0-HTLS 3

f~ ( H z )

0.202 + 4.118

0.330 + 24.40

0.041 + 0.668

CA-HTLS

- 0.236 + 4.810

0.755 + 2 7 . 9 6

0.113 + 0.583

- 0.116 + 6.524

SV-HTLS

- 0.431 5 : 4 , 7 3 6

1.894 5 : 2 7 . 9 2

- 0.090 + 0.568

- 0 . 1 1 2 -I- 6 . 3 6 6

MV-HTLS

- 0.564 _ 4.723

2.382 + 27.42

- 0 . 2 6 2 + 0.555

- 0.051 + 6.345

LS-HTLS

- 0.417 + 4.756

0.723 + 27.90

0.107 + 0.587

- 0 , 0 0 6 + 6.501

P0-HTLS

- 0.486 + 4.876

0 . 7 2 4 + 27.51

0.091 __.0 . 5 8 4

0 . 3 9 0 + 6.605

CA-HTLS

- 0.411 + 10.01

1.42 + 72.18

0 . 3 5 0 + 1.263

0 . 0 3 2 + 22.91

SV-HTLS

- 2.583 + 13.30

- 13.70 + 6 6 . 3 8

- 0 . 0 8 3 + 1.211

2.351 + 30.57

MV-HTLS

- 2 . 1 4 0 + 12.82

- 16.67 + 6 4 . 8 7

-0.214-1-1.146

1.022 + 2 9 . 8 7

LS-HTLS

- 1.623 + 13.50

- 8.11 5 : 6 9 . 6 7

0 . 1 5 4 + 1.341

P0-HTLS

0 . 1 2 8 + 13.46

11.76 5 : 7 6 . 9 5

0 . 6 1 0 + 1.495

1.310 + 28.03

CA-HTLS

1.82 + 2 0 . 6 8

- 7 . 9 0 ___ 129.3

- 0.069 + 2.336

- 0.171 + 6.581

SV-HTLS

- 10.76 + 20.73

- 2 2 . 1 2 + 115.3

- 0 . 0 5 9 + 2.163

2.225 + 6 . 1 7 2

MV-HTLS

- 13.19 + 2 0 . 1 9

- 18.31 5 : 1 1 2 . 5

- 0 . 2 9 8 + 2.061

2.840 5:6.003

- 5.31 + 21.95

- 3 0 . 1 6 + 121.1

0.063 + 2.304

1.352 + 6.493

6.49 + 2 4 . 5 8

- 28.87 + 128.8

- 0 . 4 3 9 __. 2 . 4 7 0

- 0.394 + 7.332

LS-HTLS P0--HTLS

4. Results The following examples are given to illustrate the differences in resolution enhancement that may accrue from the use of the different signal enhancement algorithms in the problem of fitting complex data to a sum of damped exponentials defined by (2). Simulation procedure

N data points, uniformly sampled atfs Hz (such that tn = n/f~s), are exactly modeled by a Kth order model

function given by (2). These data points are perturbed by white Gaussian noise whose real and imaginary components have standard deviation tr,. (Relative) root mean-squared error (RMSE), bias and standard deviation values of the estimates of all signal parameters fk, dk, ak and q~k are computed using 200 or 1000 noise realizations (excluding failures) at each considered o-~..A failure occurs if not all p e a k s are resolved within specified intervals lying symmetrically around

1.483 + 2 9 . 9 8

the exact frequencies. The signal parameter estimation algorithm considered here is HTLS as outlined in Section 2.1, where H in (7) is as square as possible, i.e. lf,l = N~ 2 and L =/~/( or ~ + 1 if N is even). If no signal enhancement is used the method is called P0-HTLS. If Cadzow's method, LS estimation, SVA or the MV estimation method is used as signal enhancement algorithm prior to HTLS, the method is called respectively CA-HTLS, LS-HTLS, SV-HTLS or MV-HTLS. Varying o-~.and number of data points N are considered, as well as a varying estimated model order Ke~t and size L × M of the Hankel data matrix. Unless stated otherwise, the estimated model order Kest is set to K and the 2 is used in the singular value exact noise variance or,, correction function of the MV estimation and the SVA method. The Cramrr-Rao (CR) lower bounds are derived from the exact parameter values and 01. The SNR at peak k is defined by 10 log(a2/(2tr 2) ). Two different examples are described now.

VoL 33, No. 3, September 1993

S. Van Huffel / Enhanced resolution and exponentials

346

Failure Rate for signal SIMU1 (using 4 kinds of P-HTLS) 35

.... ............. ~, .............•

•./

..

,....

. . . .... .

30 O.0 !

25

;

*!

o...

pt

m

&

,t

,

"~.

" . o ! !J / ~

2O

.. . . . . . -'It"

Ui~

..... I~

a~

•~

CA-HTLS: MV-HTI.S: LS-HTLS: P0-HTLS: noisestd = 2

"

15

x" /'x x

.

. ~ ....

.°~r

-'"

+ -x-. o -. . :

" -I¢""

10 0

10

i

I

i

i

i

20

30

40

50

60

70

Number of Columns for Preprocessing Fig. 5. Percentage of times (based on 500 runs) that P0-HTLS, CA-HTLS, LS-HTLS and MV-HTLS fail to resolve all peaks of Example 1 versus the number of columns M used in the ( 129 - M) × M data matrix. In CA-HTLS a 65 X 64 matrix was considered throughout. N= 128 and ~,. = 2.

EXAMPLE 1. A fifth order model function ( K = 5), given in Table 1 and shown in Fig. 1, is considered. fs = 10 kHz; 200 noise realizations are taken at each considered o-~.and the size of the Hankel matrix considered in the preprocessing procedure of L S - H T L S , SV-HTLS

and

MV-HTLS,

is

(N-K, st-4) X

( K e s t + 5 ) , e.g. for N = 1 2 8 and K~st=K a l 1 9 X 1 0 matrix H is considered ( C A - H T L S uses a 65 X 64 matrix). The halfwidth o f the frequency intervals in which all peaks should be resolved are respectively 82, 82, 82, 43 and 82 Hz where 43 and 82 are the CR lower bounds for f4 and f5 at the noise standard deviation o'v = 8.8 at which the intervals of peak 4 and 5 are coinciding. This example describes a typical in vivo

31p N M R

signal [5], where a lot of spectral peaks are to be resolved using a limited number of noisy observations. Signal Processing

Usually, the signal frequencies are known a priori, but the damping factors, amplitudes and phases need to be computed as accurately as possible in order to determine the concentrations of the chemical components of interest. The following conclusions can be drawn. The use of one of the signal enhancement algorithms described above improves the frequency resolution as illustrated in Fig. 2. It is very clear that the enhanced resolution is most pronounced for the M V estimation method M V - H T L S . The S V A method exhibits a comparable slightly worse resolution performance. Cadz o w ' s method C A - H T L S has a much higher failure rate but still performs better than the non-enhanced method P 0 - H T L S . Although L S - H T L S exhibits only a slightly better resolution than C A - H T L S , it should be preferred to C A - H T L S because o f its much better computational efficiency. Like M V - H T L S and S V -

347

S. Van Huffel/ Enhanced resolution and exponentials

Failure Rate for signal SIMUI (using 5 kinds of MV-HTLS) 45 /

40 35

Estimated Estimated Estimated Estimated Estimated

Order=K+0 Order=K+l Order=K+2 Order=K+3 Order=K+4

+

--

X

-.

0

-.

to / ¢' ff

ffff

"O

t~

/

30

t tI

ff

-. "-.

te ,s

25

m /

20

/

tt

•

•

•

it ie

e,,t,.

/

,

t

t

t

.o

.."

.'

IL

..

.

"..

...

....

15 a~

10

¢ ff

5

•

- ~ ' ~ - ~ - - ~ -...~-~.~ - ~'- ~ ''~':..T "°" I,. ......

1.2

1.4

1.6

"'"

1.8

/

.-":-6:

.....

-'.m ." . . - .y.

i

2

2.2

I

i

I

2.4

2.6

2.8

Noise Standard Deviation Fig. 6. Percentage of times that MV-HTLS fails to resolve all peaks of Example 1 versus the noise standard deviation tr,. for varying model order estimation Kest= 5, 6, 7, 8 and 9. N= 128. Here the noise variance used in the MV-HTLS singular value correctionfunction is computed from the average of the squared noise singular values. HTLS, L S - H T L S requires only 4% of the computation time needed for C A - H T L S . Furthermore, all detected failures in resolution are due to the fact that one of the overlapping peaks, especially peak 4, cannot be resolved. It is important to note that the occurrence of a failure of M V - H T L S is mostly due to the fact that the estimated peak 4 (or 5) falls outside the specified frequency interval but this peak lies still in the neighbourhood. This is in contrast to Cadzow's method C A - H T L S which has more failures. If a failure occurs then we mostly have the situation that peak 4 and 5 are estimated by only one peak, as illustrated in Fig. 1, and an additional noise peak appears far away from the exact frequencies, often around + 4000 Hz.' If we look at the accuracy o f the estimatedparameters we observe that the differences in accuracy between the different methods are small, especially for the estimates of peak 1, 2 and 3. As illustrated in Fig.

3, the damping factors and amplitudes of peak 4 and 5 and also the phase of peak 5 are better estimated by M V - H T L S , but the frequencies are estimated more accurately by C A - H T L S (mainly due to a smaller bias). Note that these results are based on average values over all good runs. This number differs from one method to the other. It certainly reflects the fact that M V - H T L S is still able to compute accurate estimates of all parameters during the 'difficult' runs in which C A - H T L S , L S - H T L S and P 0 - H T L S fail to resolve peak 4 or 5. Even if the results are averaged over the same runs, the relative differences in accuracy remain the same. This loss in accuracy is due to the occurrence of a smaller singular value gap irK/trK+ ~in the 119 × 10 data matrix H used by M V - H T L S compared to the One in the 65 × 64 matrix used in C A - H T L S , as illustrated in Fig. 4. As explained in Sections 3.2-3.3, a smaller gap means that the canonical angles between the exact signal subspace Ul and the computed one U1 enlarge, Vol. 33, No. 3, S e p t e m b e r 1993

S. Van Huffel / Enhanced resolution and exponentials

348

Failure Rate for signal SIMU1 (using 5 kinds of P-HTLS)

80

(a) MV-3CADZOW-HTLS MV-2CADZOW-HTLS MV-CADZOW-HTLS MV-HTLS CADZOW-HTLS

70

60

+ * x o .

-. . . . . . . . . . . . . . . . . . -. -. ... : "

v

50

/ :.

m

D..........

.:.

4O

..""

f

30

..

.""

e-"

jr...

,):~f"

20

10

It ..... ~':2"'~ ~ 2 ~ ' ; i ~ ' : : ~ ' : l t ~ ' ' - -

1.2

1.4

1.6

'

'

1.8

2

2.2

2.4

2.8

2.6

3

Noise Standard Deviation Frequency of Peak No.4 for signal SIMU1 (using 5 kinds of P-HTLS) 10 2

(b)

MV-3CADZOW-HTLS MV-2CADZOW-HTLS MV-CADZOW-HTLS MV-HTLS CADZOW-HTI.~ CR-bound

+ * x o .

---. -. :

o 101

10,

-4

,

,

,

,

,

-3

-2

-1

1

4

6

SNR (dB) Fig. 7. Comparison.of C A - H T L S , M V - H T L S and the iterated methods MV--CADZOW-HTLS, M V - 2 C A D Z O W - H T L S and M V - 3 C A D Z O W HTLS (i.e. M V followed by 1, 2 and 3 Cadzow iterations ) for Example 1: (a) percentage of times that each method fails to resolve all peaks versus the noise standard deviation o'v; (b) R M S E values, obtained for estimating f4 = 353 Hz versus the peak SNR. The solid line is the CR bound.

Signal Processing

S. Van Huffel / Enhanced resolution and exponentials

349

Damping factor of Peak No.4 for signal SIMU1 (using 5 kinds of P-HSVD) 101 CA-HSVD SV-HSVD MV-HSVD LS-HSVD P0-HSVD

+ * x o

---. -.

t~ taO oz~.

...

tOO m

"~*'-~Z~__._.,._-.~. ".... [..............

10-1 -L

I

i

i

i

i

i

I

i

/

-3

-2

-1

0

1

2

3

4

5

6

SNR (dB) Fig. 8. Relative RMSE values, obtained for estimating the damping factor d4 = 117 Hz of Example I by CA-HSVD, SV-HSVD, MV-HSVD, LS-HSVD and P0-HSVD, are plotted against the peak SNR and comparedwith the CramCr-Raobounds (solid line). thereby also increasing the bias of the estimates as shown in Table 2. Especially the accuracy of the most difficult peak (smallest in amplitude, overlapping .... ) is affected. For Example 1 this is peak 4. Nevertheless it was also observed that increasing the number of colunms M of our ( N - M + 1 ) × M Hankel matrix from 6 to 64 (square), first decreases and then increases- from M = 15 on - the failure rate of M V - H T L S from 10% to 20%, while the accuracy of all estimates remains more or less the same for all M values (see Fig. 5). P 0 HTLS and C A - H T L S exhibit a more or less constant failure rate and parameter accuracy. Therefore we can conclude that we better take M small ( = 2K) since a higher M does not improve the performance. Observe that M V - H T L S and S V - H T L S require the estimation of the noise variance in order to correct the singular values of the Hankel data matrix. However, no differences in accuracy or resolution have been observed between the use of the exact noise variance and the use of one of the estimators (based on the tail

of the signal or based on the average of the squared noise singular values). Overestimating the model order (taking the order Kest > 5) further improves clearly the failure rate of all methods, as illustrated in Fig. 6, but the differences in failure rate remain relatively the same. However, except for the frequency estimates of peaks 4 and 5 and the phase estimate of peak 4, the accuracy of the computed parameters does not improve and is mostly slightly worse. The differences in accuracy are smallest for C A - H T L S . Finally, we can apply the signal enhancement algorithm iteratively or combine several algorithms together. As shown in Fig. 7, applying the MV estimation method followed by a Cadzow iteration prior to HTLS, called M V - C A D Z O W - H T L S , reduces further the failure rate and also improves the parameter accuracy of M V - C A D Z O W - H T L S wherever C A HTLS performs better. In particular, the M V - C A D Z O W - H T L S frequency and phase estimates of peak 4 Vol. 33, NO. 3, S e p t e m b e r 1993

S. Van Huffel / Enhanced resolution and exponentials

350

Oril~ul(-),Nolsy(: noi~etd=0.1),Rtconstructed(-) signal SIMU5by P0-HTLS 1211

(a)

100

411

211

0.1

0.12

0.14

0.16

0.18

0.2

0.22

0.24

0.26

0.28

0.3

Frequency(Hz) Original(-),Noisy(:noisestd=0.1),Reconstructed(-)signal SIMU5 by CA-HTLS

(b) 8O

411

0.1

0.12

0.14

0.16

0.18

0.2

0.22

0 . 2 4 0.26

0.28

0.3

Frequency(Hz)

Original(-),Noisy(:n¢ise~td=O.l),Recomtructed(--)signal SIMU5 by MV-HTI.B 120

(c) 100

8O

J

401

J 0.1

0.12

0,14

0.16

0.18

0.2

L 0.22

0.24

0.26

0.28

0.3

Frequency(Hz)

Fig. 9. The reconstructed signal (dashed line) of Example 2 ( N = 2 5 ) as estimated by (a) P0-HTLS, (b) C A - H T L S and (c) M V - H T L S (a 21 × 5 matrix was used). The full line represents the exact signal, the dotted line represents the noisy signal ( o',, = 0.1 ). SignalProcessing

S. Van Huffel / Enhanced resolution and exponentials

351

Failure Rate for signal SIMU5 (using 5 kinds of P-HTLS) 60 [. "" '., ""-. 50 ,, ",.

40

v

"e

30~

CADZOW-HTLS SVA-HTLS M V-HTLS LS-HTLS P0-HTLS

-I-

--

X

=.

O

-o

~

"

',

""

",

O OO

20

", • "',) , , "',

0. •.

10

0

14

•

•

,..

".~:.. -..~ .....

I

I

I

16

18

20

~

~•=

~

22

~

26

SNR (dB) Fig. 10. Percentageof tirnes that P0-HTLS, CA-HTLS, LS-HTLS, SV-HTLS and MV-HTLS fail to resolve all peaks of Example 2 versus the SNR. N = 25 and K e s t = 2. and the M V - H T L S frequency estimates of peak 5 are now as accurate as the C A - H T L S estimates• In contrast, M V - M V - H T L S , M V - L S - H T L S and M V - S V HTLS do not improve the failure rate. This can be explained as follows• After one MV iteration the rank of the obtained Hankel matrix is more or less reduced to rank K and only slightly improves the following iterations, but the Hankel mapping operator Fh, which restores the Hankel structure by averaging along the antidiagonals, still improves the results. Since a Cadzow iteration averages over more elements (for the middle data points a Cadzow iteration averages over 64 elements compared to an MV iteration which averages only over 10 elements), it is able to produce lower standard deviations of the corrected signal data points thereby improving further the accuracy of the parameter estimates. No further improvement in failure rate nor in parameter accuracy was observed when performing more than 2 iterations• Similar experiments have been performed using the

signal parameter estimation algorithm HSVD but, except for the fact that the accuracy of the HTLS estimates is generally better than that of the HSVD estimates (see Fig. 8), the same conclusions hold.

E X A M P L E 2. A second order model function (K = 2), taken from the literature [9, 1 ] and given by ( N = 25) x, = e ( -o.01 + j 2 ~ O . 2 ) t n = peakl + peak2,

"t-

e ( -0.02 + j 2 " t r

0.22)tn

n = 0 , . . . N - 1,

describes two closely spaced peaks as shown in Fig. 9. fs = 1 Hz and 1000 noise realizations are taken at each considered cry. Unless stated otherwise, the Hankel matrix H in the preprocessing procedure of L S - H T L S , S V - H T L S and M V - H T L S has dimensions 22 X 4 (for C A - H T L S a 13 × 13 matrix is used)• The halfwidth of the frequency intervals in which the two peaks should be resolved are 0.0094 and 0.0106 Hz, respectively, and are the CR lower bounds at the noise standard deviation 0-,,=0.2, where the intervals of peakl and V o l . 33, N o . 3, S e p t e m b e r

1993

S. Van Huffel / Enhanced resolution and exponentials

352

Frequency of Peak No.1 for signal SIMU5 (using 5 kinds of P-HTLS) 10-2

CADZOW-HTLS SVA-HTLS MV-HTLS LS-HTLS P0-HTI.S CR-bound

+ * x o . -

---. -. :

~:-~' .............. ~:-':-"~ :i;.' ;. ..... o t~

m

a~

10-3 14

, 16

, 18

, 20

SNR

, 22

, 24

26

(dB)

Fig. 1l. RMSE values, obtained for estimatingf~ =0.2 Hz of Example 1 by P0-HTLS, CA-HTLS, LS-HTLS, SV-HTLS and MV-HTLS, versus the peak SNR. At low SNR the curves fall below the CR bound (solid line) because the estimates of bad runs have been eliminated. peak2 are coinciding. The SNR is defined by 10 log( 1 / (2or 2) ). Similar observations as for Example 1 have been made except that here M V - H T L S exhibits a slightly better or comparable accuracy for all parameter estimates. As shown in Fig. 10, M V - H T L S improves the resolution at all considered SNRs. C a d z o w ' s method C A HTLS has a lower failure rate than L S - H T L S but this difference vanishes when also the runs with negative damping factors ( d l < 0 and d2 < 0) are excluded (in that case the curves are relatively the same but the failure rate is higher). Moreover, L S - H T L S , M V HTLS and S V - H T L S only require 15% of the computation time needed for C A - H T L S , as experimentally observed. If we look at the accuracy of the estimated parameters, we observe that M V - H T L S performs best in all cases followed by S V - H T L S , L S - H T L S , C A HTLS and P 0 - H T L S as shown in Fig. 11. In this example the singular value gap trr/o'r+~ in the 22 x 4 data Signal Processing

matrix used by M V - H T L S is not much smaller than the one that occurred in the 13 × 13 matrix used by C A - H T L S , as shown in Fig. 12. Therefore, the increase in bias, as shown in Table 3, is not as much pronounced as for Example 1. Since M V - H T L S clearly minimizes the variance it results in an overall better performance in terms of variance and MSE. Enlarging the number of data points from N = 18 to 64 or more also enlarges the difference in gap Orr/o'r+ l (although not more than one order of magnitude) between M V - H T L S and C A - H T L S and improves considerably the accuracy of the estimates, but it does not change the relative differences in parameter accuracy. For example, for N = 6 4 and S N R = 2 0 dB, erE/or3 is about 8 for C A - H T L S compared to a gap of 2 in the 61 X 4 Hankel matrix used in M V - H T L S . Also, the resolution improves considerably with the number of data points. For instance, for SNR = 14 dB the failure rate o f M V - H T L S decreases from 80% at N = 18 to less than 5% at N = 30 and to zero at N = 35.

S. Van Huffel / Enhanced resolution and exponentials

353

Singular Values of SIMU5 (SNR=30 dB, K=2, N=L+M-l=25) 20 18 sigmal

,.....,.....*...-*.-.*.....~.....,

16

,..

..,. -....

14

/

>

/.

12 ~0 0~ +

t~

10 8 6 4 sngma2

2

, . . . . . , . . . . * ..... . . . . . . . . ,t.---.¢...-. * . - . . . ~ - . . . . * - - . ' * - . . - ' 1 . . . - * . . " * - * - - * " "

sigma3

0

5

0

10

.

. . ,.... ,.... ,... ¢¢'"" *-'" '~'"-*-"-*--'~'.-'-+'.'-'~' ".': ~

15

20

25

Number of Rows : L Fig. 12. First 3 singular values of the ( 2 6 - M) × M Hankel matrix built up with 25 data points, modeled by Example 2 and perturbed by white Gaussian noise with or,.= 0.0224 (SNR = 30 dB). Table 3 Bias 5: standard deviation values of the parameter estimates of each peak of Example 2; N= 25 and the SNR = 19 dB k

Method

1

CA-HTLS SV-HTLS MV-HTLS LS-HTLS P0-HTLS CA-HTLS SV-HTLS MV-HTLS LS-HTLS P0-HTLS

2

fk (Hz) - 0.00026 -l-0.00397 - 0.00009 + 0.00384 0.00003 5:0.00375 - 0.00017 + 0.00396 - 0.00029 + 0.00410 0.00016 5:0.00444 - 0.00008 5:0.00438 - 0.00019 + 0.00431 0.00004 + 0.00443 0.00016 5:0.00461

dk (Hz)

ak (a.u.)

- 0.00047 + 0.0307 0.00314 + 0.0276 0.00470 5:0.0259 0.00098 5:0.0292 - 0.00067 + 0.0321 - 0.00109 -t-0.0357 0.00202 + 0.0307 0.00325 5:0.0290 0.00006 + 0.0332 - 0.00069 5:0.0355

0.091 + 0.532 0.101 + 0.500 0.107 5:0.497 0.098 + 0.51 ! 0.083 + 0.525 0.0917 + 0.535 0.0771 5:0.500 0.0732 5:0.492 0.0832 -I-0.506 0.0843 -I-0.517

~bk(deg.)

-

2.286 + 29.62 2.599 5:26.27 2.348 + 25.84 2.361 5:27.84 2.439 + 30.44 1.256 5:32.56 1.237 -I-29.88 0.669 + 28.77 0.547 5:31.87 1.895 + 33.01

A slight i m p r o v e m e n t in a c c u r a c y a n d r e s o l u t i o n w a s

t h e n i n c r e a s e s c o n s i d e r a b l y f r o m Kest = 5 on, w h i l e the

o b s e r v e d w h e n e s t i m a t i n g the n o i s e v a r i a n c e f r o m the

a c c u r a c y also slightly d e c r e a s e s w i t h i n c r e a s i n g m o d e l

tail o f the signal ( m e r e l y f r o m the data p o i n t s b e t w e e n

o r d e r Kest.

t157 a n d t256) c o m p a r e d to the u s e o f the o t h e r n o i s e

Finally, i n c r e a s i n g the n u m b e r o f c o l u m n s in M in the ( 2 6 - M )

v a r i a n c e e s t i m a t o r s c o n s i d e r e d in this paper.

× M H a n k e l data m a t r i x f r o m 3 to 13

slightly

( s q u a r e ) i n c r e a s e s the failure rate o f M V - H T L S f r o m

d e c r e a s e s the failure rate o f M V - H T L S ( a f e w % ) an

2 8 % to 4 2 % , as s h o w n in Fig. 13, w h i l e the a c c u r a c y

Overestimating

the

model

order

first

Vol. 33, No. 3, S e p t e m b e r 1993

S. Van Huffel / Enhanced resolution and exponentials

354

Failure Rate for signal SIMU5 (using 4 kinds of P-HTLS) 60 •- . .

. . . ' ' ' . .

.....

...............

55

~,

50

'~

45 •

e~

CADZOW-HTLS MV-HTI,S LS-HTLS P0-HTLS SNR = 14dB, N = 25

40

a,

"'4 .j ,

+ x o .

--. -. :

~

:o ...........

",

~r

o"

,-"

, :/

~...........,~ x

,"

35

313

25

3

"', .

.

.

.

.

.

.

x- " °

i

i

i

i

i

i

i

i

i

4

5

6

7

8

9

10

11

12

13

Number of Columns for Preprocessing Fig. 13. Percentage of times that P0-HTLS, CA-HTLS, LS-HTLS and MV-HTLS fail to resolve all peaks of Example 2 versus the number of columns M used in the (26 - M) × M data matrix. In CA-HTLS a 13 × 13 matrix was considered throughout. N= 25 and SNR = 14 dB. of the estimates remains more or less constant or decreases slightly. The same behaviour was observed at other SNRs. Therefore, we can conclude that M = 2K is the best choice for this example in terms of resolution, accuracy and computational efficiency.

5. Conclusions A new signal enhancement algorithm based on minimum variance estimation is presented and shown to improve significantly the resolution of all spectral peaks present in damped exponential signals embedded in white noise. Moreover, arranging the data in a very rectangular. Hankel data matrix reduces drastically the required computation time. Depending on the number of data points and the model order, the resolution can be doubled independently of the noise variance, while the required computation time is less than 10% of the time needed to execute C a d z o w ' s signal enhancement Signal

Processing

algorithm. If the gap between the signal singular values and the noise singular values is not too small, then the resulting signal parameter estimates have comparable accuracy and are generally more accurate than those computed with the non-enhanced signal. The accuracy and failure rate can be further improved by performing an additional Cadzow iteration. Overestimating the model order also improves the failure rate but not the accuracy. The newly presented signal enhancement algorithm can also be applied (after some straightforward changes) to other problems in signal processing, such as sinusoidal modeling and transfer function modeling where the data are ordered in slightly differently structured matrices. It should be noted that the proposed signal enhancement algorithm is not restricted to signals embedded in white noise. Provided that the noise covariance matrix is known or can be estimated up to a scalar multiple (e.g. from the tail o f the signal when modeling damped

S. Van Huffel / Enhanced resolution and exponentials

exponentials), it is shown how to modify the data in order to obtain similar results.

Acknowledgments The author greatly appreciated the assistance of Hua Chen in computing the results presented in the figures and tables of this paper.

References [1] T.J. Abatzoglou and J.T. Reagan, "Statistically efficient estimation of complex frequencies of sinnsoids", IEEE Trans. Signal Process., 1992, In review. [2] J. Cadzow, "Signal enhancement: A composite property mapping algorithm", IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-36, No. 1, January 1988, pp. 49-62. [3] J.A. Cadzow and D.M. Wilkes, "Enhanced sinusoidal and exponential data modeling", in: R.J. Vaccaro, ed., SVD and Signal Processing, II: Algorithms, Analysis and Applications, Elsevier, Amsterdam, 1991, pp. 335-352. [4] H. Barkhuijsen, R. De Beer and D. Van Ormondt, "Improved algorithm for noniterative time-domain model fitting to exponentially damped magnetic resonance signals", J. Magnetic Resonance, Vol. 73, No. 3, July 1987, pp. 553-557. [5] R. De Beer and D. Van Ormondt, "Analysis of NMR data using time-domain fitting procedures", in: M. Rodin, Guest Ed., ln-vivo Magnetic Resonance Spectroscopy 1: Probeheads, Radiofrequency Pulses, Spectrum Analysis, NMR Basic Principles and Progress Series, Vol. 26, Springer, Berlin, 1992, pp. 201-248. [ 6] P. De Groen and B. De Moor, "The fit of a sum of exponentials to noisy data", J. Comput.Appl. Math., Vol. 20, 1987, pp. 175187. [7] B. De Moor, "The singular value decomposition and long and short spaces of noisy matrices", ESAT-SISTA Report 199038, ESAT Lab., K.U. Leuven, Leuven, Belgium, 1990; also IEEE Trans. Signal Process., to appear. [8] G.H. Golub and C.F. Van Loan, Matrix Computations, 2nd Edition, Johns Hopkins Univ. Press, Baltimore, MD, 1989.

355

[ 9 ] Y. Hua and T.P. Sarkar, "Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise", IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-38, No. 5, May 1990, pp. 814-824. [ 10] I. Kirsteins, Analysis of reduced rank interference cancellation, Ph.D. Thesis, Univ. of Rhode Island, Kingston, RI, 1990. [ 11 ] R. Kumaresan and D.W. Tufts, "Estimating the parameters of exponentially damped sinusoids and pole-zero modeling in noise", IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-30, No. 6, December 1982, pp. 833-840. [ 12] S.Y. Kung, K.S. Arun and D.V. Bhaskar Rao, "State-space and singular value decomposition-based approximation methods for the harmonic retrieval problem", J. Opt. Soc. Amer., Vol. 73, No. 12, December 1983, pp. 1799-1811. [ 13 ] F. Li, A unified performance analysis of subspace-based DOA estimation algorithms, Ph.D. Thesis, Univ. of Rhode Island, Kingston, RI, 1990. [14] F. Li, R.J. Vaccaro and D.W. Tufts, "Unified performance analysis of subspace-based estimation algorithms", Proc. IEEE lnternat. Conf. Acoust. Speech Signal Process. (ICASSP), Albuquerque, NM, April 1990, pp. 2575-2578. [15] The Mathworks, Inc., PRO-MATLABfor Unix Computers, Version 3.5e, South Natick, MA 01760, 1991. [ 16] D. Tufts, R. Kumaresan and I. Kirsteins, "Data adaptive signal estimation by singular value decomposition of a data matrix", Proc. IEEE, Vol. 70, No. 6, June 1982, pp. 684---685. [ 17] D.W. Tufts and A.A. Shah, "Estimation of a signal waveform from noisy data using low-rank approximation to a data matrix", IEEE Trans. Signal Process., to appear. [ 18] S. Van Huffel and J. Vandewalle, "Comparison of total least squares and instrumental variable methods for parameter estimation of transfer function models", lnternat. J. Control, Vol. 50, No. 4, 1989, pp. 1039-1056. [19] S. Van Huffel and J. Vandewalle, The Total Least Squares Problem: Computational Aspects and Analysis, Frontiers in Applied Mathematics series, Vol. 9, SIAM, Philadelphia, PA, 1991. [20] S. van Huffel, L. Aerts, J. Bervoets, J. Vandewalle, C. Decanniere and P. Van Hecke, "Improved quantitative timedomain analysis of NMR data by total least squares", in: J. Vandewalle, R. Boite, M. Moonen and A. Oosterlinck, eds., Signal Processing VI: Theories and Applications. Proc. 6th European Signal Processing Conf. (EUSIPCO), Brussels, Belgium, 24-27 August 1992, Elsevier (North-Holland), Amsterdam, 1992, Vol. III, pp. 1721-1724. [21 ] M. Viberg and B. Ottersten, "Sensor array processing based on subspace fitting, IEEE Trans. Signal Process., Vol. SP-39, No. 5, May 1991, pp. 1110--1121.

Vol. 33, No. 3, September 1993