Matched-Filter Bank Interpretation of Some Spectral ... - CiteSeerX
Matched-Filter Bank Interpretation of Some Spectral Estimators Petre Stoicay
Andreas Jakobssony
Jian Liz
September 26, 1997 Abstract We make use of a matched- lter bank (MAFI) approach to derive spectral estimators for stationary signals with mixed spectra. We show that the Capon spectral estimator as well as the more recently introduced APES (amplitude and phase estimation) method are members of the MAFI class. A seemingly new MAFI-based spectral estimator is shown to reduce to APES as well. The MAFI interpretation of Capon and APES, along with some additional analysis, provides further insights into the properties of and the relationship between these two spectral estimation methods.
This work was supported in part by the Senior Individual Grant Programme of the Swedish Foundation for Strategic Research, the Swedish National Board for Technical Development (NUTEK), the NSF Grant MIP-9457388 and the ARPA Grant MDA-972-93-1-0015. y Petre Stoica and Andreas Jakobsson are with the Systems and Control Group, Department of Technology, Uppsala University, P.O. Box 27, S-751 03 Uppsala, Sweden. z Jian Li is with the Department of Electrical Engineering, University of Florida, Gainesville, FL 32611, USA.
1
1 Introduction The derivation of the Capon spectral estimator (see, e.g., [1, 2, 3]) is quite related to that of a matched lter bank processor. However the fact that Capon is indeed a matched- lter bank spectral estimator does not appear to be widely known: to the best of our knowledge only the recent reference [4] addressed this connection. In the present paper we rst introduce a matched- lter bank (MAFI) approach to spectral estimation for stationary signals with mixed (i.e. both continuous and discrete) spectra. Then we show that Capon is a member of the MAFI class of spectral estimators. We also show that the more recently introduced APES (amplitude and phase estimation) method ([5]) belongs to the MAFI class as well. The belonging of APES to MAFI was conjectured in [5]. The MAFI approach to spectral estimation may be used to devise new spectral estimation methods. Interestingly enough, a most reasonable implementation of a seemingly novel MAFI-based spectral estimator, derived in this paper, is also shown to reduce to the APES method. The MAFI interpretation of Capon and APES, along with some additional analysis, provides further insights into the properties of and the relation between these two spectral estimators. In particular, we prove that the APES spectral estimate is lower bounded by the Capon estimate at all frequencies. Additionally, we show by means of a higher-order expansion technique that the Capon estimator underestimates the true spectrum whereas the APES spectral estimator is unbiased (to within a second-order approximation). These theoretical results, supplemented with practical evidence that the Capon estimator is indeed biased downwards in samples of practical length whereas APES is nearly unbiased (see [5] and also Section 6 of this paper), are believed to provide a compelling reason for preferring APES to Capon.
2 The lter bank approach Let fx(t); t = 1; 2; : : : ; N g denote the available sample of the stationary signal the spectrum of which is to be estimated, where N denotes the number of data samples. The lter bank approach basically reduces the problem of estimating the spectrum of x(t) to that of estimating the amplitude (or power) of a sinusoidal signal buried in colored noise (see, e.g., [5, 3]). More exactly, in the aforementioned approach x(t) is additively decomposed as: x(t) = (!)ei!t + (t); t = 1; 2; : : : ; N; ! 2 (0; 2] (2.1) where (!) denotes the (complex-valued) amplitude of the sinusoidal signal referred to above and (t) is the noise (or residual) term, assumed to be zero-mean. The three main steps of any lter bank spectral estimator can be succinctly described as follows: Step (a) Pass the data fx(t)g through a bandpass lter with (varying) center frequency !, to enhance the sinusoidal component in (2.1). Step (b) Estimate (!) from the ltered data. Let ^ (!) denote the estimated amplitude. Step (c) Take ^(!), for ! 2 (0; 2], as an estimate of the complex spectrum of x(t); or j ^(!)j2 (appropriately normalized) as an estimate of the power spectrum. We assume that the lter used in Step (a) has a nite impulse response (FIR). Let h
h! = h1 : : : hL
i
(2.2) denote the vector of FIR coecients, where L denotes the lter's length and the superscript stands for the conjugate transpose operator. Observe that the dependence of the vector in (2.2) on the center
2
frequency ! is stressed by notation. By using (2.2) we can write the lter output as: z (t) = h! y(t) (2.3) where h iT y(t) = x(t) : : : x(t + L ? 1) (t = 1; : : : ; M ) (2.4) Hereafter, the superscript T denotes the transpose and M =N ?L+1 (2.5) The choice of L above should be done by a compromise between resolution and statistical stability: the larger L the better the resolution but the worse the statistical stability. It is not dicult to understand why this is so. The time-bandwidth product type of result implies that the lter can be made narrower in the frequency domain, and hence the corresponding spectral resolution enhanced, as L increases. On the other hand, M decreases with increasing L and therefore we have fewer samples for the estimation of (!) in Step (b), which should lead to a poorer statistical accuracy. Next we discuss the choice of h! . In the Capon method h! is chosen such that the power of the lter output is minimized subject to the constraint that the frequency ! is passed undistorted: min h Rh subject to h! a(!) = 1 (2.6) h ! ! where
!
h
a(!) = 1 ei! : : : ei(L?1)!
and
iT
(2.7)
R = E [y(t)y (t)]
(2.8) Hereafter E denotes the statistical expectation operator. The solution to the lter design problem in (2.6) is given by [1, 2] : R?1 a(!) (2.9) h! = a (!)R?1 a(!) assuming that R is invertible. Of course, (2.9) is not directly implementable since R is not available. The most commonly used implementable form is obtained by using the sample covariance matrix, M 1 X y(t)y (t) (2.10) R^ = M t=1
in lieu of R in (2.9): R^ ?1a(!) hCapon = (2.11) ! a (!)R^ ?1 a(!) The Fourier ([6, 7, 8]) and APES ([5]) spectral estimators have not been derived in the lter bank approach framework. However they can be cast into that framework by simple manipulations. The corresponding FIR lter vectors turn out to be (see, [6, 7] and [5], respectively): a(!) hFourier = (2.12) !
and where
L
^ ?1 hAPES = Q ^(?!1)a(!) ! a (!)Q (!)a(!)
(2.13)
Q^ (!) = R^ ? Y (!)Y (!)
(2.14)
3
and where Y (!) =
M 1 X y(t)e?i!t M
(2.15)
t=1
In the next section we propose the use of the matched lter (MAFI) approach to derive the lter in Step (a) of the lter bank-based spectral estimation methodology outlined above. Even though the matched lter should be a natural choice for Step (a), apparently its use in the previous context was not directly suggested before.
3 The matched- lter bank approach By making use of the notation h
n(t) = (t) : : : (t + L ? 1)
and
iT
(3.1)
s(t) = (!)ei!t a(!)
(3.2)
we can write y(t) in (2.4) as: y(t) = s(t) + n(t) (3.3) We assume that the initial phase of the sinusoidal signal in the equation above is a random variable which is uniformly distributed over the interval (0; 2] and is independent of the noise term. Then it follows that the covariance matrix of the signal term in (3.3) is given by: E [s(t)s (t)] = j (!)j2 a(!)a (!) (3.4) Let Q(!) denote the covariance matrix (assumed to be invertible) of the noise term in Equation (3.3): Q(!) = E [n(t)n (t)] (3.5) The notation Q(!) is meant to indicate the fact that the noise (or residual) term in (3.3), and hence its covariance, depends on !. By de nition the matched lter should be designed such that the corresponding signal-to-noise ratio (SNR) in the lter's output is maximized [9]: jh! a(!)j2 (3.6) max h! h! Q(!)h! The solution to the design problem above is well known [6, 7, 8] (for completeness we provide a simplest derivation in Appendix A: Q?1 (!)a(!) (3.7) h! = a (!)Q?1 (!)a(!) The covariance matrix in (3.7) can be estimated in several ways. As shown below two natural estimators of Q(!) lead to the Capon and APES lters, respectively. A third estimate of Q(!) yields a new lter that is di erent from both the Capon and APES lters. However, as shown in the next section, the spectral estimator corresponding to the new lter turns out to be equivalent to APES as well.
3.1 Capon Since
R = j (!)j2 a(!)a (!) + Q(!)
(3.8)
we can estimate Q(!) as: Q^ (!) = R^ ? j ^(!)j2 a(!)a (!)
(3.9)
4
where ^(!) is an estimate of (!). As we show below the second term in (3.9) has no in uence on the h! in (3.7). This means that for the Q^ (!) in (3.9) there will be in fact no need to estimate (!). A straightforward use of the matrix inversion lemma (stated in Appendix B for easy reference) yields1 : 2 ^ ?1 ^ ?1 R^ ?1 a (3.10) (R^ ? j ^j2 aa )?1 a = R^ ?1a + j ^j R 2aa ^R?1 a = 1 ? j ^j a R a 1 ? j ^j2 a R^ ?1 a Since the scaling of Q^ ?1 a in the expression of the sample counterpart of h(!) in (3.7) leaves the FIR vector unchanged, it follows from (3.10) that the estimated matched lter corresponding to the Q^ in (3.9) coincides with the Capon lter. A perhaps even more direct way to see that the Capon lter is a matched lter is as follows. Consider the MAFI design problem in (3.6). As the multiplication of h! by any complex number leaves the ratio in (3.6) unchanged, it follows that the optimization problem in (3.6) is equivalent to: min h Qh subject to h! a(!) = 1 (3.11) h ! ! !
However, by using (3.8) along with the constraint in (3.11), we can write: h i h! Qh! = h! R ? j (!)j2 a(!)a (!) h! = h! Rh! ? j (!)j2 The conclusion is that the Capon design problem is nothing but a matched- lter design.
(3.12)
3.2 APES By combining (3.2) and (3.3) we obtain: y(t) = [ (!)a(!)] ei!t + n(t) (3.13) The least-squares (LS) estimate of the vector (!)a(!) in (3.13), which ignores the fact that a(!) is known, is given by the normalized Fourier transform (see, e.g., [10]): M X 4 Y (!) y(t)e?i!t = [ (!d )a(!)] = M1
(3.14)
t=1
Inserting (3.14) into (3.9) yields the following estimate of Q(!): Q^ (!) = R^ ? Y (!)Y (!) the use of which leads to the APES lter in (2.13), (2.14).
(3.15)
3.3 MAFI Equation (3.13) suggests another way to estimate the covariance matrix Q(!): Q^ (!) =
M 1 X y(t) ? ^(!)a(!)ei!t y(t) ? ^(!)a(!)ei!t M t=1
(3.16)
where, as before, ^(!) denotes an estimate of (!). A simple calculation shows that the previous Q^ (!) matrix can be rewritten as follows: Q^ (!) = R^ ? ^ (!)Y (!)a (!) ? ^(!)a(!)Y (!) + j ^(!)j2 a(!)a (!) (3.17) By a calculation similar to (3.10) one can verify that the last term in (3.17) can be dropped when evaluating the MAFI vector. Below we make use of the matrix inversion lemma once more to show that 1
For notational convenience we omit the dependence on !.
5
also the third term in (3.17) can be omitted. Let be the matrix made from the rst two terms in (3.17). Then, by the matrix inversion lemma, ?1 ?1 ?1 (3.18) ( ? ^aY )?1 a = ?1 + 1^ ? ^YaY
?1 a a = 1 ? ^ Y a ?1 a where the dependence on ! has been omitted to simplify the notation. A further use of the matrix inversion lemma yields
?1 a = (R^ ? ^ Y a )?1 a = # " ^ R^ ?1 Y a R^ ?1 ? 1 ^ a= = R + 1 ? ^ a R^ ?1 Y ^ ?1 ^ ?1 ^ ?1 ^ ?1 ^ ?1 = R a ? ^ (a R Y )R a^ ?+1 ^ (a R a)R Y ; (3.19) 1 ? ^ a R Y which gives the following expression for the MAFI vector: h i ^ ?1a(!) + ^ (!) (a (!)R^ ?1 a(!))R^?1 Y (!) ? (a (!)R^ ?1 Y (!))R^ ?1 a(!) R (3.20) hMAFI = ! a (!)R^ ?1 a(!) The previous lter is generally di erent from both the Capon and the APES lters. In e ect, neither of the latter two depends on an estimate of (!), unlike the former. (The estimation of (!), needed by (3.20), is discussed in the next section). Despite this fact, in the next section we show the somewhat surprising result that a most reasonable implementation of the spectral estimator corresponding to (3.20) is identical to APES!
4 Spectrum estimation It is readily checked that all lter vectors previously discussed satisfy: h! a(!) = 1 By making use of this observation and of (3.13) we obtain: h! y(t) = (!)ei!t + h! n(t) The LS estimator of (!) in the above equation is given by (see, e.g., [10]) M 1 X ^(!) = h! y(t)e?i!t = h! Y (!) M t=1
(4.1) (4.2) (4.3)
The Capon and APES estimates of the complex spectrum are obtained by inserting the corresponding lter vector expressions in (4.3). The use of the MAFI lter derived in the previous section for spectrum estimation requires an initial estimate of (!). That estimate might be obtained by any of the other methods discussed herein (i.e., Fourier, Capon or APES), but proceeding in that way would be computationally unacceptable. Actually, we can avoid the need for an initial estimate of (!) in the following way. Inserting (3.20) into (4.3) yields
2
^(a R^ ?1 a) = a R^ ?1 Y + ^ (a R^ ?1 a)(Y R^ ?1 Y ) ? a R^ ?1 Y
from which we obtain (by assuming that the ^'s in the right-hand and left-hand sides are identical): a (!)R^ ?1 Y (!) (4.4) ^MAFI (!) = 2 ? 1 ? 1 ? 1 ? 1 ^ ^ ^ ^ a (!)R a(!) ? (a (!)R a(!))(Y (!)R Y (!)) ? a (!)R Y (!)
6
Next we show that the previous spectral estimator coincides with APES: ^APES (!) = ^MAFI (!) (4.5) The APES estimate of the complex spectrum is given by (once again we omit the dependence on ! for notational convenience): a (R^ ? Y Y )?1 Y ^APES (!) = ^ = a (R ? Y Y )?1 a ?1 ?1 a R^ ?1 + R^1?YYYR^?R^1 Y Y = = a R^ ?1 + R^1??1YYYR^?R^1?Y1 a a R^ ?1 Y = (4.6) ^ ?1 2 ? 1 ? 1 ? 1 ^ ^ ^ a R a ? (a R a)(Y R Y ) ? a R Y which is identical to (4.4). Hence the equality in (4.5) is proved. For comparative purposes, note that the Capon spectral estimator is given by: a (!)R^?1 Y (!) ^Capon (!) = ^ ?1 a (!)R a(!) The next section discusses the relative merits of Capon and APES.
(4.7)
5 Comparisons and discussions 5.1 The computational issue
Let R^ ?1=2 denote a square root of the positive de nite matrix R^ ?1 , and let: (!) = R^ ?1=2 a(!) (!) = R^ ?1=2 Y (!)
The Capon and APES spectral estimators can be expressed as relatively simple functions of the (!) and (!) introduced above: (!)(!) 4 (!)(!) ^Capon (!) = (5.1) = (!) (!) k (!)k2 (!)(!) (5.2) ^APES (!) = k (!)k2 ? k (!)k2 k(!)k2 ? j (!)(!)j2 The ecient estimation of the Capon spectrum estimate is discussed, for example, in [6, 7, 8]. We do not intend to dwell into the details of that computation. The point we want to make here is that APES is only slightly more involved computationally than Capon (as should be evident by comparing (5.1) and (5.2)), and hence the ecient computational means developed for the latter can also be applied to the former.
5.2 The statistical performance issue The Capon, APES as well as the Fourier spectral estimators can be shown to have the same asymptotic variance under the following condition:
7
C: The signal x(t) can be written as in (2.1), where (t) is a zero-mean stationary random process with nite spectral density at !: " (!) < 1
(5.3)
In more exact terms, the following result holds true. Theorem 1 Under condition C and the additional assumption that (t) is circularly symmetrically distributed, the estimation errors in the Fourier, Capon and APES spectral estimators are asymptotically circularly symmetrically distributed with zero-mean and the following common variance: lim ME j ^(!) ? (!)j2 = (!) (5.4) M !1 Proof: See Appendix C. The need to enforce condition C limits, to some extent, the importance of the previous result. Indeed the assumption made in C is satis ed if (and essentially only if) the signal x(t) has a mixed spectrum and ! is the location of a spectral line. In some applications, such as target feature extraction by means of a Doppler radar, the signal has a mixed spectrum and estimating the amplitudes of the sinusoidal components is a key problem (see, e.g., [5] and the references therein). The result of Theorem 1 is relevant to such applications. In other applications, however, the main interest is in the continuous component of the spectrum. Condition C does not hold exactly if the spectrum is continuous at !. The previous result is of a somewhat limited interest also because of its asymptotic character. Indeed, in simulations with medium or small-sized samples, the spectral estimators under study have been found to behave quite di erently from one another in contradiction with what is predicted by the (asymptotic) result of Theorem 1 (see, e.g., [5]; also see the next section). In particular, in scenarios with multiple, closely spaced spectral lines the Fourier estimate of the complex spectrum has been found to be (sometimes signi cantly) less accurate than the spectrum estimates obtained by using the Capon and APES methods [5, 3]. The nite-sample analysis of the spectral estimators under discussion would be of considerable interest. However, while this is possible for the Fourier estimator, a complete analysis of the more interesting, data-adaptive Capon and APES estimators appears to be quite dicult at best. A partial analysis of the latter estimators is nevertheless possible as shown in what follows. Speci cally, in Appendix D we make use of a higher-order Taylor series expansion technique to prove that, to within a second-order approximation, and under the additional mild assumption that the third-order moments of n(t) are zero, Capon is biased downwards whereas APES is unbiased: ME [^ Capon (!) ? (!)] < 0 (! 2 [?; )) (5.5) (!) and ME [^ APES (!) ? (!)] = 0 (! 2 [?; )) (5.6) for suciently large values of M . The above result provides theoretical support to the empirically observed fact that Capon underestimates the complex spectrum, whereas the APES estimate is almost unbiased (see [5]). In fact the Capon and APES spectral estimators satisfy the following neat inequality, which points to the type of result proved above: j ^Capon (!)j j ^APES (!)j and arg(^ Capon (!)) = arg(^ APES (!)) , 8! 2 [?; ) (5.7)
8
where the modulus equality holds only for M ! 1. To check (5.7) observe that the term between square parentheses in (4.6) is positive (by the Cauchy-Schwartz inequality) and that it can be zero if and only if Y (!) is proportional to a(!) (which can hold true only when M ! 1). We believe that (5.7) along with (5.5) and (5.6) provide a theoretical motivation for preferring APES to Capon in most spectrum estimation exercises. The simulations in the next section lend further evidence to the fact that APES outperforms Capon in the nite-sample case.
6 Numerical examples In this section we study the accuracy of the APES and Capon amplitude estimates in a number of cases.
6.1 Estimation performance versus SNR We begin by studying the performance of using the APES and Capon methods for complex amplitude estimation as the SNR varies. Figure 1 shows the modulus of the true spectrum of the sinusoidal signal, which consists of four dominant spectral lines and nine small spectral lines located at the following frequencies: 0.0625, 0.0875, 0.25, 0.285, 0.33, 0.35, 0.37, 0.39, 0.41, 0.43, 0.45, 0.47, 0.49. These spectral lines all have a phase value of =4. The data sequence has N = 64 data samples and is corrupted by complex white Gaussian noise with zero mean and variance 2 . The SNR for the kth spectral line is de ned as
(6.1) SNRk = 10 log10 j k2j [dB] where k is the complex amplitude of the kth spectral line. In the following, we compare the biases and variances of the APES (solid lines in the gures) and Capon (dashed lines in the gures) estimates of the complex amplitude at the frequency of the rst dominant spectral line, as its SNR varies. Our results are obtained from 100 Monte Carlo trials in which we used 100 independent sequences of circular Gaussian white noise to generate the simulation data set. 2
Figures 2(a) and 2(b) compare the real and imaginary parts, respectively, of the biases of the Capon and APES estimators when L = 24. As seen from the gures, the APES estimates are nearly unbiased, while the Capon estimates are biased downward, as was expected from our previous theoretical analysis. Figure 2(c) and 2(d) compare the variances and the MSE, respectively, of the Capon and APES estimators when L = 24. Note that the APES estimator provides not only a smaller bias but also a lower variance than the Capon estimator.
6.2 Estimation performance versus lter length Next we study the e ect of the lter length, L, on the performance of the two estimators. The data studied is the same as in the previous example, with the SNR of the rst dominant spectral line xed at 20 dB. Our results are obtained from 100 Monte Carlo trials. Figures 3(a) and 3(b) compare the real and imaginary parts, respectively, of the biases of the Capon and APES estimators as L varies. As seen from the gures, the APES estimates are again almost unbiased for not too small lter lengths, while the Capon estimates have very large biases for large L. Figure 3(c) and 3(d) compares the variances and the MSE, respectively, of the Capon and APES estimators. Note that the variance of the APES estimate varies little for the studied lter lengths, while the variance of the Capon estimate signi cantly increases with increasing L. This type of behavior makes the
9
Modulus of True Spectrum
1
0.8
0.6
0.4
0.2
0 0
0.2
0.4
0.6
0.8
1
Frequency
Figure 1: The modulus of the true spectrum. choice of L with the Capon method quite dicult. We remind the reader that to enhance the resolution L should be chosen as large as possible (subject to L < N=2 so that R^ is not singular). While APES with this choice of L works ne, Capon does not. We further illustrate this point in the next gure. In Figures 4(a) and 4(b), one realization of the modulus of the estimated spectra obtained by using APES and Capon, respectively, is compared with the true spectrum (dotted lines) for L = 15. Figures 4(c) and 4(d) are similar to Figures 4(a) and 4(b) except that the lter length is now increased to L = 24. As can be seen from the gures, APES gives much more accurate peak amplitude estimates, at the price of slightly wider spectral peaks. Note that, as L increases, the resolution of both methods improves, as expected; however, unlike APES, the performance of the Capon amplitude estimates detoriates significantly with increasing L. Note also that the APES estimates, at all evaluated frequencies, are larger than the Capon estimates, as expected. In Table 1 we show the phase estimates at the frequencies of the four dominant spectral lines. In accordance with what was expected from (5.7), the APES and Capon phase estimates coincide.
10
0.02
0.02
0
0
−0.02
APES Capon Bias
Bias
−0.02
−0.04
APES Capon −0.04
−0.06
−0.06
−0.08
−0.08
−0.1
2
4
6
8
10
12
14
16
18
−0.1
20
2
4
6
8
10
SNR
12
14
16
18
20
SNR
(a)
(b)
−1
−1
10
10 APES Capon
APES Capon
−2
−2
10 MSE
Variance
10
−3
−3
10
10
−4
10
−4
2
4
6
8
10
12
14
16
18
10
20
SNR
2
4
6
8
10
12
14
16
18
20
SNR
(c)
(d)
Figure 2: Bias and variance of the complex amplitude estimate at the frequency of the rst spectral line as a function of SNR when L = 24. (a) Real part of the bias. (b) Imaginary part of bias. (c) Variance. (d) MSE.
11
0.15
0.15 APES Capon
0.05
0.05
0
0
−0.05
−0.05
−0.1
−0.1
−0.15
−0.15
−0.2
−0.2
−0.25
−0.25
−0.3
−0.3 5
10
15
20
25
APES Capon
0.1
Bias
Bias
0.1
30
5
10
15
L
20
25
(a)
(b)
0
0
10
10 APES Capon
APES Capon
−1
−1
10
MSE
Variance
10
−2
10
−3
−2
10
−3
10
10
−4
10
30
L
−4
5
10
15
20
25
10
30
L
5
10
15
20
25
30
L
(c)
(d)
Figure 3: Bias and variance of the complex amplitude estimate at the frequency of the rst spectral line as a function of L when the SNR of the rst spectral line is 20 dB. (a) Real part of the bias. (b) Imaginary part of bias. (c) Variance. (d) MSE.
12
1 Modulus of Estimated Spectrum
Modulus of Estimated Spectrum
1
0.8
0.6
0.4
0.2
0 0
0.2
0.4
0.6
0.8
0.4
0.2
0.2
0.4
Frequency
Frequency
(a)
(b)
0.6
0.8
1
0.6
0.8
1
1 Modulus of Estimated Spectrum
Modulus of Estimated Spectrum
0.6
0 0
1
1
0.8
0.6
0.4
0.2
0 0
0.8
0.2
0.4
0.6
0.8
1
0.8
0.6
0.4
0.2
0 0
0.2
0.4
Frequency
Frequency
(c)
(d)
Figure 4: Modulus of the estimated spectrum (solid lines) compared with the true spectrum (dotted lines) when the SNR of the rst spectral line is 20 dB. (a) APES with L = 15. (b) Capon with L = 15. (c) APES with L = 24. (d) Capon with L = 24.
13
Phase of spectral line no. 1 2 3 4
True Estimated values values APES, L=15 Capon, L=15 APES, L=24 Capon, L=24 0:25 0:2457 0:2457 0:2477 0:2477 0:25 0:3506 0:3506 0:3408 0:3408 0:25 0:2513 0:2513 0:2507 0:2507 0:25 0:2493 0:2493 0:2454 0:2454
Table 1: Phase estimates at the frequencies of the four dominant spectral lines.
7 Conclusions The Capon spectral estimator has been shown to be a matched- lter bank method. In fact, the Capon's lter design was shown to be precisely a matched lter design. The APES spectral estimator has also been shown to be a member of the class of matched- lter bank methods. In e ect, APES was proved to coincide with two seemingly di erent matched- lter bank spectral estimation methods. The matched- lter bank interpretation of Capon and APES provides a (partial) theoretical explanation for the good, empirically observed performance of these two spectral estimators. For the moment it is an open issue whether other matched- lter bank spectral estimation methods with enhanced properties exist. The matched- lter bank interpretation has also provided further insight into the properties of and relationships between Capon and APES. In particular the signal model used by the two methods was shown to be adequate only at the frequencies corresponding to the spectral lines in a mixed spectrum. The model in question is not valid at the frequencies of the continuous spectral component. While APES was shown by means of simulations to yield accurate estimates of the latter component as well [5], there is a clear need for a better theoretical understanding of why and how APES works in such cases (The performance of the Capon estimate of the continuous spectral component was poor in our simulations reported in [5]). Derivation of matched- lter bank methods speci cally designed for continuous spectra is also an interesting open issue. The casting of Capon and APES into the matched- lter bank framework, along with the related analysis, has also led to a neat result concerning the relationship between these two methods: the Capon spectral estimate is always smaller than APES's, for all frequencies. This theoretical observation lends support to the empirical nding that Capon often underestimates the true spectrum whereas APES is nearly unbiased. In fact by making use of an expansion technique we were able to prove that, to within a second-order approximation, Capon is indeed biased downwards whereas APES is unbiased. As APES is only slightly more involved than Capon, from a computational standpoint, preferring the former to the latter appears to be a logical conclusion of both the theoretical and empirical analysis in this paper.
14
Appendix A. Proof of (3.7) Let Q1=2 denote a Hermitian (for notational simplicity) square root of Q. Then, by the Cauchy-Schwartz inequality, h Q1=2 Q?1=2 a(! ) 2 ! a (!)Q?1 a(!) h Qh !
!
where the equality is achieved if and only if h! = Q?1 a(!) for any scalar 2 C. As is made transparent in Section 4, it is often required that h! passes the frequency ! undistorted, that is h! a(!) = 1. This requirement can be met by choosing as = ? 1 1= a (!)Q a(!) , which yields (3.7).
Appendix B. The matrix inversion lemma Let R be an n n matrix, and let b and c be n 1 vectors. Then, under the assumption that the inverse appearing below exists, ?1 ?1
bc R (R ? bc )?1 = R?1 + R 1 ? c R?1 b
(see, e.g., [6, 10]).
Appendix C. Proof of Theorem 1 To simplify the notation we omit the dependence on ! whenever there is no possibility for confusion. From (3.13) we have that Y = a + (C.1) where M X = 1 n(t)e?i!t (C.2) M t=1
The rst and second-order moments of are easy to calculate: E () = 0 E (T ) = 0
(C.3) (C.4)
(by the circularly symmetric distribution assumption) and E ( ) =
M M X M 1 X ?i!(t?s) = 1 X (M ? j j)Rn ( )e?i! R ( t ? s ) e n M2 M2 =?M
t=1 s=1
(C.5)
where fRn ( )g is the covariance sequence of n(t). It follows from (C.5) and condition C that: 1 X lim ME ( ) = Rn ( )e?i! = n (!) = (!)aa M !1 =?1
(C.6)
where the last equality follows from standard results on the transfer of spectral densities through linear systems. Among others, the previous calculations imply that, as M ! 1, Y tends to a (in the mean square sense). Hence hCapon and hAPES have the same limit as M increases without bound. ! ! Let h denote a generic FIR vector and let h1 denote the deterministic vector that is the limit of (the
15
possibly random) h when M goes to in nity. Observe that for all methods under study the associated h and h1 vectors satisfy: h a = 1 and h1 a = 1 (C.7) By using this observation along with (C.1) we obtain: 4 h Y = + h ^ =
(C.8) Because tends to zero as M ! 1 (as shown previously), it follows from (C.8) that the estimation error in ^ can asymptotically be written as (to within a rst-order approximation): ^ ? ' h1 (C.9) From (C.9) and (C.3)-(C.6) we readily derive that ME (^ ? )2 = 0 as M ! 1 and lim ME j ^ ? j2 = h1 Mlim ME ( )h1 = (!) jh1 aj2 = (!) M !1 !1 and the proof is concluded.
Appendix D. Proof of (5.5), (5.6) Once again we omit the dependence on ! for notational convenience.
Proof of (5.5) By using (2.9) and (C.8) we obtain a R^ ?1 (D.1) ^Capon ? = ^ ?1 aR a In what follows we use the symbol to denote an \asymptotic equality" that holds to within a secondorder approximation. A straightforward manipulation of (D.1) yields: a (R^ ?1 ? R?1) a R?1 + ^ ?1 = ^Capon ? = a R^ ?1 a aR a a R^ ?1 (R ? R^ )R?1 ?1 1 ? 1 + 1 = + a R a R^ ?1 a a R^ ?1 a a R?1a a R?1 a ?1 ^ ?1 ?1 ?1 ?1 ^ ?1 ? a R (aR R??R1 a)R + aaRR?1a + a R ^ ?a1 (R ??R1 )a (a R a)(a R a) ? 1 ?1 ?1 ^ ?1 ? 1 ^ ? 1 ? a R (aR R??R1 a)R + aaRR?1a + (a R )(aaRR?1(aR)2? R)R a which, in turn, implies that " # ^ ?1 )(a R?1 a) + (a R?1 RR ^ ?1 a)(a R?1 ) ? ( a R?1 RR E [^ Capon ? ] E (D.2) (a R?1 a)2 Next we note that M X 1 ^ R = aei!t + n(t) a e?i!t + n (t) = M t=1
M X n(t)n (t) = j j2 aa + a + a + M1 t=1
16
and also remind the reader of the de nition (C.2) of and the assumption that n(t) has zero third-order moments. By making use of these facts, along with (D.2), we can write: (a R?1 a)2 E (^ Capon ? ) E ?(a R?1 a)(a R?1 a R?1 + a R?1 a R?1 )+ + (a R?1 )( a R?1 a R?1a + a R?1 a R?1 a) = = E ? (a R?1 a)(a R?1 )2 ? (a R?1 a)2 ( R?1 )+ o + (a R?1 )2 (a R?1 a) + (a R?1 a) a R?1 2 = n
o
= (a R?1 a)E a R?1 2 ? (a R?1a)( R?1 ) (D.3) The quantity between curly parentheses in (D.3) is negative and thus is its expectation (by the CauchySchwartz inequality). Hence (5.5) follows.
Proof of (5.6) The APES estimate of the complex spectrum is given by (cf. (4.5)): a R^ ?1 Y ^APES = ^ ?1 a R a? where 2 = (a R^ ?1a)(Y R^ ?1 Y ) ? a R^ ?1Y
(D.4) (D.5)
By using (C.1) in (D.5) we obtain: = (a R^ ?1a) j j2 a R^ ?1a + a R^ ?1 + R^ ?1a + R^ ?1 ? j j2 (a R^ ?1 a)2 ?
2
? (a R^ ? a)( R^ ? a) ? (a R^ ? )(a R^ ? a) ? a R^ ? = 1
1
1
1
1
2 4 = (a R^ ?1a)( R^ ?1 ) ? a R^ ?1 (a R?1a)( R?1 ) ? a R?1 2 = Next note that 1 = ^1?1 + ^ ?11 ? ^1?1 = ^1?1 + ^ ?1 ^ ?1 ? 1 ^ a R a? a R a a R a ? a R a a R a (a R a)(a R a ? ) Consequently, a R^ ?1 ( a + ) ? (a R^ ?1 a) + a R^ ?1 + ^APES ? = = ^ ?1 = a R^ ?1 a ? a R a? ^ ?1 ^ ?1 = a R ^?1+ + ^(?a1 R ^+?1 ) aR a (a R a)(a R a ? ) ? 1 ^ (D.6) a R^ ?1 + a R a a R^ ?1 a The rst term in (D.6) is equal to (^ Capon ? ) and hence it follows from the rst part of this appendix that the expectation of this term satis es (to within a second-order approximation): # " ? E () a R^ ?1 (D.7) E ^ ?1 ^ ?1 aR a aR a
Combining (D.6) and (D.7) yields: E (^ APES ? ) 0 and the proof is almost concluded. It only remains to motivate the normalizing factor M used in both (5.5) and (5.6). However the factor in question follows relatively easily from the fact that is O(1=M 1=2 ) (in the mean square sense) and this implies that the second-order approximations previously used for both E (^ Capon ? ) and E (^ APES ? ) are O(1=M ).
17
References [1] J. Capon, \High Resolution Frequency Wave Number Spectrum Analysis", Proc. IEEE, 57:1408{ 1418, 1969. [2] R. Lacoss, \Data adaptive spectral analysis methods", Geophysics, 36:134{148, 1971. [3] M.A. Lagunas, M.E. Santamaria, A. Gasull, and A. Moreno, \Maximum likelihood lters in spectral estimation problems", Signal Processing, 10(1):19{34, 1986. [4] M. Wax, \Model-based processing in sensor arrays", In S. Haykin, editor, Advances in Spectral Analysis and Array Processing, volume III, pages 1{47. Prentice-Hall, Englewood Cli s, N.J., 1995. [5] J. Li and P. Stoica, \Adaptive Filtering Approach to Spectral Estimation and SAR Imaging", IEEE Trans. on Sig. Proc., 44(6):1469{1484, June 1996. [6] S.M. Kay, Modern Spectral Estimation: Theory and Application, Prentice-Hall, Englewood Cli s, N.J., 1988. [7] S.L. Marple, Jr., Digital Spectral Analysis, Prentice-Hall, Englewood Cli s, N.J., 1987. [8] P. Stoica and R. Moses, Introduction to Spectral Analysis, Prentice Hall, Upper Saddle River, NJ, 1997. [9] B. Picinbono, Random signals and systems, Prentice-Hall, Englewood Cli s, N.J., 1993. [10] T. Soderstrom and P. Stoica, System Identi cation, Prentice Hall International, London, UK, 1989.
18
Andreas Jakobssony
Jian Liz
September 26, 1997 Abstract We make use of a matched- lter bank (MAFI) approach to derive spectral estimators for stationary signals with mixed spectra. We show that the Capon spectral estimator as well as the more recently introduced APES (amplitude and phase estimation) method are members of the MAFI class. A seemingly new MAFI-based spectral estimator is shown to reduce to APES as well. The MAFI interpretation of Capon and APES, along with some additional analysis, provides further insights into the properties of and the relationship between these two spectral estimation methods.
This work was supported in part by the Senior Individual Grant Programme of the Swedish Foundation for Strategic Research, the Swedish National Board for Technical Development (NUTEK), the NSF Grant MIP-9457388 and the ARPA Grant MDA-972-93-1-0015. y Petre Stoica and Andreas Jakobsson are with the Systems and Control Group, Department of Technology, Uppsala University, P.O. Box 27, S-751 03 Uppsala, Sweden. z Jian Li is with the Department of Electrical Engineering, University of Florida, Gainesville, FL 32611, USA.
1
1 Introduction The derivation of the Capon spectral estimator (see, e.g., [1, 2, 3]) is quite related to that of a matched lter bank processor. However the fact that Capon is indeed a matched- lter bank spectral estimator does not appear to be widely known: to the best of our knowledge only the recent reference [4] addressed this connection. In the present paper we rst introduce a matched- lter bank (MAFI) approach to spectral estimation for stationary signals with mixed (i.e. both continuous and discrete) spectra. Then we show that Capon is a member of the MAFI class of spectral estimators. We also show that the more recently introduced APES (amplitude and phase estimation) method ([5]) belongs to the MAFI class as well. The belonging of APES to MAFI was conjectured in [5]. The MAFI approach to spectral estimation may be used to devise new spectral estimation methods. Interestingly enough, a most reasonable implementation of a seemingly novel MAFI-based spectral estimator, derived in this paper, is also shown to reduce to the APES method. The MAFI interpretation of Capon and APES, along with some additional analysis, provides further insights into the properties of and the relation between these two spectral estimators. In particular, we prove that the APES spectral estimate is lower bounded by the Capon estimate at all frequencies. Additionally, we show by means of a higher-order expansion technique that the Capon estimator underestimates the true spectrum whereas the APES spectral estimator is unbiased (to within a second-order approximation). These theoretical results, supplemented with practical evidence that the Capon estimator is indeed biased downwards in samples of practical length whereas APES is nearly unbiased (see [5] and also Section 6 of this paper), are believed to provide a compelling reason for preferring APES to Capon.
2 The lter bank approach Let fx(t); t = 1; 2; : : : ; N g denote the available sample of the stationary signal the spectrum of which is to be estimated, where N denotes the number of data samples. The lter bank approach basically reduces the problem of estimating the spectrum of x(t) to that of estimating the amplitude (or power) of a sinusoidal signal buried in colored noise (see, e.g., [5, 3]). More exactly, in the aforementioned approach x(t) is additively decomposed as: x(t) = (!)ei!t + (t); t = 1; 2; : : : ; N; ! 2 (0; 2] (2.1) where (!) denotes the (complex-valued) amplitude of the sinusoidal signal referred to above and (t) is the noise (or residual) term, assumed to be zero-mean. The three main steps of any lter bank spectral estimator can be succinctly described as follows: Step (a) Pass the data fx(t)g through a bandpass lter with (varying) center frequency !, to enhance the sinusoidal component in (2.1). Step (b) Estimate (!) from the ltered data. Let ^ (!) denote the estimated amplitude. Step (c) Take ^(!), for ! 2 (0; 2], as an estimate of the complex spectrum of x(t); or j ^(!)j2 (appropriately normalized) as an estimate of the power spectrum. We assume that the lter used in Step (a) has a nite impulse response (FIR). Let h
h! = h1 : : : hL
i
(2.2) denote the vector of FIR coecients, where L denotes the lter's length and the superscript stands for the conjugate transpose operator. Observe that the dependence of the vector in (2.2) on the center
2
frequency ! is stressed by notation. By using (2.2) we can write the lter output as: z (t) = h! y(t) (2.3) where h iT y(t) = x(t) : : : x(t + L ? 1) (t = 1; : : : ; M ) (2.4) Hereafter, the superscript T denotes the transpose and M =N ?L+1 (2.5) The choice of L above should be done by a compromise between resolution and statistical stability: the larger L the better the resolution but the worse the statistical stability. It is not dicult to understand why this is so. The time-bandwidth product type of result implies that the lter can be made narrower in the frequency domain, and hence the corresponding spectral resolution enhanced, as L increases. On the other hand, M decreases with increasing L and therefore we have fewer samples for the estimation of (!) in Step (b), which should lead to a poorer statistical accuracy. Next we discuss the choice of h! . In the Capon method h! is chosen such that the power of the lter output is minimized subject to the constraint that the frequency ! is passed undistorted: min h Rh subject to h! a(!) = 1 (2.6) h ! ! where
!
h
a(!) = 1 ei! : : : ei(L?1)!
and
iT
(2.7)
R = E [y(t)y (t)]
(2.8) Hereafter E denotes the statistical expectation operator. The solution to the lter design problem in (2.6) is given by [1, 2] : R?1 a(!) (2.9) h! = a (!)R?1 a(!) assuming that R is invertible. Of course, (2.9) is not directly implementable since R is not available. The most commonly used implementable form is obtained by using the sample covariance matrix, M 1 X y(t)y (t) (2.10) R^ = M t=1
in lieu of R in (2.9): R^ ?1a(!) hCapon = (2.11) ! a (!)R^ ?1 a(!) The Fourier ([6, 7, 8]) and APES ([5]) spectral estimators have not been derived in the lter bank approach framework. However they can be cast into that framework by simple manipulations. The corresponding FIR lter vectors turn out to be (see, [6, 7] and [5], respectively): a(!) hFourier = (2.12) !
and where
L
^ ?1 hAPES = Q ^(?!1)a(!) ! a (!)Q (!)a(!)
(2.13)
Q^ (!) = R^ ? Y (!)Y (!)
(2.14)
3
and where Y (!) =
M 1 X y(t)e?i!t M
(2.15)
t=1
In the next section we propose the use of the matched lter (MAFI) approach to derive the lter in Step (a) of the lter bank-based spectral estimation methodology outlined above. Even though the matched lter should be a natural choice for Step (a), apparently its use in the previous context was not directly suggested before.
3 The matched- lter bank approach By making use of the notation h
n(t) = (t) : : : (t + L ? 1)
and
iT
(3.1)
s(t) = (!)ei!t a(!)
(3.2)
we can write y(t) in (2.4) as: y(t) = s(t) + n(t) (3.3) We assume that the initial phase of the sinusoidal signal in the equation above is a random variable which is uniformly distributed over the interval (0; 2] and is independent of the noise term. Then it follows that the covariance matrix of the signal term in (3.3) is given by: E [s(t)s (t)] = j (!)j2 a(!)a (!) (3.4) Let Q(!) denote the covariance matrix (assumed to be invertible) of the noise term in Equation (3.3): Q(!) = E [n(t)n (t)] (3.5) The notation Q(!) is meant to indicate the fact that the noise (or residual) term in (3.3), and hence its covariance, depends on !. By de nition the matched lter should be designed such that the corresponding signal-to-noise ratio (SNR) in the lter's output is maximized [9]: jh! a(!)j2 (3.6) max h! h! Q(!)h! The solution to the design problem above is well known [6, 7, 8] (for completeness we provide a simplest derivation in Appendix A: Q?1 (!)a(!) (3.7) h! = a (!)Q?1 (!)a(!) The covariance matrix in (3.7) can be estimated in several ways. As shown below two natural estimators of Q(!) lead to the Capon and APES lters, respectively. A third estimate of Q(!) yields a new lter that is di erent from both the Capon and APES lters. However, as shown in the next section, the spectral estimator corresponding to the new lter turns out to be equivalent to APES as well.
3.1 Capon Since
R = j (!)j2 a(!)a (!) + Q(!)
(3.8)
we can estimate Q(!) as: Q^ (!) = R^ ? j ^(!)j2 a(!)a (!)
(3.9)
4
where ^(!) is an estimate of (!). As we show below the second term in (3.9) has no in uence on the h! in (3.7). This means that for the Q^ (!) in (3.9) there will be in fact no need to estimate (!). A straightforward use of the matrix inversion lemma (stated in Appendix B for easy reference) yields1 : 2 ^ ?1 ^ ?1 R^ ?1 a (3.10) (R^ ? j ^j2 aa )?1 a = R^ ?1a + j ^j R 2aa ^R?1 a = 1 ? j ^j a R a 1 ? j ^j2 a R^ ?1 a Since the scaling of Q^ ?1 a in the expression of the sample counterpart of h(!) in (3.7) leaves the FIR vector unchanged, it follows from (3.10) that the estimated matched lter corresponding to the Q^ in (3.9) coincides with the Capon lter. A perhaps even more direct way to see that the Capon lter is a matched lter is as follows. Consider the MAFI design problem in (3.6). As the multiplication of h! by any complex number leaves the ratio in (3.6) unchanged, it follows that the optimization problem in (3.6) is equivalent to: min h Qh subject to h! a(!) = 1 (3.11) h ! ! !
However, by using (3.8) along with the constraint in (3.11), we can write: h i h! Qh! = h! R ? j (!)j2 a(!)a (!) h! = h! Rh! ? j (!)j2 The conclusion is that the Capon design problem is nothing but a matched- lter design.
(3.12)
3.2 APES By combining (3.2) and (3.3) we obtain: y(t) = [ (!)a(!)] ei!t + n(t) (3.13) The least-squares (LS) estimate of the vector (!)a(!) in (3.13), which ignores the fact that a(!) is known, is given by the normalized Fourier transform (see, e.g., [10]): M X 4 Y (!) y(t)e?i!t = [ (!d )a(!)] = M1
(3.14)
t=1
Inserting (3.14) into (3.9) yields the following estimate of Q(!): Q^ (!) = R^ ? Y (!)Y (!) the use of which leads to the APES lter in (2.13), (2.14).
(3.15)
3.3 MAFI Equation (3.13) suggests another way to estimate the covariance matrix Q(!): Q^ (!) =
M 1 X y(t) ? ^(!)a(!)ei!t y(t) ? ^(!)a(!)ei!t M t=1
(3.16)
where, as before, ^(!) denotes an estimate of (!). A simple calculation shows that the previous Q^ (!) matrix can be rewritten as follows: Q^ (!) = R^ ? ^ (!)Y (!)a (!) ? ^(!)a(!)Y (!) + j ^(!)j2 a(!)a (!) (3.17) By a calculation similar to (3.10) one can verify that the last term in (3.17) can be dropped when evaluating the MAFI vector. Below we make use of the matrix inversion lemma once more to show that 1
For notational convenience we omit the dependence on !.
5
also the third term in (3.17) can be omitted. Let be the matrix made from the rst two terms in (3.17). Then, by the matrix inversion lemma, ?1 ?1 ?1 (3.18) ( ? ^aY )?1 a = ?1 + 1^ ? ^YaY
?1 a a = 1 ? ^ Y a ?1 a where the dependence on ! has been omitted to simplify the notation. A further use of the matrix inversion lemma yields
?1 a = (R^ ? ^ Y a )?1 a = # " ^ R^ ?1 Y a R^ ?1 ? 1 ^ a= = R + 1 ? ^ a R^ ?1 Y ^ ?1 ^ ?1 ^ ?1 ^ ?1 ^ ?1 = R a ? ^ (a R Y )R a^ ?+1 ^ (a R a)R Y ; (3.19) 1 ? ^ a R Y which gives the following expression for the MAFI vector: h i ^ ?1a(!) + ^ (!) (a (!)R^ ?1 a(!))R^?1 Y (!) ? (a (!)R^ ?1 Y (!))R^ ?1 a(!) R (3.20) hMAFI = ! a (!)R^ ?1 a(!) The previous lter is generally di erent from both the Capon and the APES lters. In e ect, neither of the latter two depends on an estimate of (!), unlike the former. (The estimation of (!), needed by (3.20), is discussed in the next section). Despite this fact, in the next section we show the somewhat surprising result that a most reasonable implementation of the spectral estimator corresponding to (3.20) is identical to APES!
4 Spectrum estimation It is readily checked that all lter vectors previously discussed satisfy: h! a(!) = 1 By making use of this observation and of (3.13) we obtain: h! y(t) = (!)ei!t + h! n(t) The LS estimator of (!) in the above equation is given by (see, e.g., [10]) M 1 X ^(!) = h! y(t)e?i!t = h! Y (!) M t=1
(4.1) (4.2) (4.3)
The Capon and APES estimates of the complex spectrum are obtained by inserting the corresponding lter vector expressions in (4.3). The use of the MAFI lter derived in the previous section for spectrum estimation requires an initial estimate of (!). That estimate might be obtained by any of the other methods discussed herein (i.e., Fourier, Capon or APES), but proceeding in that way would be computationally unacceptable. Actually, we can avoid the need for an initial estimate of (!) in the following way. Inserting (3.20) into (4.3) yields
2
^(a R^ ?1 a) = a R^ ?1 Y + ^ (a R^ ?1 a)(Y R^ ?1 Y ) ? a R^ ?1 Y
from which we obtain (by assuming that the ^'s in the right-hand and left-hand sides are identical): a (!)R^ ?1 Y (!) (4.4) ^MAFI (!) = 2 ? 1 ? 1 ? 1 ? 1 ^ ^ ^ ^ a (!)R a(!) ? (a (!)R a(!))(Y (!)R Y (!)) ? a (!)R Y (!)
6
Next we show that the previous spectral estimator coincides with APES: ^APES (!) = ^MAFI (!) (4.5) The APES estimate of the complex spectrum is given by (once again we omit the dependence on ! for notational convenience): a (R^ ? Y Y )?1 Y ^APES (!) = ^ = a (R ? Y Y )?1 a ?1 ?1 a R^ ?1 + R^1?YYYR^?R^1 Y Y = = a R^ ?1 + R^1??1YYYR^?R^1?Y1 a a R^ ?1 Y = (4.6) ^ ?1 2 ? 1 ? 1 ? 1 ^ ^ ^ a R a ? (a R a)(Y R Y ) ? a R Y which is identical to (4.4). Hence the equality in (4.5) is proved. For comparative purposes, note that the Capon spectral estimator is given by: a (!)R^?1 Y (!) ^Capon (!) = ^ ?1 a (!)R a(!) The next section discusses the relative merits of Capon and APES.
(4.7)
5 Comparisons and discussions 5.1 The computational issue
Let R^ ?1=2 denote a square root of the positive de nite matrix R^ ?1 , and let: (!) = R^ ?1=2 a(!) (!) = R^ ?1=2 Y (!)
The Capon and APES spectral estimators can be expressed as relatively simple functions of the (!) and (!) introduced above: (!)(!) 4 (!)(!) ^Capon (!) = (5.1) = (!) (!) k (!)k2 (!)(!) (5.2) ^APES (!) = k (!)k2 ? k (!)k2 k(!)k2 ? j (!)(!)j2 The ecient estimation of the Capon spectrum estimate is discussed, for example, in [6, 7, 8]. We do not intend to dwell into the details of that computation. The point we want to make here is that APES is only slightly more involved computationally than Capon (as should be evident by comparing (5.1) and (5.2)), and hence the ecient computational means developed for the latter can also be applied to the former.
5.2 The statistical performance issue The Capon, APES as well as the Fourier spectral estimators can be shown to have the same asymptotic variance under the following condition:
7
C: The signal x(t) can be written as in (2.1), where (t) is a zero-mean stationary random process with nite spectral density at !: " (!) < 1
(5.3)
In more exact terms, the following result holds true. Theorem 1 Under condition C and the additional assumption that (t) is circularly symmetrically distributed, the estimation errors in the Fourier, Capon and APES spectral estimators are asymptotically circularly symmetrically distributed with zero-mean and the following common variance: lim ME j ^(!) ? (!)j2 = (!) (5.4) M !1 Proof: See Appendix C. The need to enforce condition C limits, to some extent, the importance of the previous result. Indeed the assumption made in C is satis ed if (and essentially only if) the signal x(t) has a mixed spectrum and ! is the location of a spectral line. In some applications, such as target feature extraction by means of a Doppler radar, the signal has a mixed spectrum and estimating the amplitudes of the sinusoidal components is a key problem (see, e.g., [5] and the references therein). The result of Theorem 1 is relevant to such applications. In other applications, however, the main interest is in the continuous component of the spectrum. Condition C does not hold exactly if the spectrum is continuous at !. The previous result is of a somewhat limited interest also because of its asymptotic character. Indeed, in simulations with medium or small-sized samples, the spectral estimators under study have been found to behave quite di erently from one another in contradiction with what is predicted by the (asymptotic) result of Theorem 1 (see, e.g., [5]; also see the next section). In particular, in scenarios with multiple, closely spaced spectral lines the Fourier estimate of the complex spectrum has been found to be (sometimes signi cantly) less accurate than the spectrum estimates obtained by using the Capon and APES methods [5, 3]. The nite-sample analysis of the spectral estimators under discussion would be of considerable interest. However, while this is possible for the Fourier estimator, a complete analysis of the more interesting, data-adaptive Capon and APES estimators appears to be quite dicult at best. A partial analysis of the latter estimators is nevertheless possible as shown in what follows. Speci cally, in Appendix D we make use of a higher-order Taylor series expansion technique to prove that, to within a second-order approximation, and under the additional mild assumption that the third-order moments of n(t) are zero, Capon is biased downwards whereas APES is unbiased: ME [^ Capon (!) ? (!)] < 0 (! 2 [?; )) (5.5) (!) and ME [^ APES (!) ? (!)] = 0 (! 2 [?; )) (5.6) for suciently large values of M . The above result provides theoretical support to the empirically observed fact that Capon underestimates the complex spectrum, whereas the APES estimate is almost unbiased (see [5]). In fact the Capon and APES spectral estimators satisfy the following neat inequality, which points to the type of result proved above: j ^Capon (!)j j ^APES (!)j and arg(^ Capon (!)) = arg(^ APES (!)) , 8! 2 [?; ) (5.7)
8
where the modulus equality holds only for M ! 1. To check (5.7) observe that the term between square parentheses in (4.6) is positive (by the Cauchy-Schwartz inequality) and that it can be zero if and only if Y (!) is proportional to a(!) (which can hold true only when M ! 1). We believe that (5.7) along with (5.5) and (5.6) provide a theoretical motivation for preferring APES to Capon in most spectrum estimation exercises. The simulations in the next section lend further evidence to the fact that APES outperforms Capon in the nite-sample case.
6 Numerical examples In this section we study the accuracy of the APES and Capon amplitude estimates in a number of cases.
6.1 Estimation performance versus SNR We begin by studying the performance of using the APES and Capon methods for complex amplitude estimation as the SNR varies. Figure 1 shows the modulus of the true spectrum of the sinusoidal signal, which consists of four dominant spectral lines and nine small spectral lines located at the following frequencies: 0.0625, 0.0875, 0.25, 0.285, 0.33, 0.35, 0.37, 0.39, 0.41, 0.43, 0.45, 0.47, 0.49. These spectral lines all have a phase value of =4. The data sequence has N = 64 data samples and is corrupted by complex white Gaussian noise with zero mean and variance 2 . The SNR for the kth spectral line is de ned as
(6.1) SNRk = 10 log10 j k2j [dB] where k is the complex amplitude of the kth spectral line. In the following, we compare the biases and variances of the APES (solid lines in the gures) and Capon (dashed lines in the gures) estimates of the complex amplitude at the frequency of the rst dominant spectral line, as its SNR varies. Our results are obtained from 100 Monte Carlo trials in which we used 100 independent sequences of circular Gaussian white noise to generate the simulation data set. 2
Figures 2(a) and 2(b) compare the real and imaginary parts, respectively, of the biases of the Capon and APES estimators when L = 24. As seen from the gures, the APES estimates are nearly unbiased, while the Capon estimates are biased downward, as was expected from our previous theoretical analysis. Figure 2(c) and 2(d) compare the variances and the MSE, respectively, of the Capon and APES estimators when L = 24. Note that the APES estimator provides not only a smaller bias but also a lower variance than the Capon estimator.
6.2 Estimation performance versus lter length Next we study the e ect of the lter length, L, on the performance of the two estimators. The data studied is the same as in the previous example, with the SNR of the rst dominant spectral line xed at 20 dB. Our results are obtained from 100 Monte Carlo trials. Figures 3(a) and 3(b) compare the real and imaginary parts, respectively, of the biases of the Capon and APES estimators as L varies. As seen from the gures, the APES estimates are again almost unbiased for not too small lter lengths, while the Capon estimates have very large biases for large L. Figure 3(c) and 3(d) compares the variances and the MSE, respectively, of the Capon and APES estimators. Note that the variance of the APES estimate varies little for the studied lter lengths, while the variance of the Capon estimate signi cantly increases with increasing L. This type of behavior makes the
9
Modulus of True Spectrum
1
0.8
0.6
0.4
0.2
0 0
0.2
0.4
0.6
0.8
1
Frequency
Figure 1: The modulus of the true spectrum. choice of L with the Capon method quite dicult. We remind the reader that to enhance the resolution L should be chosen as large as possible (subject to L < N=2 so that R^ is not singular). While APES with this choice of L works ne, Capon does not. We further illustrate this point in the next gure. In Figures 4(a) and 4(b), one realization of the modulus of the estimated spectra obtained by using APES and Capon, respectively, is compared with the true spectrum (dotted lines) for L = 15. Figures 4(c) and 4(d) are similar to Figures 4(a) and 4(b) except that the lter length is now increased to L = 24. As can be seen from the gures, APES gives much more accurate peak amplitude estimates, at the price of slightly wider spectral peaks. Note that, as L increases, the resolution of both methods improves, as expected; however, unlike APES, the performance of the Capon amplitude estimates detoriates significantly with increasing L. Note also that the APES estimates, at all evaluated frequencies, are larger than the Capon estimates, as expected. In Table 1 we show the phase estimates at the frequencies of the four dominant spectral lines. In accordance with what was expected from (5.7), the APES and Capon phase estimates coincide.
10
0.02
0.02
0
0
−0.02
APES Capon Bias
Bias
−0.02
−0.04
APES Capon −0.04
−0.06
−0.06
−0.08
−0.08
−0.1
2
4
6
8
10
12
14
16
18
−0.1
20
2
4
6
8
10
SNR
12
14
16
18
20
SNR
(a)
(b)
−1
−1
10
10 APES Capon
APES Capon
−2
−2
10 MSE
Variance
10
−3
−3
10
10
−4
10
−4
2
4
6
8
10
12
14
16
18
10
20
SNR
2
4
6
8
10
12
14
16
18
20
SNR
(c)
(d)
Figure 2: Bias and variance of the complex amplitude estimate at the frequency of the rst spectral line as a function of SNR when L = 24. (a) Real part of the bias. (b) Imaginary part of bias. (c) Variance. (d) MSE.
11
0.15
0.15 APES Capon
0.05
0.05
0
0
−0.05
−0.05
−0.1
−0.1
−0.15
−0.15
−0.2
−0.2
−0.25
−0.25
−0.3
−0.3 5
10
15
20
25
APES Capon
0.1
Bias
Bias
0.1
30
5
10
15
L
20
25
(a)
(b)
0
0
10
10 APES Capon
APES Capon
−1
−1
10
MSE
Variance
10
−2
10
−3
−2
10
−3
10
10
−4
10
30
L
−4
5
10
15
20
25
10
30
L
5
10
15
20
25
30
L
(c)
(d)
Figure 3: Bias and variance of the complex amplitude estimate at the frequency of the rst spectral line as a function of L when the SNR of the rst spectral line is 20 dB. (a) Real part of the bias. (b) Imaginary part of bias. (c) Variance. (d) MSE.
12
1 Modulus of Estimated Spectrum
Modulus of Estimated Spectrum
1
0.8
0.6
0.4
0.2
0 0
0.2
0.4
0.6
0.8
0.4
0.2
0.2
0.4
Frequency
Frequency
(a)
(b)
0.6
0.8
1
0.6
0.8
1
1 Modulus of Estimated Spectrum
Modulus of Estimated Spectrum
0.6
0 0
1
1
0.8
0.6
0.4
0.2
0 0
0.8
0.2
0.4
0.6
0.8
1
0.8
0.6
0.4
0.2
0 0
0.2
0.4
Frequency
Frequency
(c)
(d)
Figure 4: Modulus of the estimated spectrum (solid lines) compared with the true spectrum (dotted lines) when the SNR of the rst spectral line is 20 dB. (a) APES with L = 15. (b) Capon with L = 15. (c) APES with L = 24. (d) Capon with L = 24.
13
Phase of spectral line no. 1 2 3 4
True Estimated values values APES, L=15 Capon, L=15 APES, L=24 Capon, L=24 0:25 0:2457 0:2457 0:2477 0:2477 0:25 0:3506 0:3506 0:3408 0:3408 0:25 0:2513 0:2513 0:2507 0:2507 0:25 0:2493 0:2493 0:2454 0:2454
Table 1: Phase estimates at the frequencies of the four dominant spectral lines.
7 Conclusions The Capon spectral estimator has been shown to be a matched- lter bank method. In fact, the Capon's lter design was shown to be precisely a matched lter design. The APES spectral estimator has also been shown to be a member of the class of matched- lter bank methods. In e ect, APES was proved to coincide with two seemingly di erent matched- lter bank spectral estimation methods. The matched- lter bank interpretation of Capon and APES provides a (partial) theoretical explanation for the good, empirically observed performance of these two spectral estimators. For the moment it is an open issue whether other matched- lter bank spectral estimation methods with enhanced properties exist. The matched- lter bank interpretation has also provided further insight into the properties of and relationships between Capon and APES. In particular the signal model used by the two methods was shown to be adequate only at the frequencies corresponding to the spectral lines in a mixed spectrum. The model in question is not valid at the frequencies of the continuous spectral component. While APES was shown by means of simulations to yield accurate estimates of the latter component as well [5], there is a clear need for a better theoretical understanding of why and how APES works in such cases (The performance of the Capon estimate of the continuous spectral component was poor in our simulations reported in [5]). Derivation of matched- lter bank methods speci cally designed for continuous spectra is also an interesting open issue. The casting of Capon and APES into the matched- lter bank framework, along with the related analysis, has also led to a neat result concerning the relationship between these two methods: the Capon spectral estimate is always smaller than APES's, for all frequencies. This theoretical observation lends support to the empirical nding that Capon often underestimates the true spectrum whereas APES is nearly unbiased. In fact by making use of an expansion technique we were able to prove that, to within a second-order approximation, Capon is indeed biased downwards whereas APES is unbiased. As APES is only slightly more involved than Capon, from a computational standpoint, preferring the former to the latter appears to be a logical conclusion of both the theoretical and empirical analysis in this paper.
14
Appendix A. Proof of (3.7) Let Q1=2 denote a Hermitian (for notational simplicity) square root of Q. Then, by the Cauchy-Schwartz inequality, h Q1=2 Q?1=2 a(! ) 2 ! a (!)Q?1 a(!) h Qh !
!
where the equality is achieved if and only if h! = Q?1 a(!) for any scalar 2 C. As is made transparent in Section 4, it is often required that h! passes the frequency ! undistorted, that is h! a(!) = 1. This requirement can be met by choosing as = ? 1 1= a (!)Q a(!) , which yields (3.7).
Appendix B. The matrix inversion lemma Let R be an n n matrix, and let b and c be n 1 vectors. Then, under the assumption that the inverse appearing below exists, ?1 ?1
bc R (R ? bc )?1 = R?1 + R 1 ? c R?1 b
(see, e.g., [6, 10]).
Appendix C. Proof of Theorem 1 To simplify the notation we omit the dependence on ! whenever there is no possibility for confusion. From (3.13) we have that Y = a + (C.1) where M X = 1 n(t)e?i!t (C.2) M t=1
The rst and second-order moments of are easy to calculate: E () = 0 E (T ) = 0
(C.3) (C.4)
(by the circularly symmetric distribution assumption) and E ( ) =
M M X M 1 X ?i!(t?s) = 1 X (M ? j j)Rn ( )e?i! R ( t ? s ) e n M2 M2 =?M
t=1 s=1
(C.5)
where fRn ( )g is the covariance sequence of n(t). It follows from (C.5) and condition C that: 1 X lim ME ( ) = Rn ( )e?i! = n (!) = (!)aa M !1 =?1
(C.6)
where the last equality follows from standard results on the transfer of spectral densities through linear systems. Among others, the previous calculations imply that, as M ! 1, Y tends to a (in the mean square sense). Hence hCapon and hAPES have the same limit as M increases without bound. ! ! Let h denote a generic FIR vector and let h1 denote the deterministic vector that is the limit of (the
15
possibly random) h when M goes to in nity. Observe that for all methods under study the associated h and h1 vectors satisfy: h a = 1 and h1 a = 1 (C.7) By using this observation along with (C.1) we obtain: 4 h Y = + h ^ =
(C.8) Because tends to zero as M ! 1 (as shown previously), it follows from (C.8) that the estimation error in ^ can asymptotically be written as (to within a rst-order approximation): ^ ? ' h1 (C.9) From (C.9) and (C.3)-(C.6) we readily derive that ME (^ ? )2 = 0 as M ! 1 and lim ME j ^ ? j2 = h1 Mlim ME ( )h1 = (!) jh1 aj2 = (!) M !1 !1 and the proof is concluded.
Appendix D. Proof of (5.5), (5.6) Once again we omit the dependence on ! for notational convenience.
Proof of (5.5) By using (2.9) and (C.8) we obtain a R^ ?1 (D.1) ^Capon ? = ^ ?1 aR a In what follows we use the symbol to denote an \asymptotic equality" that holds to within a secondorder approximation. A straightforward manipulation of (D.1) yields: a (R^ ?1 ? R?1) a R?1 + ^ ?1 = ^Capon ? = a R^ ?1 a aR a a R^ ?1 (R ? R^ )R?1 ?1 1 ? 1 + 1 = + a R a R^ ?1 a a R^ ?1 a a R?1a a R?1 a ?1 ^ ?1 ?1 ?1 ?1 ^ ?1 ? a R (aR R??R1 a)R + aaRR?1a + a R ^ ?a1 (R ??R1 )a (a R a)(a R a) ? 1 ?1 ?1 ^ ?1 ? 1 ^ ? 1 ? a R (aR R??R1 a)R + aaRR?1a + (a R )(aaRR?1(aR)2? R)R a which, in turn, implies that " # ^ ?1 )(a R?1 a) + (a R?1 RR ^ ?1 a)(a R?1 ) ? ( a R?1 RR E [^ Capon ? ] E (D.2) (a R?1 a)2 Next we note that M X 1 ^ R = aei!t + n(t) a e?i!t + n (t) = M t=1
M X n(t)n (t) = j j2 aa + a + a + M1 t=1
16
and also remind the reader of the de nition (C.2) of and the assumption that n(t) has zero third-order moments. By making use of these facts, along with (D.2), we can write: (a R?1 a)2 E (^ Capon ? ) E ?(a R?1 a)(a R?1 a R?1 + a R?1 a R?1 )+ + (a R?1 )( a R?1 a R?1a + a R?1 a R?1 a) = = E ? (a R?1 a)(a R?1 )2 ? (a R?1 a)2 ( R?1 )+ o + (a R?1 )2 (a R?1 a) + (a R?1 a) a R?1 2 = n
o
= (a R?1 a)E a R?1 2 ? (a R?1a)( R?1 ) (D.3) The quantity between curly parentheses in (D.3) is negative and thus is its expectation (by the CauchySchwartz inequality). Hence (5.5) follows.
Proof of (5.6) The APES estimate of the complex spectrum is given by (cf. (4.5)): a R^ ?1 Y ^APES = ^ ?1 a R a? where 2 = (a R^ ?1a)(Y R^ ?1 Y ) ? a R^ ?1Y
(D.4) (D.5)
By using (C.1) in (D.5) we obtain: = (a R^ ?1a) j j2 a R^ ?1a + a R^ ?1 + R^ ?1a + R^ ?1 ? j j2 (a R^ ?1 a)2 ?
2
? (a R^ ? a)( R^ ? a) ? (a R^ ? )(a R^ ? a) ? a R^ ? = 1
1
1
1
1
2 4 = (a R^ ?1a)( R^ ?1 ) ? a R^ ?1 (a R?1a)( R?1 ) ? a R?1 2 = Next note that 1 = ^1?1 + ^ ?11 ? ^1?1 = ^1?1 + ^ ?1 ^ ?1 ? 1 ^ a R a? a R a a R a ? a R a a R a (a R a)(a R a ? ) Consequently, a R^ ?1 ( a + ) ? (a R^ ?1 a) + a R^ ?1 + ^APES ? = = ^ ?1 = a R^ ?1 a ? a R a? ^ ?1 ^ ?1 = a R ^?1+ + ^(?a1 R ^+?1 ) aR a (a R a)(a R a ? ) ? 1 ^ (D.6) a R^ ?1 + a R a a R^ ?1 a The rst term in (D.6) is equal to (^ Capon ? ) and hence it follows from the rst part of this appendix that the expectation of this term satis es (to within a second-order approximation): # " ? E () a R^ ?1 (D.7) E ^ ?1 ^ ?1 aR a aR a
Combining (D.6) and (D.7) yields: E (^ APES ? ) 0 and the proof is almost concluded. It only remains to motivate the normalizing factor M used in both (5.5) and (5.6). However the factor in question follows relatively easily from the fact that is O(1=M 1=2 ) (in the mean square sense) and this implies that the second-order approximations previously used for both E (^ Capon ? ) and E (^ APES ? ) are O(1=M ).
17
References [1] J. Capon, \High Resolution Frequency Wave Number Spectrum Analysis", Proc. IEEE, 57:1408{ 1418, 1969. [2] R. Lacoss, \Data adaptive spectral analysis methods", Geophysics, 36:134{148, 1971. [3] M.A. Lagunas, M.E. Santamaria, A. Gasull, and A. Moreno, \Maximum likelihood lters in spectral estimation problems", Signal Processing, 10(1):19{34, 1986. [4] M. Wax, \Model-based processing in sensor arrays", In S. Haykin, editor, Advances in Spectral Analysis and Array Processing, volume III, pages 1{47. Prentice-Hall, Englewood Cli s, N.J., 1995. [5] J. Li and P. Stoica, \Adaptive Filtering Approach to Spectral Estimation and SAR Imaging", IEEE Trans. on Sig. Proc., 44(6):1469{1484, June 1996. [6] S.M. Kay, Modern Spectral Estimation: Theory and Application, Prentice-Hall, Englewood Cli s, N.J., 1988. [7] S.L. Marple, Jr., Digital Spectral Analysis, Prentice-Hall, Englewood Cli s, N.J., 1987. [8] P. Stoica and R. Moses, Introduction to Spectral Analysis, Prentice Hall, Upper Saddle River, NJ, 1997. [9] B. Picinbono, Random signals and systems, Prentice-Hall, Englewood Cli s, N.J., 1993. [10] T. Soderstrom and P. Stoica, System Identi cation, Prentice Hall International, London, UK, 1989.
18
