High-resolution frequency estimation

Two tones closer than one bin of a 64-point DFT blur into one lump. Here three methods read the same 64 samples in turn: the periodogram, Burg’s AR model and MUSIC. Watch the diamonds at the peaks against the two dotted lines, where the true tones are.

Closer than one bin

Two stars close together in the night sky can look like one blur in a small telescope. A bigger telescope splits them. But if you know that each star is a point of light, you can fit two points to the blur and read where they sit. This page does the same for two tones.

Here is the record. I take $N=64$ samples of

where $v[n]$ is white noise of nominal RMS 0.1, as in Random processes (24.2). It is a seeded draw (seed 254), so it is the same on every visit. This draw measures RMS 0.0950, a power of 0.0090. Each tone has power 0.5, so each stands 17.4 dB above the noise.

One bin of a 64-point DFT is $2\pi/N=2\pi/64=0.03125\pi$. The tones are $0.03\pi$ apart, which is 0.96 of a bin. In “Two tones 100 Hz apart, more and more bins” of The DFT (13.2), two tones needed a bin spacing no larger than their gap. Here that asks for about $2\pi/(0.03\pi)\approx67$ samples, and I have 64.

Padding does not rescue me. In “Padding cannot split two tones” of Zero-padding and resolution (15.3), each tone left a lump as wide as the record allows, and padding only drew that lump more finely. So the periodogram $P_N(e^{j\Omega})=\frac1N\lvert X(e^{j\Omega})\rvert^2$ of The periodogram (25.1), computed at 8193 frequencies from 0 to π, shows one lump. Its top is at 0.221π, between the two tones.

A periodogram assumes nothing about the signal, and that is also its limit. The two methods on this page assume a model: a known number of tones in white noise. In return they can split tones closer than one bin. From statistics I assume only what 24.2 used: white noise is uncorrelated from one sample to the next, so $R_v[\ell]=\sigma_v^2\,\delta[\ell]$.

Poles instead of bins

In “A smooth model against a ragged periodogram” of Parametric models and linear prediction (25.3), an AR model gave the spectrum $E_p/\lvert A_p(e^{j\Omega})\rvert^2$. It peaks wherever $A_p(e^{j\Omega})$ comes close to 0, which happens next to a pole near the unit circle. A sharp peak needs a pole near the circle, not a long record. So two pairs of poles can make two peaks closer than a bin.

The question is how to find the coefficients from 64 samples. In “One order at a time”, 25.3 estimated the autocorrelation and ran the Levinson–Durbin recursion on it. That estimate divides by $N$ and treats the samples outside the record as zeros, and its DTFT is exactly the periodogram. So it carries the same blur: on this record, order 16 by that route still gives one peak, at 0.223π.

Burg’s method, from J. P. Burg (1967), keeps the recursion and changes one step. It still adds one reflection coefficient $k_m$ per order, with $E_m=E_{m-1}(1-k_m^2)$. But it takes each $k_m$ from the data themselves, not from an autocorrelation.

At order $m$ it predicts each sample from the $m$ samples before it, and also from the $m$ samples after it. Then it picks the $k_m$ that makes the squares of both prediction errors, added up, as small as possible. It only uses products of samples inside the record, so nothing is filled with zeros. And no $\lvert k_m\rvert$ can exceed 1, so the model is stable, by 25.3’s rule.

As a check, on 25.3’s 4096 samples of AR(4) noise, Burg’s order 4 returns the true coefficients to within 0.014.

Tones in noise are not exactly an AR process, so the model needs more than the four poles of two tones. On this record order 8 still gives one peak, at 0.215π. Order 16 gives two, at 0.201π and 0.231π, against the true 0.200π and 0.230π.

Four above the floor

The second method looks at short pieces of the record side by side. Cut the record into overlapping pieces of $m$ samples. Each piece is a snapshot, written as a column:

A snapshot runs forward in time, oldest sample first. With $m=16$ and 64 samples, $n$ goes from 0 to 48, so there are 49 snapshots.

For one snapshot, the products of every pair of its samples fill a 16 × 16 table, $\mathbf{x}_n\mathbf{x}_n^\top$. Average the 49 tables, and you have an estimate of the correlation matrix:

Its entry in row $i$ and column $k$, counting from 0, is an average of $x[n+i]\,x[n+k]$. So it estimates $R_x[i-k]$ of 24.2, as in 25.3’s matrix $\mathbf{R}$, but made directly from products of samples.

A real process has the same correlations played backward, because $R_x[-\ell]=R_x[\ell]$. So the true matrix does not change if I reverse the order of its rows and of its columns. The estimate changes a little, so I average it with its reversed copy, which uses every snapshot twice, once read forward and once backward. This is forward–backward averaging, and from here on $\mathbf{R}$ is the averaged matrix.

Directions that are only stretched

Now one idea from linear algebra. An eigenvector $\mathbf{v}_i$ of $\mathbf{R}$ is a direction that $\mathbf{R}$ only stretches, and the stretch is its eigenvalue $\lambda_i$:

You met eigenvectors in “Arrows in, the same arrows out” of The DFT as a matrix (13.5). The one fact I need is this: a symmetric $m\times m$ matrix such as $\mathbf{R}$ has $m$ eigenvectors that are orthogonal, at right angles to each other, with real eigenvalues. A library routine such as NumPy’s eigh finds them, each scaled to length 1, and I number them from the largest eigenvalue down.

Why should tones stand out among them? A cosine is two complex exponentials, $\cos\theta=\tfrac12(e^{j\theta}+e^{-j\theta})$, as in “Two points spinning opposite ways” of Complex exponentials & phasors (3.4). So two real tones are four exponentials, at $\Omega=\pm0.20\pi$ and $\pm0.23\pi$. A snapshot of $e^{j\Omega n}$ is the number $e^{j\Omega n}$ times one fixed column:

Its conjugate transpose $\mathbf{f}(\Omega)^{\mathsf H}$ is a row of 13.5’s matrix $\mathbf{F}$, with $2\pi k/N$ replaced by any $\Omega$. Call the four frequencies $\Omega_1$ to $\Omega_4$. For a long record, the correlation matrix of the two tones plus white noise is

The $\tfrac14$ is each exponential’s power, the square of its amplitude $\tfrac12$. The noise adds $\sigma_v^2$ on the diagonal only, because white noise is uncorrelated at every lag but 0.

Now take any column $\mathbf{v}$ orthogonal to all four, so that $\mathbf{f}(\Omega_k)^{\mathsf H}\mathbf{v}=0$ for each $k$. Every term of the sum then sends $\mathbf{v}$ to 0, and

So $\mathbf{v}$ is an eigenvector with eigenvalue $\sigma_v^2$. In 16 dimensions, four columns leave room for $16-4=12$ such directions. The other four eigenvectors lie among the exponentials, and their eigenvalues add the tones’ power to $\sigma_v^2$.

The four eigenvectors with the largest eigenvalues span the signal subspace, which means all their combinations. The other twelve span the noise subspace. So I expect four eigenvalues well above the noise power and twelve equal to it.

Here are this record’s eigenvalues in decibels, $10\log_{10}\lambda_i$. No eigenvalue is negative, because $\mathbf{R}$ is an average of tables $\mathbf{x}_n\mathbf{x}_n^\top$.

The four tallest bars are 8.509, 7.426, 0.428 and 0.253. The other twelve range from 0.0025 to 0.0142 and average 0.0101, which is −20.0 dB. That is close to the measured noise power, 0.0090. They are not all equal, because 49 snapshots are a small sample: the exact matrix of the model, with $\sigma_v^2=0.0090$, has all twelve at 0.0090.

Why are the third and fourth so much smaller than the first two? Over 16 samples, one tone slips only $15\times0.03\pi=0.45\pi$ against the other, less than a quarter turn. So the two tones look almost alike, and most of their power lies along two directions they share. The model’s exact matrix agrees: 0.484 and 0.268 for the third and fourth.

This is how to read the number of tones from data. Count the eigenvalues that stand clearly above the floor, and halve the count, because each real tone is two exponentials. Here there are four, so two tones.

Two tones inside one bin

The twelve noise eigenvectors are orthogonal to every tone’s column $\mathbf{f}(\Omega_k)$. Turn this around: try every frequency $\Omega$, and ask how much of $\mathbf{f}(\Omega)$ lies in the noise subspace. Because the $\mathbf{v}_i$ have length 1 and are orthogonal, that part has the squared length

Each term $\mathbf{f}(\Omega)^{\mathsf H}\mathbf{v}_i$ is the DTFT of the 16 entries of $\mathbf{v}_i$, the sum of “Samples, turned and added” in The DTFT (12.2). So each noise eigenvector, used as a filter of 16 taps, nearly blocks both tones. On this record the sum is 0.00015 at 0.20π and 0.00079 at 0.23π, against 12.0 on average from 0 to π.

MUSIC (multiple signal classification, from R. O. Schmidt, 1986) draws the reciprocal:

Where the sum nearly vanishes, the curve shoots up. This curve is a pseudospectrum: its height says how nearly $\mathbf{f}(\Omega)$ misses the noise subspace, not how much power sits at $\Omega$.

The picture at the top of the page shows the three estimates of this one record in turn: the periodogram, Burg’s AR(16), and MUSIC with $m=16$ and a signal subspace of four dimensions. They measure different things, so each curve is drawn in decibels relative to its own peak. After the clip, three buttons switch between them.

Notice where the diamonds land at the end. MUSIC’s peaks sit at 0.201π and 0.232π, against the true 0.200π and 0.230π: off by 0.03 and 0.07 of a bin. Between them its curve falls to −20.9 dB.

MUSIC has two relatives worth naming. Pisarenko’s method (V. F. Pisarenko, 1973) is the smallest case: $m$ is one more than the number of exponentials, so the noise subspace is a single eigenvector. With $m=5$ on this record it finds one tone, at 0.217π, and puts its other two frequencies at 0 and π: five samples are too short to tell these tones apart.

ESPRIT (R. Roy and T. Kailath, 1989) uses the signal subspace instead. Moving a snapshot on by one sample turns each exponential by $e^{j\Omega}$, and ESPRIT reads those turns from the signal eigenvectors. This page leaves its algebra there.

What the model costs

Both methods bought their sharpness with assumptions. MUSIC needs the number of tones, or a clear gap in the eigenvalues to read it from. It needs white noise, because the step to $\mathbf{R}\mathbf{v}=\sigma_v^2\mathbf{v}$ used the noise’s $\sigma_v^2\mathbf{I}$. If the signal is not a few tones in white noise, the model is wrong, and sharp peaks can mislead.

Here is what a wrong count does on this record. Tell MUSIC there is one tone, two signal dimensions, and it finds one peak, at 0.215π. Tell it three tones, six dimensions, and it still finds 0.201π and 0.232π. Here, too few dimensions merged the tones, and two too many did little harm.

Burg’s method has its condition too: the order must be high enough, and at 64 samples order 8 was not.

The periodogram has a trap of its own. The second tone’s phase, 2.0 rad, was chosen so that the periodogram shows one lump. At a phase of 1.0 rad, the same tones give two lumps of nearly equal height, 0.3 dB apart, at 0.195π and 0.236π. Both sit outside the true pair, so two lumps there are no better than one.

On that record, MUSIC still finds 0.200π and 0.232π.

Worked example

1. Resolution. One bin at $N=64$ is $2\pi/64=0.03125\pi$. The tones are $0.23\pi-0.20\pi=0.03\pi$ apart, so they are $0.03/0.03125=0.96$ of a bin apart. A periodogram would need about $2\pi/(0.03\pi)\approx67$ samples.

2. Count the tones. The eigenvalues are 8.509, 7.426, 0.428 and 0.253, then twelve averaging 0.0101. The fourth is 25.2 times the floor’s average, computed before rounding. Four exponentials make two real tones, and the floor sits near the measured noise power, 0.0090.

3. Read the frequencies. MUSIC’s peaks are at 0.201π and 0.232π, errors of 0.0008π and 0.0023π. That is 0.03 and 0.07 of a bin, where the periodogram gave one lump at 0.221π.

Where you’ll meet this

MUSIC was made for antenna arrays: Schmidt’s paper is titled “Multiple emitter location and signal parameter estimation”. Put $m$ antennas in a row, and a wave arriving at an angle reaches each one a little later than its neighbour. The $m$ antennas’ samples at one instant form a snapshot, and the phase step between neighbours plays the part of $\Omega$. Radar and sonar use this to find the directions of several sources at once.

The same methods estimate frequencies when a record must be short, in communications and in vibration analysis, where a machine or a bridge can have resonances close together. Burg’s method came from geophysics, where Burg presented it in 1967 for short seismic records. Praat, a widely used program for speech analysis, finds the resonances of the voice with it.

When the whole shape of a signal is known, not only its frequencies, a different tool finds it in noise: Matched filters and detection (26.1).

Reference card