Ruiming Guo, Thierry Blu
This article introduces a simple formula that provides the exact frequency of a pure sinusoid from just two samples of its discrete-time Fourier transform (DTFT). Even when the signal is not a pure sinusoid, this formula still works in a very good approximation (optimally after a single refinement), paving the way for the high-resolution frequency tracking of quickly varying signals or simply improving the frequency resolution of the peaks of a discrete Fourier transform (DFT).
We learn (and teach!) in college that the frequency content of a signal is encrypted in its FT and that the main frequency modes of a sampled signal are the peaks of its DFT. We also learn that the accuracy of these frequency values is limited by the inverse of the time range of the signal (Heisenberg uncertainty), which correlates nicely with the inherent frequency resolution of a DFT: ${2}{\pi}{/}{N}$, if N is the number of samples of the signal. This knowledge is so deeply rooted that it is hard to reconcile with the fact that the frequencies of a signal made of a sum of K complex exponentials can be recovered exactly from as few as 2K samples by using a two-century-old method due to Gaspard de Prony. This apparent contradiction is resolved by recognizing that Heisenberg uncertainty relies on a much weaker signal assumption (basically, that its time and frequency uncertainties are finite) than the sum-of-exponentials model. And it is only when this model is inexact that the estimated frequencies may be inaccurate, with their uncertainty now given by the Cramér-Rao lower bound, which assumes unbiased estimators and known noise statistics.
The contrast between the Fourier approach (analytic, intuitive, but Heisenberg limited) and Prony’s method (algebraic, black box, but exact) has made it difficult to envision a higher frequency resolution that would rely on the DFT. Yet, given that the DFT coefficients are just samples of a very smooth function (the DTFT), analytic considerations suggest that this function can be approximated locally by a quadratic polynomial, leading to an estimate of its off-grid peak location; as few as three samples around the max of the DFT already provide a very good estimate of this frequency [1], [2].
The main motivation for this article is to put forward a formula that provides the frequency value of the maximum of the DTFT from just two DTFT coefficients of a single-frequency signal; see the detailed setting in “Box 1: Notations.” Not only is this formula exact, but it is also very robust to the inaccuracies of the model, as we shall see later. More precisely, if we assume that the unknown frequency ${\omega}_{0}$ of the signal ${x}_{n}$ lies inside an “uncertainty band” $\left[{{\omega}_{1},{\omega}_{2}}\right]$ where ${\omega}_{2}{-}{\omega}_{1}$ is an integer multiple of ${2}{\pi}{/}{N}$, this formula specifies how the amplitudes of the DTFT of ${x}_{n}$ at ${\omega}_{1}$ and ${\omega}_{2}$ should be combined so as to recover ${\omega}_{0}$(see the visualization in Figure 1) \begin{align*} & {\omega}_{0} = \frac{{\omega}_{1} + {\omega}_{2}}{2} + {2}\arctan \\ & \left({\tan\left({\frac{{\omega}_{2}{-}{\omega}_{1}}{4}}\right)\frac{\left|{{X}{(}{\omega}_{2}{)}}\right|{-}\left|{{X}{(}{\omega}_{1}{)}}\right|}{\left|{{X}{(}{\omega}_{2}{)}}\right| + \left|{{X}{(}{\omega}_{1}{)}}\right|}}\right){.} \tag{1} \end{align*}
Figure 1. The frequency ${\omega}_{0}$ (blue) of a single-frequency signal is obtained from the DTFT values at the endpoints (red) of a frequency interval known to contain ${\omega}_{0}$ (assumption: ${\omega}_{2}{-}{\omega}_{1} = {\text{integer}}\,{\times}\,{2}{\pi}{/}{N}{)}$. Dotted line: a full period of the DTFT ${X}{(}{\omega}{)}$ of the signal.
Box 1: Notations
A proof is provided in “Box 2: Step-by-Step Didactic Proof of (1),” requiring only elementary electrical and electronics engineering math knowledge. This formula becomes even simpler when the uncertainty bandwidth is small (i.e., ${\omega}_{2}{-}{\omega}_{1}\ll{\pi}{)}$ \begin{align*} {\omega}_{0}\approx{\omega}_{1} & {\times}\frac{\left|{{X}{(}{\omega}_{1}{)}}\right|}{\left|{{X}{(}{\omega}_{1}{)}}\right| + \left|{{X}{(}{\omega}_{2}{)}}\right|} \\ & + {\omega}_{2}\times\frac{\left|{{X}{(}{\omega}_{2}{)}}\right|}{\left|{{X}{(}{\omega}_{1}{)}}\right| + \left|{{X}{(}{\omega}_{2}{)}}\right|}\end{align*} which is what intuition would suggest: a weighting of the end frequencies based on the relative magnitude of their DTFT.
Box 2: Step-by-Step Didactic Proof of (1)
Steps
In practice, this formula is most useful when the uncertainty band is smallest, i.e., ${\omega}_{2}{-}{\omega}_{1} = {2}{\pi}{/}{N}$. Then, a straightforward procedure for estimating a single frequency from N uniform signal samples is to
This works because, for a single-frequency signal, the peak of the DTFT is always within $\pm{\pi}{/}{N}$ of the maximum of the DFT. But it would also be possible to bypass the computation of the full DFT if a rough estimate of ${\omega}_{0}$ were available, e.g., when the frequency of the signal is continuously changing (tracking between successive signal windows, as in radar applications) or when the frequency is a priori known up to some perturbation (physical resonance experiments, laser-based optical measurements, etc.).
When the single-exponential model is not exact, (1) is still a very robust estimator of its frequency. This is particularly so when the uncertainty bandwidth ${\omega}_{2}{-}{\omega}_{1}$ is reduced to ${2}{\pi}{/}{N}$ an assumption that we will make from now on. Indeed, consider a noise model ${x}_{n} = {a}_{0}{\text{e}}^{jn{\omega}_{0}} + {b}_{n}$ where the complex-valued samples ${b}_{n}$ are independent realizations of a Gaussian random variable with variance ${\sigma}^{2}$—the “additive white Gaussian noise” assumption. When the number of samples N is large enough, a linearization of (1) makes it possible to calculate (tediously; not shown here) the standard deviation $\Delta{\omega}$ of the frequency estimation error. In particular, in the case where ${\omega}_{0}$ is at the center of the interval $\left[{{\omega}_{1},{\omega}_{2}}\right]$, this error behaves according to \[\Delta{\omega} = \frac{\frac{\sqrt{2}{\pi}^{2}}{4}}{N\sqrt{N}\text{SNR}} = \frac{{\pi}^{2}}{4\sqrt{6}}\Delta{\omega}_{\text{CR}} \tag{2} \] where ${\text{SNR}} = \mid{a}_{0}\mid{/}{\sigma}$ and where $\Delta{\omega}_{\text{CR}}$ is the Cramér-Rao lower bound of the problem (see [3] for a calculation). Obviously, $\Delta{\omega}$ is very close to $\Delta{\omega}_{\text{CR}}$, within less than 1%. When ${\omega}_{0}$ is closer to the extremities of the interval $\left[{{\omega}_{1},{\omega}_{2}}\right]$, $\Delta{\omega}$ deviates from the Cramér-Rao lower bound by up to 80%, still a very low error in absolute terms. In fact, a simple refinement of the trick, as depicted in Figure 2, shows how to attain this lower bound, outlining the near optimality of this procedure; no other unbiased single-frequency estimation algorithm would be able to improve this performance by more than 1%.
Figure 2. A refined trick: a double application of (1) achieves a quality that is equivalent to the best unbiased single-frequency estimation algorithm.
This is confirmed in Figure 3 by simulations that consist of 1 million tests where
Figure 3. Histograms from 1 million random tests (number of samples $N$, SNR, frequency, and noise). The error that results from using Jacobsen’s frequency estimator [1], the trick (1), or the refined trick (Figure 2) is further normalized by the Cramér-Rao lower bound of the estimation problem ${(}\Delta{\omega}_{\text{CR}} = {2}\sqrt{3}{N}{}^{{-}{3}{/}{2}}{\text{SNR}}^{{-}{1}}{)}$. The standard deviations of the three estimators are 1.5325, 1.3008, and 1.0092, respectively. PDF: probability density function.
In each test, the uncertainty band before refinement is set by ${\omega}_{1} = {\omega}_{\text{DFT}}{-}{\pi}{/}{N}$ and ${\omega}_{2} = {\omega}_{\text{DFT}} + {\pi}{/}{N}$. For comparison purposes, Figure 3 also shows the distribution of errors of Jacobsen’s estimator [1], [2], which is based on three consecutive DFT coefficients around ${\omega}_{\text{DFT}}$. The better performance of our formula is likely due to the higher SNR enjoyed by the two DFT coefficients around the maximum of the DTFT in comparison to the three DFT coefficients used by Jacobsen’s formula, one of which has a significantly lower SNR because it is further away from the DTFT peak by more than ${2}{\pi}{/}{N}$. A somewhat milder difference is also that (1) assumes an exact single-frequency model, whereas Jacobsen’s formula is but a local quadratic approximation.
Beyond just noise, when the single-frequency model is rendered inaccurate due to, e.g., quantization, sample windowing, or the addition of other sinusoidal/polynomial terms, the frequency estimation error of (1) is controlled by $\varepsilon$, the maximum error of the magnitude of the DTFT at the frequencies ${\omega}_{1}$ and ${\omega}_{2}$, according to \begin{align*}\mid{\bar{\omega}}_{0}{-}{\omega}_{0}\mid & \leq\underbrace{\frac{{4}\varepsilon\tan\left({\frac{\pi}{2N}}\right)}{{N}\mid{a}_{0}\mid}} \\ & \quad\approx\frac{{2}{\pi}\varepsilon}{{N}{}^{2}\mid{a}_{0}\mid}\,{\text{for}}\,{\text{large}}\,{N} \tag{3} \end{align*} where $\mid{a}_{0}\mid$ is the amplitude of the single-frequency “ground-truth” signal. (See the proof in the “Error Bound (3)” section in “Other Proofs.”)
Note that because it is valid for every single “noise” instance, this bound is of a very different nature than the statistical result (2)—an average over infinitely many additive white Gaussian noise realizations. Interestingly, the inaccuracies of the model outside the uncertainty band do not contribute to the estimation error, which suggests that (3) can be used to predict the accuracy of a multiple-frequency estimation problem that uses the single-frequency trick.
This formula can easily be used in a multiple-frequency scenario provided that the frequencies to estimate are sufficiently separated. A straightforward procedure consists of first locating the isolated peaks of the DFT of the signal (e.g., using MATLAB’s findpeaks function) and then applying (1) to refine each frequency individually; see an example in Figure 4. The estimation error of each frequency can be quantified by using the bound (3) where, in the absence of other noise, the data inaccuracy $\varepsilon$ in the neighborhood of that frequency is essentially caused by the tail of the DTFT of the other frequencies. An example of such a calculation is shown in the “Error Bound (4)” section in “Other Proofs,” leading to the following statement. Assume that the frequencies ${\omega}_{k}$ of the signal are distant from each other (modulo ${(}{-}{\pi},{\pi}{]}{)}$ by at least ${\delta}{\omega}{>}{\pi}{/}{N}$ and that the amplitude of the dominant sinusoid is A; then, the estimation error of any of the frequencies of the signal is bounded according to \begin{align*} & \left|{{\bar{\omega}}_{k}{-}{\omega}_{k}}\right|\leq\mathop{\underbrace{\frac{{2}{\pi}}{N}}}\limits_{{\text{DFT}}\,{\text{resolution}}} \\ & \quad\times\mathop{\underbrace{\frac{{2}{(}{K}{-}{1}{)}\tan{(}{\pi}{/}{(}{2}{N}{)}{)}}{{\pi}\sin{((}{\delta}{\omega}{-}{\pi}{/}{N}{)}{/}{2}{)}}\frac{A}{\mid{a}_{k}\mid}}}\limits_{\mbox{“}{\text{super}}{-}{\text{resolution}}\mbox{”}\,{\text{coefficient}}}{.} \tag{4} \end{align*}
Figure 4. Multiple frequencies (three complex exponentials, 12 samples) can be accurately estimated by first locating the isolated peaks (e.g., using findpeaks in MATLAB) and then applying (1) individually. The actual estimation errors of the three frequencies (from left to right) are roughly ${(}{0}{.}{01},{0}{.}{02},{0}{.}{03}{)}\,{\times}\,{2}{\pi}{/}{12}$, well below the resolution ${2}{\pi}{/}{12}$ of the DFT; for comparison, the upper bound (4) provides the much more conservative values ${(}{0}{.}{66},{0}{.}{45},{0}{.}{58}{)}\,{\times}\,{2}{\pi}{/}{12}$.
Here, ${\bar{\omega}}_{k}$, ${\omega}_{k}$, $\mid{\text{and}}\,{a}_{k}\mid$ are the estimated frequency, the ground-truth frequency, and its amplitude, respectively. Despite its coarseness (see Figure 4), this inequality already demonstrates superresolution potential since the “superresolution” coefficient is usually smaller than one and, in fact, tends to zero when N tends to infinity (for fixed ${\delta}{\omega}$).
Other Proofs
Error Bound (3)
A direct proof of this inequality uses the fact that $\mid\arctan{a}{-}\arctan{b}\mid\leq\mid{a}{-}{b}\mid$ and the triangle inequality $\mid{a} + {b}\mid\leq\mid{a}\mid + \mid{b}\mid$. More specifically, denoting by ${X}_{1}$, ${X}_{2}$ the discrete-time Fourier transform (DTFT) of the “ground-truth” signal at ${\omega}_{1}$, ${\omega}_{2}$, and by ${\varepsilon}_{1}$, ${\varepsilon}_{2}$ the errors (caused by noise or otherwise) on $\mid{X}_{1}\mid$, $\mid{X}_{2}\mid$, we have \begin{align*}\left|{{\bar{\omega}}_{0}{-}{\omega}_{0}}\right| & = \left|{{2}\arctan\left({\tan\left({\frac{\pi}{2N}}\right)\frac{\left|{{X}_{1}}\right| + {\varepsilon}_{1}{-}\left|{{X}_{2}}\right|{-}{\varepsilon}_{2}}{\left|{{X}_{1}}\right| + {\varepsilon}_{1} + \left|{{X}_{2}}\right| + {\varepsilon}_{2}}}\right)}\right. \\ & \qquad\left.{{-}{2}\arctan\left({\tan(\frac{\pi}{2N})\frac{\left|{{X}_{1}}\right|{-}\left|{{X}_{2}}\right|}{\left|{{X}_{1}}\right| + \left|{{X}_{2}}\right|}}\right)}\right| \\ & \leq{2}\tan\left({\frac{\pi}{2N}}\right)\left|{\underbrace{\frac{\left|{{X}_{1}}\right| + {\varepsilon}_{1}{-}\left|{{X}_{2}}\right|{-}{\varepsilon}_{2}}{\left|{{X}_{1}}\right| + {\varepsilon}_{1} + \left|{{X}_{2}}\right| + {\varepsilon}_{2}}{-}\frac{\left|{{X}_{1}}\right|{-}\left|{{X}_{2}}\right|}{\left|{{X}_{1}}\right| + \left|{{X}_{2}}\right|}}}\right| \\ & \qquad\qquad\qquad = \frac{{2}{\varepsilon}_{1}{(}\left|{{X}_{2}}\right| + {\varepsilon}_{2}{)}{-}{2}{\varepsilon}_{2}{(}\left|{{X}_{1}}\right| + {\varepsilon}_{1}{)}}{{(}\left|{{X}_{1}}\right| + \left|{{X}_{2}}\right|{)(}\left|{{X}_{1}}\right| + {\varepsilon}_{1} + \left|{{X}_{2}}\right| + {\varepsilon}_{2}{)}} \\ & \leq\underbrace{{2}\tan\left({\frac{\pi}{2N}}\right)\frac{\max{(}\mid{2}{\varepsilon}_{1}\mid,\mid{2}{\varepsilon}_{2}\mid{)(}\left|{{X}_{1}}\right| + {\varepsilon}_{1} + \left|{{X}_{2}}\right| + {\varepsilon}_{2}{)}}{{(}\left|{{X}_{1}}\right| + \left|{{X}_{2}}\right|{)(}\left|{{X}_{1}}\right| + {\varepsilon}_{1} + \left|{{X}_{2}}\right| + {\varepsilon}_{2}{)}}} \\ & \qquad\qquad\qquad = {4}\tan\left({\frac{\pi}{2N}}\right)\frac{\max\left({\mid{\varepsilon}_{1}\mid,\mid{\varepsilon}_{2}\mid}\right)}{\left|{{X}_{1}}\right| + \left|{{X}_{2}}\right|}\end{align*} which leads to the inequality (3) after noticing that $\mid{X}_{1}\mid + \mid{X}_{2}\mid\geq{N}\mid{a}_{0}\mid$ (because ${\omega}_{0}\,{\in}\,{[}{\omega}_{1},{\omega}_{2}{]}$).
Error Bound (4)
Denoting by ${\omega}_{1},\,{\omega}_{2},\ldots{\omega}_{K}$ the K different frequencies and ${a}_{1},\,{a}_{2},\ldots{a}_{K}$ the associated (complex-valued) amplitudes, the DTFT of the samples ${x}_{n}$ is given by \[{X}{(}{\omega}{)} = \mathop{\sum}\limits_{{k} = {1}}^{K}{{a}_{k}}\frac{{e}^{{jN}{(}{\omega}_{k}{-}{\omega}{)}}{-}{1}}{{e}^{{j}{(}{\omega}_{k}{-}{\omega}{)}}{-}{1}}{.}\]
Evaluating the estimation error of the frequency ${\omega}_{{k}_{0}}$ using (1) requires calculating the bound $\varepsilon$ in (3), i.e., the maximum error between ${X}{(}{\omega}{)}$ and the DTFT of a single-frequency model, when $\mid{\omega}{-}{\omega}_{{k}_{0}}\mid\leq{\pi}{/}{N}$ (with the hypothesis that the minimum distance between ${\omega}_{{k}_{0}}$ and the other ${\omega}_{k}$ is at least ${\delta}{\omega}{>}{\pi}{/}{N}$) \begin{align*}&\left|{X}{(}{\omega}{)}{-}{a}_{{k}_{0}}\frac{{e}^{{jN}{(}{\omega}_{{k}_{0}}{-}{\omega}{)}}{-}{1}}{{e}^{{j}{(}{\omega}_{{k}_{0}}{-}{\omega}{)}}{-}{1}}\right| = \left|\mathop{\sum}\limits_{{k}\ne{k}_{0}}{{a}_{k}}\frac{{e}^{{jN}{(}{\omega}_{k}{-}{\omega}{)}}{-}{1}}{{e}^{{j}{(}{\omega}_{k}{-}{\omega}{)}}{-}{1}}\right| \\ & \quad\leq\mathop{\sum}\limits_{{k}\ne{k}_{0}}{\mid{a}_{k}\mid\left|{\frac{\sin{(}{N}{(}{\omega}_{k}{-}{\omega}{)}{/}{2}{)}}{\sin{((}{\omega}_{k}{-}{\omega}{)}{/}{2}{)}}}\right|}\,\,{(}{\text{using triangle inequality and Euler}}{’}{\text{s formula}}{)} \\ & \quad\leq\mathop{\sum}\limits_{{k}\ne{k}_{0}}\frac{A}{\left|{\sin}{((}{\omega}_{k}{-}{\omega}{)}{/}{2}{)}\right|}\quad\quad\quad\quad{(}{\text{denoting}}\,{A} = \mathop{\max}\limits_{{k} = {1}\ldots{K}}\left|{{a}_{k}}\right|{)}\\ & \quad\leq\frac{{(}{K}{-}{1}{)}{A}}{\mathop{\min}\limits_{{k}\ne{k}_{0}}\left|{\sin{((}{\omega}_{k}{-}{\omega}{)}{/}{2}{)}}\right|} \\ & \quad\leq\frac{{(}{K}{-}{1}{)}{A}}{\sin{((}{\delta}{\omega}{-}{\pi}{/}{N}{)}{/}{2}{)}}\quad\quad\quad\quad\,\,\,{(}{\text{since}}\,\left|{{\omega}{-}{\omega}_{{k}_{0}}}\right|\,leq{\pi}{/}{N}\lt{\delta}{\omega}{).} \end{align*}
The right-hand side of the last inequality provides an upper bound for $\varepsilon$, which we can use in (3) to find \[\left|{{\bar{\omega}}_{{k}_{0}}{-}{\omega}_{{k}_{0}}}\right|\leq\frac{{2}{\pi}}{N}\times\frac{{2}{(}{K}{-}{1}{)}\tan{(}{\pi}{/}{(}{2}{N}{))}}{{\pi}\sin{((}{\delta}{\omega}{-}{\pi}{/}{N}{)}{/}{2}{)}}\frac{A}{\mid{a}_{{k}_{0}}\mid}{.}\]
Empirically, a minimum value of ${4}{\pi}{/}{N}$ for ${\delta}{\omega}$, or two DFT bins, seems to be sufficient to obtain good frequency estimates. Of course, this cheap approach to high-resolution multifrequency estimation is not optimal, yet it could be used as the starting point of any iterative algorithm designed to maximize the likelihood of the problem.
The frequency of a single complex exponential can be found exactly using the magnitude of only two samples of its DTFT, as this article shows. In the presence of noise or other inaccuracies, the trick that we provide is very robust and can even be iterated once to reach the theoretical optimum (Cramér-Rao lower bound)—up to less than 1%. The robustness of this formula makes it possible to, e.g., refine the peaks of the DFT of a signal, but we also anticipate that it can be used as a tool for high-resolution frequency estimation. For teaching purposes, we provide a step-by-step proof that requires only undergraduate signal processing knowledge.
Thierry Blu is the corresponding author.
Ruiming Guo (ruiming.guo@imperial.ac.uk) received his Ph.D. degree in electronic engineering at the Chinese University of Hong Kong (CUHK), Hong Kong Special Administrative Region, China, in 2021. He is currently a postdoctoral research associate at the Department of Electrical and Electronic Engineering, Imperial College London, SW7 2AZ London, U.K., working with Prof. Ayush Bhandari on computational imaging and modulo sampling. He worked as a postdoctoral research fellow with Prof. Thierry Blu at the Electronic Engineering (EE) Department at CUHK, from 2021 to 2022. He received the Postgraduate Student Research Excellence Award from the EE Department of CUHK in 2022. His research interests include sparse signal processing, sampling theory, inverse problems, modulo sampling, and computational sensing and imaging.
Thierry Blu (thierry.blu@m4x.org) received his Ph.D. degree from Télécom Paris (ENST) in 1996. He is a professor in the Department of Electronic Engineering, the Chinese University of Hong Kong, Sha Tin New Territories, Hong Kong Special Administrative Region, China, where he has been since 2008. He received two Best Paper Awards (2003 and 2006) and is the coauthor of a paper that received a Young Author Best Paper Award (2009), all from the IEEE Signal Processing Society. His research interests include wavelets, approximation and sampling theory, sparse representations, biomedical imaging, optics, and wave propagation. He is a Fellow of IEEE.
[1] E. Jacobsen and P. Kootsookos, “Fast, accurate frequency estimators [DSP Tips & Tricks] ,” IEEE Signal Process. Mag., vol. 24, no. 3, pp. 123–125, May 2007, doi: 10.1109/MSP.2007.361611.
[2] Ç. Candan, “A method for fine resolution frequency estimation from three DFT samples,” IEEE Signal Process. Lett., vol. 18, no. 6, pp. 351–354, Jun. 2011, doi: 10.1109/LSP.2011.2136378.
[3] P. Stoica and A. Nehorai, “MUSIC, maximum likelihood, and Cramér-Rao bound,” IEEE Trans. Acoust., Speech, Signal Process., vol. 37, no. 5, pp. 720–741, May 1989, doi: 10.1109/29.17564.
Digital Object Identifier 10.1109/MSP.2023.3311592