Frequency Estimation and the QD Method

Richard M. Todd; J.R. Cruz    School of Electrical Engineering, The University of Oklahoma, Norman, OK 73019

I Introduction to Frequency Estimation

In this chapter we discuss some recently developed algorithms for frequency estimation, that is, taking an input signal assumed to be a combination of several sinusoids and determining the frequencies of these sinusoids. Here we briefly review the better known approaches to frequency estimation, and some of the mathematics behind them.

I.A Periodograms and Blackman-Tukey Methods

The original periodogram-based frequency estimation algorithm is fairly simple. Suppose we have as input a regularly sampled signal x[n] (where n varies from 0 to N − 1. so there are N data points in total). To compute the periodogram, one simply takes the signal, applies a Discrete Fourier Transform, and take the magnitude squared of the result. This gives one an estimate of the power spectral density (PSD) ^xxf of the input signal:

^xxf=1N∑n=0N−1xne−j2πfn2

where j is the square root of minus one. If the input signal is assumed to be a sum of sinusoids, one can then derive estimates for the frequencies of these sinusoids from the peaks of the ^xxf estimate as a function of frequency — each peak will correspond to an input sinusoid.

This approach to frequency estimation seems simple enough, but it has some problems. One would expect that the variance of this estimate would go to zero as the number of input data samples available increases, i.e, as the number of samples increases, the quality of the estimate of the spectral density would get better. Surpringly, this does not happen; to first order, the variance of ^xxf does not depend on N (see [1] for details). Intuitively, this lack of decrease of the variance can be explained as follows: the spectral estimate ^xxf can be looked at as

^xxf=N∑k=0N−1h−kxk2

where

k=1Nexp−j2πfkk=−N−1,…,0

turns out to be the impulse response of a bandpass filter centered around f. So the spectral estimate turns out to be, basically, a single sample of the output of a bandpass filter centered at f. Since only one sample of the output goes into the computation of the estimate, there is no opportunity for averaging to lower the variance of the estimate.

This suggests a simple way to improve the variance of the periodogram: split the input signal into M separate pieces of length N/M, compute the periodogram of each one, and average them together. This does reduce the variance of the spectral estimates by a factor of 1/M. Alas, this benefit does not come without a cost. It turns out that the spectral estimates from the periodogram not only have a (rather nasty) variance, they also have a bias; the expected value of the estimate turns out to be the true power spectral density convolved with WB(f), the Fourier transform of the Bartlett window function

bn=1−kN|k|<N−1

As N grows larger, WB(f) approaches a delta function, so the periodogram approaches an unbiased estimate as N → ∞. But when we split the signal into M pieces of size N/M, we are now computing periodograms on segments 1/Mth the size of the original, so the resulting periodograms have considerably worse bias than the original. We have improved the variance of the periodogram by splitting the signal into M pieces, but at a cost of significantly increasing the bias, increasing the blurring effect of convolving the true PSD with the WB(f) function. As the variance improves, the resolution of the algorithm, the ability for it to detect two closely spaced sinusoids as being two separate sinusoids, goes down as well; with smaller segment sizes, the two peaks of the true PSD get blurred into one peak.

Blackman and Tukey invented another modification of the periodogram estimator. One can readily show that the original periodogram

can be written as

^xxf=∑k=−N−1N−1r^xxkexp−j2πfk

where

^xxk=1N∑n=0N−1−kx*nxn+kfork=0,…,N–1r^xx*−kfork=−N−1,…,−1

is a (biased) estimate of the autocorrelation function of the signal. Hence, the periodogram can be thought of as the Discrete Fourier Transform of the estimated autocorrelation function of the signal. The Blackman-Tukey algorithm simply modifies the periodogram by multiplying the autocorrelation by a suitable window function before taking the Fourier Transform, thus giving more weight to those autocorrelation estimates in the center and less weight to those out at the ends (near ± N), where the estimate depends on relatively few input sample values. As is shown in Kay [1], this results in a trade-off similar to that involved in the segmented periodogram algorithm; the variance improves, but at the expense of worsened bias and poorer resolution.

I.B Linear Prediction Methods

A wide variety of frequency estimation algorithms are based on the idea of linear prediction. Linear prediction, basically, is assuming that one’s signal satisfies some sort of linear difference equation, and using this model to further work with the signal. We will primarily concentrate in this section on the so-called AR (auto-regressive) models, as opposed to the MA (moving-average) model, or the ARMA model, which is a combination of AR and MA models. This is because, as we will show, the AR model is particularly applicable to the case we are interested in, the case of a sum of sinusoids embedded in noise.

An ARMA random process of orders l and m is a series of numbers y[n] which satisfy the recurrence relation

n=∑i=0mbien−i−∑i=1laiyn−i

where e[n – i] is a sequence of white Gaussian noise with zero mean and some known variance. The ai are called the auto-regressive coefficients, because they express the dependence of y[n] on previous values of y. The bi coefficients are called the moving average coefficients, because they produce a moving average of the Gaussian noise process. One can consider this process as being the sum of the output of two digital filters, one filtering the noise values e[n], and one providing feedback by operating on earlier values of the output of the ARMA process.

A MA process is just a special case of the ARMA process with a1 = … = al = 0, i.e.,

n=∑i=0mbien−i

(note that without loss of generality, we can take b0 = 1 by absorbing that term into the variance of the noise process e[n]). Similarly, an AR process is just the special case with bı = … = bm = 0, so

n=−∑i=0laiyn−i+en

(note that, again, we can always take b0 = 1.)

As mentioned above, the AR model is particularly suitable to describing a sum of sinusoids; we now show why this is so. Consider the AR process described by Eq. (9), and consider what happens when we filter it as follows:

n=yn+a1yn−1+…+a1yn−l,

that is to say, run it through a digital filter with coefficients 1, a1,…, al. Looking at the definition of the AR process y[n], one can readily see that the output x[n] of our filter is nothing but the noise sequence:

n=en∀n.

The above digital filter is called the predictor error filter, because it is the error between the actual signal y[n] and a prediction of that signal based on the previous m values of the signal. It can be shown [1] that this particular filter is optimal in the sense that if you consider all possible filters of order m which have their first coefficient set to 1, the one that produces an output signal of least total power is the one based on the AR coefficients of the original AR process. Furthermore, the optimal filter is the one which makes the output signal white noise, removing all traces of frequency dependence in the spectrum of the output signal. Now suppose we have a signal composed of a sum of m complex sinusoids plus some white Gaussian...