Parameter estimation for exponential sums by approximate Prony method

ARTICLE IN PRESS Signal Processing 90 (2010) 1631–1642

Contents lists available at ScienceDirect

Signal Processing journal homepage: www.elsevier.com/locate/sigpro

Parameter estimation for exponential sums by approximate Prony method Daniel Potts a,, Manfred Tasche b a b

Chemnitz University of Technology, Department of Mathematics, D-09107 Chemnitz, Germany University of Rostock, Institute for Mathematics, D-18051 Rostock, Germany

a r t i c l e in fo

abstract

Article history: Received 14 April 2009 Received in revised form 8 November 2009 Accepted 9 November 2009 Available online 17 November 2009

The recovery of signal parameters from noisy sampled data is a fundamental problem in digital signal processing. In this paper, we consider the following spectral analysis problem: Let f be a real-valued sum of complex exponentials. Determine all parameters of f, i.e., all different frequencies, all coefficients, and the number of exponentials from finitely many equispaced sampled data of f. This is a nonlinear inverse problem. In this paper, we present new results on an approximate Prony method (APM) which is based on [1]. In contrast to [1], we apply matrix perturbation theory such that we can describe the properties and the numerical behavior of the APM in detail. The number of sampled data acts as regularization parameter. The first part of APM estimates the frequencies and the second part solves an overdetermined linear Vandermonde-type system in a stable way. We compare the first part of APM also with the known ESPRIT method. The second part is related to the nonequispaced fast Fourier transform (NFFT). Numerical experiments show the performance of our method. & 2009 Elsevier B.V. All rights reserved.

Keywords: Spectral analysis problem Parameter estimation Exponential sum Nonequispaced fast Fourier transform Digital signal processing Approximate Prony method Matrix perturbation theory Perturbed Hankel matrix Vandermonde-type matrix ESPRIT method

1. Introduction We consider a real-valued sum of complex exponentials of the form f ðxÞ :¼

M X

rj eioj x ðx 2 RÞ;

q :¼ minfoj þ 1  oj ; j ¼ 0; . . . ; Mg:

with complex coefficients rj a0 and real frequencies oj , where

o0 ¼ 0 o o1 o    o oM o p: For simplicity, we discuss only the case r0 a0. Since f is a real function, we can assume that

rj ¼ r j ðj ¼ 0; . . . ; MÞ:

ð1:3Þ

ð1:1Þ

j ¼ M

oj ¼  oj ;

The frequencies oj (j ¼  M; . . . ; M) can be extended with the period 2p such that for example oM1 :¼ 2p þ oM and oM þ 1 :¼ 2p  oM . We introduce the separation distance q of the frequency set foj : j ¼ 0; . . . ; Mg by

ð1:2Þ

Then we have qM o p. Let N 2 N be an integer with N Z 2M þ 1. Assume that the sampled data hk :¼ f ðkÞ ðk ¼ 0; . . . ; 2NÞ are given. Since oM o p, we infer that the Nyquist condition is fulfilled (see [2, p. 183]). From the 2N þ 1 sampled data hk (k ¼ 0; . . . ; 2N) we have to recover the positive integer M, the complex coefficients rj , and the frequencies oj 2 ½0; pÞ (j ¼ 0; . . . ; M) with (1.2) such that M X

rj eioj k ¼ hk ðk ¼ 0; . . . ; 2NÞ:

j ¼ M

 Corresponding author.

E-mail addresses: [email protected] (D. Potts), [email protected] (M. Tasche). 0165-1684/$ - see front matter & 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2009.11.012

By N Z2M þ 1, the overmodeling factor ð2N þ 1Þ=ð2M þ1Þ is larger than 2. Here 2N þ1 is the number of sampled data and 2M þ1 is the number of exponential terms.

ARTICLE IN PRESS 1632

D. Potts, M. Tasche / Signal Processing 90 (2010) 1631–1642

Overmodeling means that the corresponding overmodeling factor is larger than 1. The above spectral analysis problem is a nonlinear inverse problem which can be simplified by original ideas of G.R. de Prony. But the classical Prony method is notorious for its sensitivity to noise such that numerous modifications were attempted to improve its numerical behavior. For the more general case with oj 2 C see e.g. [3] and the references therein. Our results are based on the paper [1] of Beylkin and Monzo´n. The nonlinear problem of finding the frequencies and coefficients can be split into two problems. To obtain the frequencies, we solve a singular value problem of the rectangular Hankel matrix and find the frequencies via roots of a H ¼ ðf ðkþ lÞÞ2NL;L k;l ¼ 0 convenient polynomial of degree L, where L denotes an a priori known upper bound of 2M þ 1. To obtain the coefficients, we use the frequencies to solve an overdetermined linear Vandermonde-type system. Note that the solution of the overdetermined linear Vandermondetype system is closely related to the inverse nonequispaced fast Fourier transform studied recently by the authors [4]. In contrast to [1], we present an approximate Prony method (APM) by means of matrix perturbation theory such that we can describe the properties and the numerical behavior of the APM in detail. The first part of APM estimates the frequencies and the second part solves an overdetermined linear Vandermonde-type system in a stable way. We compare the first part of APM also with the known ESPRIT method. In applications, perturbed values h~ k 2 R of the exact sampled data hk ¼ f ðkÞ are only known with the property h~ k ¼ hk þek ;

jek j r e1

ðk ¼ 0; . . . ; 2NÞ;

where the error terms ek are bounded by certain accuracy e1 40. Also even if the sampled values hk are accurately determined, we still have a small roundoff error due to the use of floating point arithmetic. Furthermore we assume that jrj jZ e1 ðj ¼ 0; . . . ; MÞ. This paper is organized as follows. In Section 2, we sketch the classical Prony method. Then in Section 3, we generalize a first APM based on ideas of Beylkin and Monzo´n [1]. Sections 4 and 5 form the core of this paper with new results on APM. Using matrix perturbation theory, we discuss the properties of small singular values and related right singular vectors of a real rectangular Hankel matrix formed by given noisy data. By means of the separation distance of the frequency set, we can describe the numerical behavior of the APM also for clustered frequencies, if we use overmodeling. Further we discuss the sensitivity of the APM to perturbation. Finally, various numerical examples are presented in Section 6.

If

The classical Prony method works with exact sampled data. Following an idea of G.R. de Prony from 1795 (see e.g. [5, pp. 303–310]), we regard the sampled data hk ¼ f ðkÞ ðk ¼ 0; . . . ; 2NÞ as solution of a homogeneous linear difference equation with constant coefficients.

rj ðeioj Þk ;

j ¼ M

with (1.2) is a solution of certain homogeneous linear difference equation with constant coefficients, then eioj (j ¼  M; . . . ; M) must be zeros of the corresponding characteristic polynomial. Thus P0 ðzÞ :¼

M Y

ðz  eioj Þ ¼ ðz  1Þ

j ¼ M

¼

2M þ1 X

M Y

ðz2  2z cosoj þ1Þ

j¼1

pk zk

ðz 2 CÞ;

ð2:1Þ

k¼0

with p2M þ 1 ¼  p0 ¼ 1 is the monic characteristic polynomial of minimal degree. With the real coefficients pk ðk ¼ 0; . . . ; 2M þ 1Þ, we compose the homogeneous linear difference equation 2M þ1 X

xl þ m pl ¼ 0

ðm ¼ 0; 1; . . .Þ;

ð2:2Þ

l¼0

which obviously has P0 as characteristic polynomial. Consequently, (2.2) has the real general solution P ioj m ðm ¼ 0; 1; . . .Þ with arbitrary coeffixm ¼ M j ¼ M rj e cients r0 2 R and rj 2 C (j ¼ 1; . . . ; M), where (1.2) is fulfilled. Then we determine rj ðj ¼ 0; . . . ; MÞ in such a way that xk  hk ðk ¼ 0; . . . ; 2NÞ. To this end, we compute the least squares solution of the overdetermined linear Vandermonde-type system M X

rj eioj k ¼ hk ðk ¼ 0; . . . ; 2NÞ:

j ¼ M

Let L 2 N be a convenient upper bound of 2M þ 1, i.e., 2M þ 1r Lr N. In applications, such an upper bound L is mostly known a priori. If this is not the case, then one can choose L ¼ N. The idea of G.R. de Prony is based on the separation of the unknown frequencies oj from the unknown coefficients rj by means of a homogeneous linear difference equation (2.2). With the 2N þ 1 sampled data hk 2 R we form the rectangular Hankel matrix H :¼ ðhl þ m Þ2NL;L 2 Rð2NL þ 1ÞðL þ 1Þ : l;m ¼ 0

ð2:3Þ

Using the coefficients pk ðk ¼ 0; . . . ; 2M þ 1Þ of (2.1), we construct the vector p :¼ ðpk ÞLk ¼ 0 , where p2M þ 2 ¼    ¼ pL :¼ 0. By S :¼ ðdkl1 ÞLk;l ¼ 0 we denote the forward shift matrix, where dk is the Kronecker symbol. Lemma 2.1. Let L; M; N 2 N with 2M þ 1r L rN be given. Furthermore let hk 2 R be given by hk ¼ f ðkÞ ¼

2. Classical Prony method

M X

hk ¼ f ðkÞ ¼

M X

rj eioj k ðk ¼ 0; . . . ; 2NÞ;

ð2:4Þ

j ¼ M

with r0 2 R\f0g, rj 2 C\f0g ðj ¼ 1; . . . ; MÞ, o0 ¼ 0 o o1 o    o oM o p, where (1.2) is fulfilled. Then the rectangular Hankel matrix (2.3) has the singular value 0, where ker H ¼ spanfp; Sp; . . . ; SL2M1 pg

ARTICLE IN PRESS D. Potts, M. Tasche / Signal Processing 90 (2010) 1631–1642

and dimðker HÞ ¼ L  2M:

For shortness the proof is omitted. The Prony method is based on the following Lemma 2.2. Let L; M; N 2 N with 2M þ1 r Lr N be given. Let (2.4) be exact sampled data with r0 2 R\f0g, rj 2 C\f0g ðj ¼ 1; . . . ; MÞ and o0 ¼ 0 o o1 o    o oM o p, where (1.2) is fulfilled. Then the following assertions are equivalent:

This follows immediately from the fact PðzÞ ¼ detðzI  UÞ. Note that we consider a rectangular Hankel matrix (2.3) with only L columns in order to determine the zeros of a polynomial (2.5) with relatively low degree L (see step 2 of Algorithm 2.3). Remark 2.5. Let N 4 2M þ 1. If one knows M or a good approximation of M, then one can use the following least squares Prony method, see e.g. [6]. Since the leading coefficient p2M þ 1 of the characteristic polynomial P0 is equal to 1, (2.2) gives rise to the overdetermined linear system 2M X

(i) The polynomial PðzÞ ¼

L X

¼  h2M þ 1 þ m

uk zk

ðz 2 CÞ;

ð2:5Þ

with real coefficients uk ðk ¼ 0; . . . ; LÞ has 2M þ1 different zeros eioj ðj ¼  M; . . . ; MÞ on the unit circle. (ii) 0 is a singular value of the real rectangular Hankel matrix (2.3) with a right singular vector u :¼ ðul ÞLl ¼ 0 2 RL þ 1 .

The simple proof is omitted. Now we formulate Lemma 2.2 as algorithm: Algorithm 2.3 (Classical Prony method). Input: L; N 2 N ðN b1; 3 r Lr N, L is an upper bound of the number of exponentials), hk ¼ f ðkÞ (k ¼ 0; . . . ; 2N), 0 o e 5 1. u ¼ ðul ÞLl ¼ 0

1. Compute a right singular vector corresponding to the singular value 0 of the exact rectangular Hankel matrix (2.3). 2. Form the corresponding polynomial (2.5) and evaluate ~ with oj 2 ½0; pÞ and (1.2) all zeros eioj ðj ¼ 0; . . . ; MÞ ~ þ 1. lying on the unit circle. Note that LZ 2M ~ 3. Compute r0 2 R and rj 2 C ðj ¼ 1; . . . ; MÞ with (1.2) as least squares solution of the overdetermined linear Vandermonde-type system

rj eioj k ¼ hk ðk ¼ 0; . . . ; 2NÞ:

ð2:6Þ

~ j ¼ M

~ with jr jr e 4. Delete all the pairs ðol ; rl Þ ðl 2 f1; . . . ; MgÞ l and denote the remaining set by fðoj ; rj Þ : j ¼ 1; . . . ; Mg ~ with M r M. Output: M 2 N, r0 2 R, rj 2 C, oj 2 ð0; pÞ ðj ¼ 1; . . . ; MÞ. Remark 2.4. We can determine all roots of the polynomial (2.5) with uL ¼ 1 by computing the eigenvalues of the companion matrix 0

0 B1 B B U :¼ B B0 B^ @ 0

0

...

0

0

...

0

1

...

0

...

1

^ 0

hl þ m pl ¼  h2M þ 1 þ m p2M þ 1

l¼0

k¼0

~ M X

1633

^

u0

1

u1 C C C LL u2 C C2R : ^ C A uL1

ðm ¼ 0; . . . ; 2N  2M  1Þ;

which can be solved by a least squares method. See also the relation to the classic Yule–Walker system [7]. 3. Approximate Prony method The classical Prony method is known to perform poorly when noisy data are given. Therefore numerous modifications were attempted to improve the numerical behavior of the classical Prony method. In practice, only perturbed values h~ k :¼ hk þ ek (k ¼ 0; . . . ; 2N) of the exact sampled data hk ¼ f ðkÞ are known. Here we assume that jek j r e1 with certain accuracy e1 4 0. Then the rectangular error Hankel matrix E :¼ ðek þ l Þ2NL;L k;l ¼ 0 has a small spectral norm satisfying pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi JEJ2 r JEJ1 JEJ1 r ðL þ 1Þð2N  Lþ 1Þe1 r ðN þ 1Þe1 : ð3:1Þ Thus the perturbed Hankel matrix can be represented by ~ :¼ ðh~ k þ l Þ2NL;L ¼ H þ E 2 Rð2NL þ 1ÞðL þ 1Þ : H k;l ¼ 0

ð3:2Þ

Using an analog of the Theorem of Weyl (see [8, p. 419]) and Lemma 2.1, we receive bounds for small singular ~ More precisely, L  2M singular values of H ~ values of H. are contained in ½0; JEJ2 , if each nonzero singular value of H is larger than 2JEJ2 . In the following we use this property and evaluate a small singular value s~ ð0 o s~ r e2 51Þ and corresponding right and left singular ~ vectors of the perturbed rectangular Hankel matrix H. Recently, a very interesting approach is described by Beylkin and Monzo´n [1], where the more general problem of approximation by exponential sums is considered. In contrast to [1], we consider a perturbed rectangular Hankel matrix (3.2). ´n [1]). Let s~ 2 Theorem 3.1 (cf. Beylkin and Monzo ð0; e2  ð0 o e2 51Þ be a small singular value of the perturbed rectangular Hankel matrix (3.2) with a right singular vector u~ ¼ ðu~ k ÞLk ¼ 0 2 RL þ 1 and a left singular vector v~ ¼ 2NL þ 1 . Assume that the polynomial ðv~ k Þ2NL k¼0 2 R ~ PðzÞ ¼

L X k¼0

u~ k zk

ðz 2 CÞ;

ð3:3Þ

ARTICLE IN PRESS 1634

D. Potts, M. Tasche / Signal Processing 90 (2010) 1631–1642

~ n 2 C ðn ¼ 1; . . . ; LÞ. Further let has L pairwise distinct zeros w K 4 2N. Then there exists a unique vector ðan ÞLn ¼ 1 2 CL such that h~ k ¼

~ kn þ s~ dk an w

ðk ¼ 0; . . . ; 2NÞ;

K where the vector ðdk ÞK1 k ¼ 0 2 C is defined as follows. Let L1 X

L X

xk þ n u~ n ¼ 0

u~ k e2pikl=K

ðl ¼ 0; . . . ; K  1Þ;

ð3:4Þ

L X

xk ¼ 2NL X

d^ l :¼

v~ k e2pikl=K

ðl ¼ 0; . . . ; K  1Þ;

if u^ l a0; if u^ l ¼ 0:

v^ l =u^ l 1

ð3:6Þ

1 1 KX d^ e2pikl=K Kl¼0 l

ðk ¼ 0; . . . ; K  1Þ:

ð3:7Þ

The vector ðan ÞLn ¼ 1 can be computed as solution of the linear Vandermonde system L X

ðk ¼ 0; 1; . . .Þ

ð3:5Þ

Then dk :¼

~ kn þ s~ dk an w

n¼1

k¼0

(

ð3:11Þ

and a special solution s~ ðdk Þ1 k ¼ 0 of (3.10). The characteri~ stic equation of (3.10) is PðzÞ ¼ 0 and has the pairwise ~ n 2 C ðn ¼ 1; . . . ; LÞ by assumption. different solutions w Thus the general solution of (3.10) reads as follows:

k¼0

v^ l :¼

ðk ¼ 0; 1; . . .Þ

n¼0

L X n¼1

u^ l :¼

linear difference equation

~ kn ¼ h~ k  s~ dk an w

ðk ¼ 0; . . . ; L  1Þ:

ð3:8Þ

with arbitrary constants an 2 C. By the initial conditions xk ¼ h~ k ðk ¼ 0; . . . ; L  1Þ, the constants an are the unique solutions of the linear Vandermonde system (3.8). Note that ðxk Þ1 k ¼ 0 is not an K-periodic sequence in general. 3. By u~ L a0, we can extend the given data vector ðh~ l Þ2N l¼0 to a sequence ðh~ l Þ1 l ¼ 0 in the following way: L1 1 X s~ ~ h~ 2N þ k :¼  h~ 2NL þ k þ n u~ n þ v u~ L 2NL þ k u~ L

ðk ¼ 1; 2; . . .Þ:

n¼0

Thus we receive L X

h~ k þ n u~ n ¼ s~ v~ k

ðk ¼ 0; 1; . . .Þ;

n¼1

n¼0

Furthermore, if

i.e., the sequence ðh~ l Þ1 l ¼ 0 is also a solution of (3.10) with the initial conditions xk ¼ h~ k (k ¼ 0; . . . ; L  1). Since this initial value problem is uniquely solvable, it follows for all k ¼ 0; 1 . . . that

jv^ l j r g ju^ l j

ðl ¼ 0; . . . ; K  1Þ

for certain constant g 40, then K 1 X

jdk j2 r g2 ;

h~ k ¼ xk ¼

k¼0

L X

~ kn þ s~ dk : an w

n¼1

    L X  ~ kn  r e1 þ ge2 an w hk   

~ u~ ¼ s~ v~ (and H ~ T v~ ¼ s~ u) ~ it follows that Proof. 1. From H

4. By means of discrete Fourier transform of length K (DFT(K)), we can construct a special K-periodic solution s~ ðdk Þ1 k ¼ 0 of (3.10). Then for the K-periodic sequences ~ 1 ~ 1 ðdk Þ1 k ¼ 0 , ðu k Þk ¼ 0 , and ðv k Þk ¼ 0 we obtain

L X

L X

ðk ¼ 0; . . . ; 2NÞ:

n¼1

h~ l þ m u~ l ¼ s~ v~ m

ðm ¼ 0; . . . ; 2N  LÞ:

ð3:9Þ

l¼0

n¼0

We set u~ n :¼ 0 ðn ¼ Lþ 1; . . . ; K  1Þ and extend the vector ~ 1 ~ ðu~ n ÞK1 n ¼ 0 to the K-periodic sequence ðu n Þn ¼ 0 by u n þ jK :¼ u~ n ðj 2 N; n ¼ 0; . . . ; K  1Þ. Analogously, we set v~ n :¼ 0 ðn ¼ 2N  Lþ 1; . . . ; K  1Þ and form the K-periodic se~ ~ quence ðv~ n Þ1 n ¼ 0 by v n þ jK :¼ v n (j 2 N; n ¼ 0; . . . ; K  1). Then we consider the inhomogeneous linear difference equation with constant coefficients L X

dk þ n u~ n ¼

xk þ n u~ n ¼ s~ v~ k

ðk ¼ 0; 1; . . .Þ

ð3:10Þ

n¼0

under the initial conditions xk ¼ h~ k ðk ¼ 0; . . . ; L  1Þ. By assumption, the polynomial P~ possesses L pairwise different zeros such that u~ L a0. Hence the linear difference equation (3.10) has the order L. Thus the above initial value problem of (3.10) is uniquely solvable. 2. As known, the general solution of (3.10) is the sum of the general solution of the corresponding homogeneous

K 1 X

dk þ n u~ n ¼ v~ k

ðk ¼ 0; . . . ; K  1Þ;

ð3:12Þ

n¼0

where !1

K 1 X

dk þ n u~ n

n¼0

k¼0

is the K-periodic correlation of ðu~ k Þ1 k ¼ 0 with respect to . Introducing the following entries (3.4), (3.5), and ðdk Þ1 k¼0 d^ l :¼

K 1 X

dk e2pikl=K

ðl ¼ 0; . . . ; K  1Þ;

k¼0

from (3.12) it follows by DFT(K) that d^ l u^ l ¼ v^ l

ðl ¼ 0; . . . ; K  1Þ:

Note that u^ l ¼ 0 implies v^ l ¼ 0. Hence we obtain (3.6) with jd^ l j r g ðl ¼ 0; . . . ; K  1Þ. In the case L ¼ N, one can show that g ¼ 1. Using inverse DFT(K), we receive (3.7). By

ARTICLE IN PRESS D. Potts, M. Tasche / Signal Processing 90 (2010) 1631–1642

but rather in approximate representations     M X   ioj k  h~ k  r e j  r e~ ðk ¼ 0; . . . ; 2NÞ    j ¼ M

Parseval equation we see that K 1 X

jdk j2 ¼

k¼0

1 1 KX jd^ j2 r g2 Kl¼0 l

and hence jdk j r g ðk ¼ 0; . . . ; K  1Þ. Then for k ¼ 0; . . . ; K  1 we can estimate     L X  k ~ an w n  hk    n¼1 r jhk  h~ k j þjh~ k 

L X n¼1

&

Remark 3.2. This Theorem 3.1 yields a different repre~ n and s~ remain sentation for each K 4 2N even though w the same. If K is chosen as power of 2, then the entries u^ l , v^ l and dk ðl; k ¼ 0; . . . ; K  1Þ can be computed by fast Fourier transforms. Note that the least squares solution ðb~ n ÞLn ¼ 1 of the overdetermined linear Vandermonde-type system L X

~ kn ¼ h~ k b~ n w

ðk ¼ 0; . . . ; 2NÞ

ð3:13Þ

n¼1

has an error with Euclidean norm less than ge2 , since 2 2     2N  L 2N  L X X X X ~ ~ k k ~ ~ ~ h b h  r  a w w  k  k n n n n     n¼1 n¼1 k¼0

¼

k¼0

2N X

2

js~ dk j r s~

k¼0

2

K1 X

jdk j2 r ðge2 Þ2 :

k¼0

By Remark 3.2 an algorithm of the first APM (cf. [1]) reads as follows: Algorithm 3.3 (APM 1). Input: L, N 2 N (3r L rN, L is an upper bound of the number of exponentials), h~ k 2 R ðk ¼ 0; . . . ; 2NÞ, accuracies e1 , e2 4 0. 1. Compute a small singular value s~ 2 ð0; e2  and corresponding right and left singular vectors u~ ¼ 2NL þ 1 of the perturbed ðu~ l ÞLl ¼ 0 2 RL þ 1 , v~ ¼ ðv~ l Þ2NL l¼0 2 R rectangular Hankel matrix (3.2). ~ n 2 C (n ¼ 1; . . . ; L) of the corre2. Determine all zeros w sponding polynomial (3.3). Assume that all the zeros of P~ are simple. 3. Determine the least squares solution ðb~ n ÞLn ¼ 1 2 CL of the overdetermined linear Vandermonde-type system (3.13). 4. Denote by ðwj ; rj Þ ðj ¼  M; . . . ; MÞ all that pairs ~ k j  1 and ~ k ; b~ k Þ ðk ¼ 1; . . . ; LÞ with the properties jw ðw jb~ k j Z e1 . Assume that arg wj 2 ½0; pÞ and wj ¼ w j ðj ¼ 0; . . . ; MÞ. Set oj :¼ arg wj 4 0 ðj ¼ 1; . . . ; MÞ and o0 :¼ 0. Output: M 2 N, rj 2 C, oj 2 ½0; pÞ ðj ¼ 0; . . . ; MÞ. Similarly as in [1], we are not interested in exact representations of the sampled values h~ k ¼

L X n¼1

~ kn b~ n w

ðk ¼ 0; . . . ; 2NÞ

for very small accuracy e~ 4 0 and minimal number M of nontrivial summands. This fact explains the name of our method too.

4. Approximate Prony method based on NFFT

~ kn j r e1 þ gs~ r e1 þ ge2 : an w

This completes the proof.

1635

Now we present a second new approximate Prony method by means of matrix perturbation theory. First we introduce the rectangular Vandermonde-type matrix 2 Cð2N þ 1Þð2M þ 1Þ : V :¼ ðeikoj Þ2N;M k ¼ 0;j ¼ M

ð4:1Þ

Note that V is also a nonequispaced Fourier matrix (see [9,10,4]). We discuss the properties of V. Especially, we show that V is left-invertible and estimate the spectral norm of its left inverse. For these results, the separation distance (1.3) plays a crucial role. Lemma 4.1. Let M, N 2 N with N Z 2M þ 1 given. Furthermore let o0 ¼ 0 o o1 o    o oM o p with a separation distance q4

1 p2 pffiffiffi N þ1 3

ð4:2Þ

be given. Let (1.2) be fulfilled. Let D :¼ diagð1  jkj= ðN þ1ÞÞN k ¼ N be a diagonal matrix. Then for arbitrary r 2 C2M þ 1 , the following inequality: 

N þ1 

p4 3q2 ðN þ 1Þ

  JrJ2 r JD1=2 VrJ22 r N þ 1 þ

p4

 JrJ2

3q2 ðN þ 1Þ

is fulfilled. Further L :¼ ðVH DVÞ1 VH D

ð4:3Þ

is a left inverse of (4.1) and the squared spectral norm of L can be estimated by JLJ22 r

3q2 ðN þ 1Þ 3q2 ðN þ1Þ2  p4

:

ð4:4Þ

Proof. 1. The rectangular Vandermonde-type matrix V 2 Cð2N þ 1Þð2M þ 1Þ has full rank 2M þ 1, since its submatrix ðeikoj Þ2M;M k ¼ 0;j ¼ M is a regular Vandermonde matrix. Hence we infer that VH DV is Hermitian and positive definite such that all eigenvalues of VH DV are positive. 2. We introduce the 2pperiodic Feje´r kernel FN by   N X jkj FN ðxÞ :¼ 1 eikx : Nþ1 k ¼ N Then we obtain ðVH DVÞj;l ¼ eiNoj FN ðoj  ol ÞeiNol

ðj; l ¼  M; . . . ; MÞ

for the ðj; lÞth entry of the matrix VH DV. We use Gershgorin’s Disk Theorem (see [8, p. 344]) such that for an arbitrary eigenvalue l of the matrix VH DV we preserve

ARTICLE IN PRESS 1636

D. Potts, M. Tasche / Signal Processing 90 (2010) 1631–1642

obtain for the squared spectral norm of the left inverse L (see also [4, Theorem 4.2])

the estimate (see also [11, Theorem 4.1]) jl  FN ð0Þj ¼ jl  N  1 9 8 > > = > ; :j ¼ M

JLJ22 r JðVH DVÞ1 VH D1=2 J22 JD1=2 J22 r lmin r

ð4:5Þ

jal

This completes the proof.

3. Now we estimate the right-hand side of (4.5). As known, the Feje´r kernel FN can be written in the form   1 sinððN þ 1Þx=2Þ 2 0 rFN ðxÞ ¼ : N þ1 sinðx=2Þ Thus we obtain the estimate

N4

M 1 X jFN ðoj  ol Þj r jsinððoj  ol Þ=2Þj2 N þ 1 j ¼ M j ¼ M

 p4

:

&

p2 q

;

ð4:7Þ

the squared spectral norm of (4.3) satisfies

jal

for l 2 fM; . . . ; Mg. In the case l ¼ 0, we use (1.2), qM o p, oj Z jq ðj ¼ 1; . . . ; MÞ and sinx Z

2

p

x

jsinðoj =2Þj2 ¼ 2

M X

ðsinðoj =2ÞÞ2 r 2p2

j¼1

j ¼ M ja0

M X

o2 j r

j¼1

jsinððoj  ol Þ=2Þj2 o

j ¼ M jal

M 2p2 X p4 j2 o 2 : q2 j ¼ 1 3q

p4

Hence it follows in each case that

p4

:

ð4:6Þ

4. Let lmin and lmax be the smallest and largest eigenvalue of VH DV, respectively. Using (4.6), we receive 4

Nþ1 

p

3q2 ðN þ1Þ

r N þ 1þ

r lmin rN þ 1 r lmax

p4 3q2 ðN þ1Þ

;

where by assumption 4

Nþ1 

p

3q2 ðN þ1Þ

p2

p2

4 pffiffiffi : N 3ðN þ 1Þ

Thus the assumptions of Lemma 4.1 are fulfilled. Using the upper estimate (4.4) of JLJ22 and (4.7), we obtain the result. &

3q2

jsinððoj  ol Þ=2Þj ¼ jsinð 7 p þðoj  ol Þ=2Þj:

3q2 ðN þ 1Þ

2p p2 þ1o o N: q q

Further from (4.7) it follows that q4

in case l 2 f 7 1; . . . ; 7 Mg. Note that we use the 2pperiodization of the frequencies oj ðj ¼  M; . . . ; MÞ such that

jl  N  1j o

3 : 2N þ2

2M þ 1o

By similar arguments, we obtain M X

JLJ22 r

Proof. By (4.7) and q M o p, we see immediately that

ðx 2 ½0; p=2Þ

and then we estimate the above sum by M X

3q2 ðN þ 1Þ2

Corollary 4.2. Let M 2 N be given. Furthermore let o0 ¼ 0 o o1 o    o oM o p with separation distance q. Let (1.2) be fulfilled. Then for each N 2 N with

M X jal

3q2 ðN þ 1Þ

1

4 0:

Using the variational characterization of the Rayleigh–Ritz ratio for the Hermitian matrix VH DV (see [8, p. 176]), we obtain for arbitrary r 2 C2M þ 1 that

lmin JrJ2 r JD1=2 VrJ2 r lmax JrJ2 :

Corollary 4.3. Let M 2 N be given. Furthermore let o0 ¼ 0 o o1 o    o oM o p with separation distance q. Let (1.2) be fulfilled. Then the squared spectral norm of the Vandermonde-type matrix V is bounded, i.e.,   2p p 1 þ ln : JVJ22 r 2N þ1 þ q q Proof. We follow the lines of the proof of Lemma 4.1 with the trivial diagonal matrix D :¼ diagð1ÞN k ¼ N . Instead of the Feje´r kernel FN we use the modified Dirichlet kernel DN ðxÞ :¼

2N X

eikx ¼ eiNx

k¼0

sinðð2N þ1Þx=2Þ sinðx=2Þ

and we obtain for the ðj; lÞth entry of VH V ðVH VÞj;l ¼ DN ðoj  ol Þ

ðj; l ¼  M; . . . ; MÞ:

We proceed with jl  DN ð0Þj ¼ jl  ð2N þ1Þj 8 9 > > M > :j ¼ M ; jal

From LV ¼ ðVH DVÞ1 VH D1=2 D1=2 V ¼ I it follows that (4.3) is a left inverse of V and that the pffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffi singular values of D1=2 V lie in ½ lmin ; lmax . Hence we

use the estimate M X j ¼ M jal

jDN ðoj  ol Þj r

M X j ¼ M jal

jsinððoj  ol Þ=2Þj1

ARTICLE IN PRESS D. Potts, M. Tasche / Signal Processing 90 (2010) 1631–1642

and infer M X

1637

and hence

jsinðoj =2Þj1 ¼ 2

M X

ðsinðoj =2ÞÞ1 r2p

j¼1

j ¼ M ja0

M X

o1 j

~ 2 þ2 r20 Pð1Þ

j¼1

~ ioj Þj2 r ðjs~ j þ JEJ2 Þ2 JL1 J2 : jrj j2 jPðe 2

j¼1

  M 2p X 2p 2p p ð1 þ lnMÞ o 1 þ ln : r j1 o q j¼1 q q q

This completes the proof.

&

Lemma 4.5. If the assumptions of Theorem 4.4 are fulfilled with sufficiently small accuracies e1 , e2 40 and if d 4 0 is the smallest singular value a0 of H, then

Finally we obtain for the largest eigenvalue of VH V   2p p 1 þln lmax ðVH VÞ r2N þ 1 þ q q and hence the assertion.

M X

~ r Ju~  P uJ

&

Let L, M, N 2 N with 2M þ 1r L rN be given, where L is even. Now we consider a special submatrix V1 of V defined by 2 Cð2NL þ 1Þð2M þ 1Þ : V1 :¼ ðeikoj Þ2NL;M k ¼ 0;j ¼ M

Using Lemma 4.1, the matrix V1 has a left inverse 1 H V1 D1 : L1 :¼ ðVH 1 D1 V1 Þ

d

;

ð4:9Þ

where P is the orthogonal projector of RL þ 1 onto ker H. Furthermore the polynomial (3.3) has zeros close to eioj ðj ¼  M; . . . ; MÞ, where ~ 2 þ2 Pð1Þ

The matrix V1 has full rank 2M þ 1 (see step 1 of Lemma 4.1). Hence VH 1 D1 V1 is Hermitian and positive definite, where  NL=2 jkj : D1 :¼ diag 1  N þ 1  L=2 k ¼ N þ L=2

e2 þðN þ 1Þ e1

    e2 þ ðN þ1Þe1 2 ~ ioj Þj2 r 2N  Lþ1 þ 2p 1 þln p jPðe : q q d j¼1 M X

Proof. 1. Let u~ and v~ be normed right and left singular ~ with respect to the singular value s~ 2 ð0; e2  vectors of H ~ u~ ¼ s~ v. ~ By assumption, d 4 0 is the smallest such that H 2 singular value a0 of the exact Hankel matrix H. Then d is the smallest eigenvalue a0 of the symmetric matrix HT H. Using the Rayleigh–Ritz Theorem (see [8, pp. 176–178]), we receive uT HT Hu uT u u?ker H

d2 ¼ min ua0 Theorem 4.4. Let L, M, N 2 N with 2M þ1 r Lr N be given, where L is even. Furthermore let hk 2 R be given as in (2.4) with r0 2 R\f0g, rj 2 C\f0g ðj ¼ 1; . . . ; M), and o0 ¼ 0 o o1 o    o oM o p. Let (1.2) be fulfilled. Assume that the perturbed Hankel matrix (3.2) has s~ 2 ð0; e2  as singular value with the corresponding normed right and left singular 2NL þ 1 . vectors u~ ¼ ðu~ n ÞLn ¼ 0 2 RL þ 1 and v~ ¼ ðv~ m Þ2NL m¼0 2 R ~ Let P~ be the polynomial (3.3) related to u. ~ ioj Þ ðj ¼  M; . . . ; MÞ fulfill the estimate Then the values Pðe ~ 2 þ2 r20 Pð1Þ

M X

ðm ¼ 0; . . . ; 2N  LÞ:

l¼0

Applying (2.4) and h~ k ¼ hk þ ek , we obtain for m ¼ 0; . . . ; 2N  L M X j ¼ M

~ ioj Þeimoj ¼ s~ v~ m  rj Pðe

L X

el þ m u~ l :

ð4:8Þ

l¼0

2. Using the matrix-vector notation of (4.8) with the rectangular Vandermonde-type matrix V1 , we receive ~ ioj ÞÞM ~ ~ ~ V1 ðrj Pðe j ¼ M ¼ s v  E u: 1 H The matrix V1 has a left inverse L1 ¼ ðVH V 1 D1 1 D1 V 1 Þ such that

~ ioj ÞÞM ~ ~ ~ ðrj Pðe j ¼ M ¼ s L 1 v  L 1 E u

Thus the following estimate:

dJu~  uJ r JHðu~  uÞJ is true for all u 2 RL þ 1 with u~  u ? ker H. Especially for ~ we see that u~  P u~ ? ker H and hence by (3.1) u ¼ P u, r s~ þ ðN þ1Þe1

~ u~ ¼ s~ v, ~ i.e., Proof. 1. By assumption we have H h~ l þ m u~ l ¼ s~ v~ m

JHuJ JHðu~  uÞJ ¼ min : u~ ua0 JuJ Ju~  uJ u?ker H u~ u?ker H

d ¼ min ua0

~  EÞ uJ ~ rJH uJ ~ ¼ JðH ~ ¼ Js~ v~  E uJ ~ dJu~  P uJ

~ ioj Þj2 r ðe2 þ JEJ2 Þ2 JL1 J2 : jrj j2 jPðe 2

j¼1

L X

and hence

such that (4.9) follows. 2. Thereby u ¼ P u~ is a right singular vector of H with respect to the singular value 0. Thus the corresponding polynomial (2.5) has the values eioj ðj ¼  M; . . . ; MÞ as zeros by Lemma 2.2. By (4.9), the coefficients of P differ ~ Consequently, the only a little from the coefficients of P. zeros of P~ lie nearby the zeros of P, i.e., P~ has zeros close to eioj ðj ¼  M; . . . ; MÞ (see [8, pp. 539–540]). By JV1 J2 ¼ JVH 1 J2 , (4.9), and Corollary 4.3, we obtain the estimate M X j ¼ M

~ ioj Þj2 ¼ jPðe

M X

~ ioj Þj2 jPðeioj Þ  Pðe

j ¼ M

H 2 ~ 2 ~ 2 ¼ JVH 1 ðu  uÞJ r JV1 J2 Ju  uJ    2p p e2 þðN þ 1Þe1 2 ð1 þ ln Þ r 2N  L þ1 þ : q q d

This completes the proof.

&

ARTICLE IN PRESS 1638

D. Potts, M. Tasche / Signal Processing 90 (2010) 1631–1642

Corollary 4.6. Let L, M 2 N with LZ 2M þ 1 be given, where L is even. Furthermore let o0 ¼ 0 o o1 o    o oM o p with separation distance q be given. Let (1.2) be fulfilled. Further let N 2 N with  2  p L þ ;L N 4 max q 2 be given. Assume that the perturbed rectangular Hankel matrix (3.2) has a singular value s~ 2 ð0; e2  with corresponding normed right and left singular vectors u~ ¼ ðu~ n ÞLn ¼ 0 2 þ1 RL þ 1 and v~ ¼ ðv~ n Þ2NL 2 R2NL þ 1 . n¼0 ~ ioj Þj ðj ¼  M; . . . ; MÞ of (3.3) can be Then the values jPðe estimated by ~ 2 þ2 r20 Pð1Þ

M X j¼1

~ ioj Þj2 r jrj j2 jPðe

3 ðe2 þ ðN þ 1Þe1 Þ2 ; 2N  L þ2

~ ioj Þj ðj ¼  M; . . . ; MÞ are small. Further i.e., the values jPðe ~ the polynomial P has zeros close to eioj ðj ¼  M; . . . ; MÞ. Proof. The assertion follows immediately by the assumption on N, the estimate (3.1), Corollary 4.2, and Theorem 4.4. & Now we can formulate a second approximate Prony method. Algorithm 4.7 (APM 2). Input: L, N 2 N ð3 r Lr N, L is an upper bound of the number of exponentials), h~ k ¼ f ðkÞ þ ek ðk ¼ 0; . . . ; 2NÞ with jek j r e1 , accuracies e1 , e2 , e3 , e4 . ð1Þ 1. Compute a right singular vector u~ ¼ ðu~ ð1Þ ÞLl ¼ 0 correl sponding to a singular value sð1Þ 2 ð0; e2  of the perturbed rectangular Hankel matrix (3.2). P ð1Þ 2. Form the corresponding polynomial P~ ðzÞ :¼ L k¼0

ð1Þ

io zk and evaluate all zeros rjð1Þ e j ðj ¼ 1; . . . ; Mð1Þ Þ u~ ð1Þ k

2 ð0; pÞ, (1.2) and jrjð1Þ  1j r e4 . Note that with oð1Þ j L Z2Mð1Þ þ 1. ð2Þ ÞLl ¼ 0 corre3. Compute a right singular vector u~ ¼ ðu~ ð2Þ l ð2Þ sponding to a singular value s 2 ð0; e2  ðsð1Þ asð2Þ Þ of the perturbed rectangular Hankel matrix (3.2). P ð2Þ 4. Form the corresponding polynomial P~ ðzÞ :¼ L k¼0

ð2Þ

zk and evaluate all zeros rkð2Þ eiok ðk ¼ 1; . . . ; Mð2Þ Þ u~ ð2Þ k 2 ð0; pÞ, (1.2) and jrkð2Þ  1j r e4 . Note that with oð2Þ k ð2Þ L Z2M þ 1. 5. Determine all frequencies ~ o~ l :¼ 12ðoð1Þ þ oð2Þ Þ ðl ¼ 1 . . . ; MÞ; jðlÞ kðlÞ if there exist indices jðlÞ 2 f1; . . . ; M ð1Þ g and kðlÞ 2 ð1Þ  oð2Þ j r e3 . Replace rjðlÞ f1; . . . ; Mð2Þ g such that joð1Þ jðlÞ kðlÞ ð2Þ ~ þ 1. and rkðlÞ by 1. Note that LZ 2M ~ with (1.2) as ~ ~ 6. Compute r 0 2 R and r j 2 C ðj ¼ 1; . . . ; MÞ least squares solution of the overdetermined linear Vandermonde-type system ~ M X ~ j ¼ M

r~ j eio~ j k ¼ h~ k ðk ¼ 0; . . . ; 2NÞ

ð4:10Þ

with the diagonal preconditioner D ¼ diagð1  jkj= ðN þ 1ÞÞN k ¼ N . For very large M and N use the CGNR method (conjugate gradient on the normal equations, see e.g. [12, p. 288]), where the multiplication of the ~ M Vandermonde-type matrix V~ :¼ ðeiko~ j Þ2N; ~ is reak ¼ 0;j ¼ M lized in each iteration step by the nonequispaced fast Fourier transform NFFT, see [9,10]. ~ with jr~ j r e1 ~ l ; r~ j Þ ðl 2 f1; . . . ; MgÞ 7. Delete all the pairs ðo l and denote the remaining frequency set by foj : ~ j ¼ 1; . . . ; Mg with M r M. 8. Repeat step 6 and solve the overdetermined linear PM ioj k ¼ h~ k ðk ¼ 0; Vandermonde-type system j ¼ M rj e . . . ; 2NÞ with respect to the new frequency set foj : j ¼ 1; . . . ; Mg again. Output: M 2 N, r0 2 R, rj 2 C ðrj ¼ r j Þ, oj 2 ð0; pÞ ðj ¼ 1; . . . ; MÞ. ð1Þ

We can determine the zeros of the polynomials P~ and ð2Þ by computing the eigenvalues of the related compaP~ nion matrices, see Remark 2.4. Note that by Corollary 4.6, the correct frequencies oj can be approximately computed by any right singular vector corresponding to a small singular value of (3.2). In general, the additional roots are not nearby for different right singular vectors. In many applications, Algorithm 4.7 can be essentially simplified by using only one right singular vector (cf. steps 1–2), omitting the steps 3–5, and determining the frequencies, the coefficients and the number M in the steps 6–8. We use this simplified variant of Algorithm 4.7 in our numerical examples of Section 6. Note that more advanced methods for computing the unknown frequencies were developed, such as matrix pencil methods. In [13,14] a relationship between the matrix pencil methods and several variants of the ESPRIT method [15,16] is derived showing comparable performance. Now we replace the steps 1–5 of Algorithm 4.7 by the least squares ESPRIT method [17, p. 493]. This leads to the following Algorithm 4.8 (APM based on ESPRIT). Input: L, N, P 2 N ð3 r Lr N, L is an upper bound of the number of exponentials, P þ 1 is upper bound of the number of positive frequencies, P þ 1 rLÞ, h~ k ¼ f ðkÞ þek ðk ¼ 0; . . . ; 2NÞ with jek jr e1 . 1. Form the perturbed rectangular Hankel matrix (3.2). 2. Compute the singular value decomposition of (3.2) with a diagonal matrix S 2 Rð2NL þ 1ÞðL þ 1Þ with nonnegative diagonal elements in decreasing order, and unitary matrices L 2 Cð2NL þ 1Þð2NL þ 1Þ and U :¼ ~ ¼ LSUH . ðuk;l ÞLk;l ¼ 0 2 CðL þ 1ÞðL þ 1Þ such that H and U2 :¼ 3. Form the matrices U1 :¼ ðuk;l ÞL1;P k;l ¼ 0 . ðuk þ 1;l ÞL1;P k;l ¼ 0 4. Compute the matrix P :¼ Uy1 U2 2 CðP þ 1ÞðP þ 1Þ , where y 1 H U1 is the Moore–Penrose pseudoinverse U1 :¼ ðUH 1 U1 Þ of U1 . 5. Compute the eigenvalues r~ k eio~ k ðk ¼ 1; . . . ; P þ 1Þ of the ~ k 2 ð0; pÞ and r~ k  1 ðk ¼ 1; . . . ; P þ 1Þ. matrix P, where o 6. Compute r~ 0 2 R and r~ j 2 C ðj ¼ 1; . . . ; P þ 1Þ with (1.2) as least squares solution of the overdetermined linear

ARTICLE IN PRESS D. Potts, M. Tasche / Signal Processing 90 (2010) 1631–1642

Vandermonde-type system PX þ1

r~ j eio~ j k ¼ h~ k ðk ¼ 0; . . . ; 2NÞ;

ð4:11Þ

j ¼ P1

with the diagonal preconditioner D ¼ diagð1  jkj= ðN þ1ÞÞN k ¼ N . For large P and N, use the CGNR method, where the multiplication of the Vandermonde-type þ1 is realized in each iteration matrix ðeiko~ j Þk2N;P ¼ 0;j ¼ P1 step by NFFT (see [9,10]). ~ l ; r~ j Þ ðl 2 f1; . . . ; P þ 1gÞ with 7. Delete all the pairs ðo jr~ l j r e1 and denote the remaining frequency set by foj : j ¼ 1; . . . ; Mg with M rP þ 1. 8. Repeat step 6 and solve the overdetermined linear PM io j k ¼ h~ k ðk ¼ 0; Vandermonde-type system j ¼ M rj e . . . ; 2NÞ with respect to the new frequency set foj : j ¼ 1; . . . ; Mg again. Output: M 2 N, r0 2 R, rj 2 C ðrj ¼ r j Þ, oj 2 ð0; pÞ ðj ¼ 1; . . . ; MÞ. Since a good estimate for the number of positive frequencies is often unknown in advance, we can simply choose P ¼ L  1 in our numerical experiments. Clearly, we can choose the parameter L similarly as suggested in [14] based on the accuracy e2 as in Algorithm 4.7. Furthermore we would like to point out that in Algorithm ~ related to the small 4.8 the right singular vectors of H singular values are discarded, otherwise the Algorithm 4.7 uses right singular vectors related to the small singular values in order to compute the roots of the corresponding polynomials. Note that Algorithm 4.8 requires more arithmetical operations than Algorithm 4.7, due to the evaluation of the Moore–Penrose pseudoinverse of U1 in step 4 of Algorithm 4.8. 5. Sensitivity analysis of the approximate Prony method

~ where (4.3) is a left inverse of V. If the and q~ ¼ L h, assumptions of Corollary 4.2 are fulfilled, then for each N 2 N with N 4 p2 q1 the condition number kðVÞ :¼ JLJ2 JVJ2 is uniformly (with respect to N) bounded by sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi 3 p kðVÞ r 3 þ 1 þ ln : ð5:2Þ p q Furthermore, the following stability inequalities are fulfilled: rffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi 3 ~ r 3e ; ð5:3Þ Jh  hJ Jq  q~ J r 1 2N þ 2 ~ Jq  q~ J Jh  hJ r kðVÞ : JqJ JhJ

This provides (5.2), since by assumption ðN þ 1Þq 4 p2 . 2. The inequality (5.3) follows immediately from Jq  ~ and Corollary 4.2. The estimate (5.4) q~ J rJLJ2 Jh  hJ arises from (5.3) by multiplication with JVqJ JhJ1 ¼ 1. & Now we consider the general case, where the weighted least squares problem (5.1) is solved for perturbed data ~ j ðj ¼  M; h~ k ðk ¼ 0; . . . ; 2NÞ and computed frequencies o . . . ; MÞ. Theorem 5.2. Assume that jh~ k  hk j r e1 ðk ¼ 0; . . . ; 2NÞ ~ j  oj j r d ðj ¼ 0; . . . ; MÞ. Let V be given by (4.1) and jo and let : V~ :¼ ðeiko~ j Þ2N;M k ¼ 0;j ¼ M

~ ¼ min: JD1=2 ðV~ q~  hÞJ

N 4 p2 maxfq1 ; q~

Here h~ ¼ ðh~ k Þ2N k ¼ 0 is the perturbed data vector and V~ ¼ ðeiko~ j Þ2N;M k ¼ 0;j ¼ M is the Vandermonde-type matrix with the computed ~ 0 oo ~ 1o  ~ j ðj ¼  M; . . . ; MÞ, where 0 ¼ o frequencies o ~j¼ o ~ j ðj ¼ 1; . . . ; MÞ and q~ is the separa~ M o p, o oo ~j: tion distance of the computed frequency set fo ~ M . Note that V~ has ~ M þ 1 ¼ 2p  o j ¼ 0; . . . ; M þ1g with o full rank and that the unique solution of (5.1) is given by H ~ 1 V~ H D of V. ~ q~ ¼ L~ h~ with the left inverse L~ ¼ ðV~ D VÞ We begin with a normwise perturbation result, if all ~ j ðj ¼ frequencies are exactly determined, i.e., oj ¼ o M; . . . ; MÞ, and if a perturbed data vector h~ is given. ~ j ðj ¼  M; . . . ; MÞ and that Lemma 5.1. Assume that oj ¼ o jh~ k  hk j r e1 . Let V be given by (4.1). Further let V q ¼ hao

ð5:4Þ

Proof. 1. The condition number kðVÞ of the rectangular Vandermonde-type matrix V is defined as the number JLJ2 JVJ2 . Note that kðVÞ does not coincide with the condition number of V related to the spectral norm, since the left inverse L is not the Moore–Penrose pseudoinverse of V by DaI. Applying Corollaries 4.2 and 4.3, we receive that sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   3p p 1 þ ln : kðVÞ r 3 þ ðN þ 1Þq q

In this section, we discuss the sensitivity of step 6 of the Algorithms 4.7 and 4.8, respectively. Assume that ~ ¼ M. We solve the overdetermined N Z 2M þ 1 and M ~ as linear Vandermonde-type system (4.10) with M ¼ M weighted least squares problem ð5:1Þ

1639

~ where L~ ¼ ðV~ H D VÞ ~ 1 V~ H D is Further let Vq ¼ h and q~ ¼ L~ h, ~ a left inverse of V. If the assumptions of Corollary 4.2 are fulfilled, then for each N 2 N with 1

g

ð5:5Þ

the following estimate: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi Jq  q~ J r ð6N þ 6Þ ð2M þ 1Þ JhJd þ 3e1 is fulfilled. Proof. 1. Using the matrix-vector notation, we can write (2.6) in the form V q ¼ h. The overdetermined linear ~ Using the left ~ ¼ M reads V~ q~ ¼ h. system (4.10) with M ~ the solution q~ of the weighted least squares inverse V, problem ~ ¼ min JD1=2 ðV~ q~  hÞJ ~ Thus it follows that is q~ ¼ L~ h. ~ ¼ JL~ V~ q  L~ V q  Lð ~ h~  hÞJ Jq  q~ J ¼ Jq  L~ hJ ~ ~ ~ ~ r JLJ JV  VJ JqJ þ JLJ Jh  hJ: 2

2

2

ARTICLE IN PRESS 1640

D. Potts, M. Tasche / Signal Processing 90 (2010) 1631–1642

2. Now we estimate the squared spectral norm of V~  V by means of the Frobenius norm 2N M X X

~ 2¼ ~ 2 rJV  VJ JV  VJ 2 F

jeikoj  eiko~ j j2

k ¼ 0 j ¼ M

¼

M X

2N X

j ¼ M

k¼0

!

Example 6.1. We sample the trigonometric sum

~ j ÞkÞ ½2  2cosððoj  o

f1 ðxÞ :¼ 14  8cosð0:453xÞ þ 9sinð0:453xÞ þ 4cosð0:979xÞ þ8sinð0:979xÞ  2cosð0:981xÞ þ 2cosð1:847xÞ 3sinð1:847xÞ þ0:1cosð2:154xÞ  0:3sinð2:154xÞ ð6:1Þ

  M X ~ j ÞÞ sinðð2N þ1=2Þðoj  o ¼ : 4N þ1  ~ j Þ=2ÞÞ sinððoj  o j ¼ M Using the special property of the Dirichlet kernel !   sinðð4N þ1Þx=2Þ 1 1 3 N 1 þ þ Z 4N þ 1 þ  2N þ x2 ; sinðx=2Þ 3 2 6 24 we infer ~ 2r JV  VJ 2

M X j ¼ M

!   1 1 3 N 1 ð2N þ 1=2Þ3 ð2M þ 1Þ 2 2N þ   d2 r d : 3 2 6 24 3

ð5:6Þ Thus we receive   1 1 3=2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~ ~ 2 Jh~  hJ: Jq  q~ J r pffiffiffi 2N þ 2M þ 1JLJ JqJd þJLJ 2 3 From Corollary 4.2 it follows that for each N 2 N with (5.5) rffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 ~ 2r JLJ : 2N þ 2 Finally we use pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi Jh~  hJ r 2N þ 1Jh~  hJ1 r 2N þ 1e1 ; JqJ ¼ JLhJ r JLJ2 JhJ r

rffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 JhJ 2N þ 2

and obtain the result.

stant points of the interval ½0; 2N, where f~ is the recovered exponential sum. For increasing number of sampled data, we obtain a very precise parameter estimation. With other words, N acts as regularization parameter.

at the equidistant nodes x ¼ k ðk ¼ 0; . . . ; 2NÞ. Then we apply the Algorithms 4.7 and 4.8 with exact sampled data hk ¼ f1 ðkÞ, i.e., ek ¼ 0 ðk ¼ 0; . . . ; 2NÞ. Thus the accuracy e1 can be chosen as the unit roundoff e1 ¼ 253  1:11  1016 (see [20, p. 45]). Note that the frequencies o2 ¼ 0:979 and o3 ¼ 0:981 are very close, so that q ¼ 0:002. Nevertheless for small N and for the accuracies e2 ¼ e3 ¼ e4 ¼ 107 , we obtain very precise results, see Table 1. Here we assume that P is unknown and use the estimate P ¼ L  1 in Algorithm 4.8. Example 6.2. We use again the function (6.1). Now we consider noisy sampled data h~ k ¼ f1 ðkÞ þek ðk ¼ 0; . . . ; 2NÞ, where ek is a random error with jek j r e1 ¼ 5  105 , more precisely we add pseudorandom values drawn from the standard uniform distribution on ½5  105 ; 5  105 . Now we choose e4 o 5  105 . We see that for noisy sampled data increasing N improves the results, see Table 2. The last column shows the signal-to-noise ratio 2N SNR :¼ 10 log10 ðJðf1 ðkÞÞ2N k ¼ 0 J=Jðek Þk ¼ 0 JÞ. Example 6.3. As in [21], we sample the 24-periodic trigonometric sum   px px px 5px f2 ðxÞ :¼ 2cos þ200cos þ2cos þ2cos 6 4 2 6

&

By Theorem 5.2, we see that we have to compute the frequencies very carefully. This is the reason why we repeat the computation of the frequencies in the steps 3–4 of Algorithm 4.7. 6. Numerical examples Finally, we apply the Algorithms 4.7 and 4.8 to various examples. We have implemented our algorithms in MATLAB with IEEE double precision arithmetic. There exist a variety of algorithms to recover the frequencies oj like ESPRIT [15,16] and more advanced methods like total least squares phased averaging ESPRIT [18,19]. We emphasize that the aim of this paper is not a widespread numerical comparison with all the existing methods. We demonstrate that a straightforward application of the ESPRIT method [17, p. 493] leads to similar results. M Let x :¼ ðoj ÞM j ¼ 0 and q :¼ ðrj Þj ¼ 0 . We compute the ~ j j, where o ~ j are frequency error eðxÞ :¼ maxj ¼ 0;...;M joj  o computed by Algorithms 4.7 and 4.8, respectively. Furthermore let eðqÞ :¼ maxj ¼ 0;...;M jrj  r~ j j be the error of the coefficients. Finally we determine the error eðf Þ :¼ maxjf ðxÞ  f~ ðxÞj=maxjf ðxÞj computed on 10 000 equidi-

at the equidistant nodes x ¼ k ðk ¼ 0; . . . ; 2NÞ. In [21], only the dominant frequency o2 ¼ p=4 is iteratively computed via zeros of Szego¨ polynomials. For N ¼ 50, o2 is determined to 2 correct decimal places. For N ¼ 200 and 500, the frequency o2 is computed to 4 correct decimal places in [21]. Now we apply Algorithms 4.7 and 4.8 with the same parameters as in Example 6.1, but we use only the 19 sampled data f2 ðkÞ ðk ¼ 0; . . . ; 18Þ. In this case we are able to compute all frequencies and all coefficients very accurate, see Table 3. Table 1 Results for f1 from Example 6.1. N

L

eðxÞ

eðqÞ

eðf1 Þ

Algorithm 12 20 20

4.7 11 11 19

1.303e  09 2.057e  11 6.505e  12

2.624e  06 5.527e  08 4.370e  13

3.257e  14 1.703e  13 1.202e  14

Algorithm 12 20 20

4.8 11 11 19

6.971e  10 2.199e  11 6.926e  12

1.390e  06 5.387e  08 9.948e  14

5.383e  14 9.778e  14 1.457e  14

ARTICLE IN PRESS D. Potts, M. Tasche / Signal Processing 90 (2010) 1631–1642

Table 2 Results for f1 from Example 6.2.

Table 4 Results for f3 from Example 6.4.

eðxÞ

eðqÞ

eðf1 Þ

SNR

Algorithm 4.7 20 11 60 11 60 19 60 29

9.876e  01 2.557e  02 1.098e  03 4.135e  04

3.992e  02 9.801e  04 1.017e  04 3.744e  06

1.338e  04 2.807e  04 1.771e  05 5.130e  06

52.64 52.61 52.53 52.81

Algorithm 4.8 20 11 60 11 60 19 60 29

4.006e 01 5.197e  02 2.536e  04 5.493e  04

5.283e  04 4.762e  04 3.421e  05 2.362e  06

3.070e  05 3.589e  04 4.991e  06 6.440e  06

52.85 52.60 52.94 52.54

N

L

Table 3 Results for f2 from Example 6.3. N

L

Algorithm 4.7 9 8 50 8 50 49 Algorithm 4.8 9 8 50 8 50 49

eðxÞ

eðqÞ

eðf2 Þ

4.996e  15 7.883e  15 9.548e  15

1.037e  12 2.728e  12 1.536e  11

7.123e  15 3.045e  14 1.664e  13

5.662e  15 8.105e  15 2.442e  15

1.847e  13 4.066e  13 1.423e  11

f3 ðxÞ :¼ 34 þ 300cosðpx=4Þ þ cosðpx=2Þ and sample f3 at the nodes x ¼ k ðk ¼ 0; . . . ; 2NÞ. Note that in [22], uniform noise in the range ½0; 3 is added and this experiment is repeated 500 times. Then the frequencies are recovered very accurately. In this case, the noise exceeds the lowest coefficient of f3 . Now we use Algorithms 4.7 and 4.8, where we add noise in the range ½0; R to the sampled data f3 ðkÞ ðk ¼ 0; . . . ; 2NÞ. For R ¼ 1 and 3 we obtain that the SNR is approximately 28 and 24, respectively. Of course, we can compute the three frequencies by choosing L ¼ 3 and select all frequencies, i.e., we choose e1 ¼ 0. However, in this case the error eðqÞ has a high oscillation depending on the noise, therefore we repeat the experiment also 50 times and present the averages of the errors in Table 4. Example 6.5. Finally, we consider the function 40 X

L

R

eðxÞ

eðqÞ

eðf3 Þ

Algorithm 32 64 128 256 512 1024

4.7 5 5 5 5 5 5

1 1 1 1 1 1

4.718e  03 3.438e  03 1.672e  03 1.205e  03 1.192e  03 6.521e  04

6.263e  01 6.028e  01 6.686e  01 6.134e  01 6.512e  01 6.505e  01

3.468e  03 3.448e  03 3.595e  03 3.901e  03 4.519e  03 4.601e  03

32 64 128 256 512 1024

5 5 5 5 5 5

3 3 3 3 3 3

4.336e  02 2.257e  02 1.711e  02 1.089e  02 7.551e  03 6.261e  03

1:716e þ 00 1:969e þ 00 1:405e þ 00 1:503e þ 00 1:487e þ 00 1:710e þ 00

9.667e  03 1.037e  02 1.058e  02 1.154e  02 1.011e  02 1.118e  02

Algorithm 32 64 128 256 512 1024

4.8 5 5 5 5 5 5

1 1 1 1 1 1

3.992e  03 2.953e  03 2.068e  03 1.329e  03 8.220e  04 7.259e  04

5.531e  01 6.538e  01 5.835e  01 6.833e  01 7.226e  01 5.172e  01

3.348e  03 3.501e  03 3.537e  03 3.818e  03 4.140e  03 4.543e  03

32 64 128 256 512 1024

5 5 5 5 5 5

3 3 3 3 3 3

3.849e  02 2.254e  02 1.745e  02 1.152e  02 9.114e  03 7.110e  03

1:618e þ 00 1:835e þ 00 1:680e þ 00 2:025e þ 00 2:159e þ 00 1:866e þ 00

1.047e  02 1.054e  02 1.073e  02 1.091e  02 1.158e  02 1.102e  02

N

6.1500e  15 8.6300e  14 1.0765e  13

Example 6.4. A different method for finding the frequencies based on detection of singularities is suggested in [22], where in addition this method is also suitable for more difficult problems, such as finding singularities of different orders and estimating the local Lipschitz exponents. Similarly as in [22], we consider the 8-periodic test function

f4 ðxÞ :¼

1641

cosðoj xÞ;

j¼1

where we choose random frequencies oj drawn from the standard uniform distribution on ð0; pÞ. We sample the function f4 at the equidistant nodes x ¼ k ðk ¼ 0; . . . ; 2NÞ. For the results of Algorithms 4.7 and 4.8 see Table 5. Note

Table 5 Results for f4 from Example 6.5. eðxÞ

eðqÞ

eðf4 Þ

Algorithm 4.7 170 84 170 169 512 84 512 511 1024 84 1024 500 1024 1000

1.4325e  11 2.5629e  10 4.6367e  08 1.1102e  14 1.0560e  10 1.9984e  14 1.2254e  14

5.8613e  10 3.2084e  07 8.3982e  06 7.8241e  12 1.7671e  07 2.9797e  11 1.0122e  11

6.2146e  10 9.7832e 12 4.8472e  05 3.0388e  11 1.1728e  07 7.8334e  11 8.5558e 11

Algorithm 4.8 170 84 170 169 512 84 512 509 1024 1000

7.5935e  08 3.1123e  11 2.7270e  05 2.7420e  12 7.1054e  15

6.3152e  06 2.7379e  09 1.3895e  02 1.0926e  09 2.5256e  09

2.4869e  08 7.0488e  12 3.1438e  02 3.8504e  11 7.9915e  11

N

L

that without the diagonal preconditioner D the overdetermined linear Vandermonde-type system is not numerically solvable due to the ill-conditioning of the Vandermonde-type matrix. In summary, we obtain very accurate results already for relatively few sampled data. We can analyze both periodic and nonperiodic functions without preprocessing (as filtering or windowing). APM works correctly for noisy sampled data and for clustered frequencies assumed that the separation distance of the frequencies is not too small and the number 2N þ 1 of sampled data is sufficiently large. We can also use the ESPRIT method in order to determine the frequencies. Overmodeling improves the

ARTICLE IN PRESS 1642

D. Potts, M. Tasche / Signal Processing 90 (2010) 1631–1642

results. The numerical examples show that one can choose the number 2N þ 1 of sampled data less than expected by the theoretical results. We have essentially improved the stability of our APM by using a weighted least squares method. Our numerical examples confirm that the proposed APM is robust with respect to noise. Acknowledgements The first named author gratefully acknowledges the support by the German Research Foundation within the Project KU 2557/1-1. Moreover, the authors would like to thank the referees for their valuable suggestions. References [1] G. Beylkin, L. Monzo´n, On approximations of functions by exponential sums, Appl. Comput. Harmon. Anal. 19 (2005) 17–48. [2] W.L. Briggs, V.E. Henson, The DFT, SIAM, Philadelphia, 1995. [3] J.M. Papy, L. De Lathauwer, S. Van Huffel, Exponential data fitting using multilinear algebra: the single-channel and multi-channel case, Numer. Linear Algebra Appl. 12 (2005) 809–826. [4] D. Potts, M. Tasche, Numerical stability of nonequispaced fast Fourier transforms, J. Comput. Appl. Math. 222 (2008) 655–674. [5] S.L. Marple Jr., Digital spectral analysis with applications, in: Prentice Hall Signal Processing Series, Prentice Hall Inc., Englewood Cliffs, 1987. [6] G.H. Golub, P. Milanfar, J. Varah, A stable numerical method for inverting shape from moments, SIAM J. Sci. Comput. 21 (1999) 1222–1243. [7] P.L. Dragotti, M. Vetterli, T. Blu, Sampling moments and reconstructing signals of finite rate of innovation: Shannon meets StrangFix, IEEE Trans. Signal Process. 55 (2007) 1741–1757. [8] R.A. Horn, C.R. Johnson, Matrix Analysis, Cambridge University Press, Cambridge, 1985.

[9] D. Potts, G. Steidl, M. Tasche, Fast Fourier transforms for nonequispaced data: a tutorial, in: J.J. Benedetto, P.J.S.G. Ferreira (Eds.), Modern Sampling Theory: Mathematics and Applications, ¨ Birkhauser, Boston, 2001, pp. 247–270. [10] J. Keiner, S. Kunis, D. Potts, NFFT 3.0, C subroutine library /http:// www.tu-chemnitz.de/potts/nfftS, 2006. [11] S. Kunis, D. Potts, Stability results for scattered data interpolation by trigonometric polynomials, SIAM J. Sci. Comput. 29 (2007) 1403–1419. ˚ Bjorck, ¨ [12] A. Numerical Methods for Least Squares Problems, SIAM, Philadelphia, 1996. [13] Y. Hua, T.K. Sarkar, Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise, IEEE Trans. Acoust. Speech Signal Process. 38 (1990) 814–824. [14] T.K. Sarkar, O. Pereira, Using the matrix pencil method to estimate the parameters of a sum of complex exponentials, IEEE Antennas Propagation Mag. 37 (1995) 48–55. [15] R. Roy, T. Kailath, ESPRIT-estimation of signal parameters via rotational invariance techniques, IEEE Trans. Acoust. Speech Signal Process. 37 (1989) 984–994. [16] R. Roy, T. Kailath, ESPRIT-estimation of signal parameters via ¨ rotational invariance techniques, in: L. Auslander, F.A. Grunbaum, J.W. Helton, T. Kailath, P. Khargoneka, S. Mitter (Eds.), Signal Processing, Part II: Control Theory and Applications, Springer, New York, 1990, pp. 369–411. [17] D.G. Manolakis, V.K. Ingle, S.M. Kogon, Statistical and Adaptive Signal Processing, McGraw-Hill, Boston, 2005. [18] P. Strobach, Fast recursive low-rank linear prediction frequency estimation algorithms, IEEE Trans. Signal Process. 44 (1996) 834–847. [19] P. Strobach, Total least squares phased averaging and 3-D ESPRIT for joint azimuth-elevation-carrier estimation, IEEE Trans. Signal Process. 59 (2001) 54–62. [20] N.J. Higham, Accuracy and Stability of Numerical Algorithms, SIAM, Philadelphia, 1996. ˚ [21] W.B. Jones, O. Njastad, H. Waadeland, Application of Szego+ polynomials to frequency analysis, SIAM J. Math. Anal. 25 (1994) 491–512. [22] H.N. Mhaskar, J. Prestin, On local smoothness classes of periodic functions, J. Fourier Anal. Appl. 11 (2005) 353–373.