Finding Frequencies Without an FFT: Sinusoids as Autoregressive Models
Suppose you have three consecutive samples \(y_0, y_1, y_2\) of a cosine, and you know nothing else: not its frequency, amplitude or phase. With no noise, you can still recover the frequency exactly:
\[ \cos\omega=\frac{y_0+y_2}{2y_1}. \]No FFT, no long recording. (We just need \(y_1\neq 0\).) This works because every cosine obeys the same simple recurrence,
\[ y_t=2\cos(\omega)\,y_{t-1}-y_{t-2}, \qquad t\geq 2, \]and setting \(t=2\) gives the formula above. The cosine satisfies an AR(2) relation: each sample is a fixed linear combination of the previous two, with a coefficient that depends only on the frequency.
The same is true of any sum of sinusoids, with one more past sample for each component:
A sum of \(n\) complex sinusoids is AR(\(n\))
For constant amplitudes \(A_i\) and frequencies \(\omega_i\), let
\[ y_t=\sum_{i=1}^{n}A_i e^{j\omega_i t}, \qquad t=0,1,2,\ldots \]Then every sample is a fixed linear combination of the previous \(n\):
\[ y_t=\sum_{k=1}^{n}b_k\,y_{t-k},\qquad t\geq n, \]where the coefficients \(b_k\) come from expanding the polynomial
\[ P(z)=\prod_{i=1}^{n}\bigl(z-e^{j\omega_i}\bigr)=z^n-\sum_{k=1}^{n}b_k z^{n-k}. \]In particular, the roots of \(P(z)\) are \(e^{j\omega_i}\): their angles are the frequencies.
This result is useful because we can run it backwards. Given samples of an unknown signal, fit the coefficients \(b_k\), find the roots of \(P(z)\), and read off the frequencies. This is Prony's method, published by Gaspard de Prony in 1795. An FFT can only separate frequencies about \(f_s/N\) apart, where \(N\) is the number of samples; on clean data, Prony's method recovers \(n\) frequencies exactly from \(2n\) samples. The same idea underlies linear prediction in speech coding and modern high-resolution frequency estimators.
I first worked this out from scratch in 2014, and only afterwards found out that Prony had got there about 220 years earlier. This post is a cleaned-up version of those old notes.
We will build up the proof one sinusoid at a time: one, then two, then \(n\). Throughout, we assume there is no noise: an idealisation that lets us see the relationship exactly.
One Sinusoid
Consider one complex sinusoid:
\[ y_t=A_1e^{j\omega_1 t}. \]Here \(j^2=-1\), and \(\omega_1\) is the angular frequency in radians per sample. The amplitude \(A_1=|A_1|e^{j\phi_1}\) is a complex constant, so its magnitude sets the size of the oscillation and its angle sets the starting phase. For example, the real part of \(y_t\) is \(|A_1|\cos(\omega_1 t+\phi_1)\).
The key fact is that delaying a complex sinusoid just multiplies it by a constant. For \(t\geq 1\),
\[ y_{t-1}=A_1e^{j\omega_1(t-1)}=e^{-j\omega_1}A_1e^{j\omega_1t}=e^{-j\omega_1}y_t. \]Let \(a_1=e^{j\omega_1}\). Multiplying both sides by \(a_1\) gives
\[ y_t=a_1y_{t-1}. \tag{1}\label{eq:ar1} \]This is an autoregressive relation of order one, or AR(1): the current sample is a fixed multiple of the previous one. Equivalently, the combination \(y_t-a_1y_{t-1}\) is zero. Subtracting \(a_1\) times the previous sample cancels the sinusoid.
Two Sinusoids
Now add a second sinusoid:
\[ y_t=A_1e^{j\omega_1t}+A_2e^{j\omega_2t}. \]Apply the same subtraction as before, with \(a_1=e^{j\omega_1}\) and \(a_2=e^{j\omega_2}\). Call the result \(v_t\):
\[ \begin{aligned} v_t&=y_t-a_1y_{t-1}\\ &=\underbrace{\bigl(A_1-a_1e^{-j\omega_1}A_1\bigr)}_{=\,0}e^{j\omega_1t} +\bigl(A_2-a_1e^{-j\omega_2}A_2\bigr)e^{j\omega_2t}\\ &=\bar A_2\,e^{j\omega_2t}, \qquad \bar A_2=\bigl(1-a_1/a_2\bigr)A_2. \end{aligned} \]The first sinusoid has cancelled. What remains is one sinusoid at \(\omega_2\), with a new constant amplitude \(\bar A_2\). So the one-sinusoid argument applies to \(v_t\), and subtracting \(a_2\) times its previous sample cancels it too:
\[ v_t-a_2v_{t-1}=0, \qquad t\geq 2. \]Substitute \(v_t=y_t-a_1y_{t-1}\) to write this in terms of the original signal:
\[ \begin{aligned} 0&=(y_t-a_1y_{t-1})-a_2(y_{t-1}-a_1y_{t-2})\\ &=y_t-(a_1+a_2)y_{t-1}+a_1a_2y_{t-2}. \end{aligned} \]Rearranging,
\[ y_t=(a_1+a_2)y_{t-1}-a_1a_2y_{t-2}, \qquad t\geq 2. \tag{2}\label{eq:ar2} \]Two complex sinusoids give an AR(2) relation. This matches the result at the top: \((z-a_1)(z-a_2)=z^2-(a_1+a_2)z+a_1a_2\), so \(b_1=a_1+a_2\) and \(b_2=-a_1a_2\).
Check: one real cosine
Now we can explain the cosine example from the start. Write the cosine as two complex sinusoids:
\[ \cos(\omega t)=\tfrac12 e^{j\omega t}+\tfrac12 e^{-j\omega t}. \]This is the two-sinusoid case with \(\omega_1=\omega\) and \(\omega_2=-\omega\). Therefore
\[ a_1+a_2=e^{j\omega}+e^{-j\omega}=2\cos\omega, \qquad a_1a_2=e^{j\omega}e^{-j\omega}=1, \]and \(\eqref{eq:ar2}\) becomes \(y_t=2\cos(\omega)\,y_{t-1}-y_{t-2}\). The amplitudes never entered, so the same recurrence holds for any \(B\cos(\omega t+\phi)\).
Any Number of Sinusoids
With \(n\) sinusoids, we do exactly what we did for two: cancel one sinusoid, notice that we are left with the same problem with one fewer, and repeat.
Cancel one sinusoid
To keep track of how many sinusoids are left, write the signal with a superscript:
\[ y_t^{[n]}=\sum_{i=1}^{n}A_i e^{j\omega_i t}, \qquad a_i=e^{j\omega_i}. \]The superscript \([n]\) says that \(n\) sinusoids remain. As before, delay by one sample, multiply by \(a_1\), and subtract. Each term is delayed separately, and delaying term \(i\) multiplies it by \(e^{-j\omega_i}=1/a_i\):
\[ \begin{aligned} y_t^{[n-1]}&=y_t^{[n]}-a_1y_{t-1}^{[n]}\\ &=\sum_{i=1}^{n}\Bigl(1-\frac{a_1}{a_i}\Bigr)A_i e^{j\omega_i t}\\ &=\sum_{i=2}^{n}\bar A_i\, e^{j\omega_i t}, \qquad \bar A_i=\Bigl(1-\frac{a_1}{a_i}\Bigr)A_i. \end{aligned} \]For \(i=1\), the factor is \(1-a_1/a_1=0\), so the first sinusoid disappears. The other \(n-1\) sinusoids keep their frequencies; only their amplitudes change. In the two-sinusoid section, this was \(v_t\).
Repeat until nothing is left
Notice that \(y_t^{[n-1]}\) is the same kind of signal we started with, just with one fewer sinusoid. So we can apply the same step to it, this time with \(a_2\) to cancel the sinusoid at \(\omega_2\). In general, the \(r\)th step cancels the sinusoid at \(\omega_r\):
\[ y_t^{[n-r]}=y_t^{[n-r+1]}-a_r\,y_{t-1}^{[n-r+1]}, \qquad r=1,\ldots,n. \tag{3}\label{eq:cancellation} \]Each step uses one more past sample, so after \(r\) steps the result is valid for \(t\geq r\). After all \(n\) steps, no sinusoids are left:
\[ y_t^{[0]}=0, \qquad t\geq n. \]Unwind back to the original signal
We now have a statement about \(y_t^{[0]}\), but we want one about the original signal \(y_t^{[n]}\). To unwind the steps, it helps to write the delay as an operator. Let \(z^{-1}\) shift a sequence back by one sample:
\[ z^{-1}y_t=y_{t-1}, \qquad z^{-2}y_t=y_{t-2}. \]Each step in \(\eqref{eq:cancellation}\) then becomes a single factor:
\[ y_t^{[n-r]}=(1-a_rz^{-1})\,y_t^{[n-r+1]}. \]Start from \(y_t^{[0]}=0\) and substitute one step at a time:
\[ \begin{aligned} y_t^{[0]}&=0\\ (1-a_nz^{-1})\,y_t^{[1]}&=0\\ (1-a_nz^{-1})(1-a_{n-1}z^{-1})\,y_t^{[2]}&=0\\ &\;\;\vdots\\ \prod_{i=1}^{n}(1-a_iz^{-1})\,y_t^{[n]}&=0. \end{aligned} \]The last line is about the original signal. Each factor \(1-a_iz^{-1}\) removes the sinusoid at \(\omega_i\), and the order of the factors does not matter: we could have cancelled the sinusoids in any order.
Expand the product
To read off the recurrence, multiply out the factors. For two sinusoids,
\[ (1-a_1z^{-1})(1-a_2z^{-1})=1-(a_1+a_2)z^{-1}+a_1a_2z^{-2}, \]which is exactly \(\eqref{eq:ar2}\): the coefficient of \(z^{-k}\) multiplies \(y_{t-k}\). With all \(n\) factors,
\[ \begin{aligned} \prod_{i=1}^{n}(1-a_iz^{-1}) &=1-\Bigl(\sum_{i}a_i\Bigr)z^{-1} +\Bigl(\sum_{i<k}a_ia_k\Bigr)z^{-2}-\cdots +(-1)^n\Bigl(\prod_{i}a_i\Bigr)z^{-n}\\ &=1-\sum_{k=1}^{n}b_k z^{-k}. \end{aligned} \]The \(b_k\) are just the coefficients of this expansion: \(b_1=\sum_i a_i\), \(b_2=-\sum_{i<k}a_ia_k\), and so on, down to \(b_n=(-1)^{n+1}\prod_i a_i\). Writing \(y_t\) for \(y_t^{[n]}\) again and turning the delays back into samples,
\[ \Bigl(1-\sum_{k=1}^{n}b_k z^{-k}\Bigr)y_t=y_t-\sum_{k=1}^{n}b_k y_{t-k}=0. \]Therefore,
\[ y_t=\sum_{k=1}^{n}b_k y_{t-k}, \qquad t\geq n. \tag{4}\label{eq:arn} \]This is the AR(\(n\)) relation from the top of the post.
The roots are the frequencies
Finally, treat the product as an ordinary polynomial and multiply by \(z^n\) to clear the negative powers. Each factor \(1-a_iz^{-1}\) becomes \(z-a_i\):
\[ P(z)=z^n\prod_{i=1}^{n}(1-a_iz^{-1})=\prod_{i=1}^{n}(z-a_i)=z^n-\sum_{k=1}^{n}b_k z^{n-k}. \tag{5}\label{eq:polynomial} \]Its roots are \(a_i=e^{j\omega_i}\), so the angle of each root is one of the frequencies. One sinusoid gives AR(1), two give AR(2), and \(n\) give AR(\(n\)), and in every case the frequencies are hidden in the roots.
Estimating the Frequencies
In practice, we do not know the frequencies or the coefficients. We start with samples \(y_0,\ldots,y_{N-1}\), choose the number of components \(n\), and write \(\eqref{eq:arn}\) once for each \(t=n,\ldots,N-1\):
\[ \begin{bmatrix} y_{n-1}&y_{n-2}&\cdots&y_0\\ y_n&y_{n-1}&\cdots&y_1\\ \vdots&\vdots& &\vdots\\ y_{N-2}&y_{N-3}&\cdots&y_{N-n-1} \end{bmatrix} \begin{bmatrix}b_1\\b_2\\\vdots\\b_n\end{bmatrix} =\begin{bmatrix}y_n\\y_{n+1}\\\vdots\\y_{N-1}\end{bmatrix}. \]This is linear in the unknown \(b_k\). For exact data with distinct frequencies and nonzero amplitudes, \(2n\) samples give enough equations; with more, we solve by least squares. Then:
- Find the roots \(\hat a_i\) of \(P(z)\) in \(\eqref{eq:polynomial}\).
- Read off the frequencies: \[ \hat\omega_i=\arg(\hat a_i), \qquad \hat f_i=\frac{f_s}{2\pi}\hat\omega_i, \] where \(f_s\) is the sampling frequency. Frequencies are only identified modulo \(2\pi\) radians per sample (or modulo \(f_s\) in hertz).
- With the frequencies known, \(y_t=\sum_i A_i e^{j\hat\omega_i t}\) is linear in the amplitudes. Solve it for \(A_i\); their magnitudes and angles give the amplitudes and phases.
Practical Considerations
Real signals such as audio tones or mechanical vibrations are sums of real sines and cosines. Each real component is two complex sinusoids, as in the cosine check, so \(m\) real components need an AR(\(2m\)) model, with real coefficients.
The derivation also assumes an exact sum of sinusoids and a known number of components. Noise makes the recurrence approximate, and the estimated roots can move off the unit circle. Closely spaced frequencies make the problem ill-conditioned, so more samples do not automatically remove the sensitivity to noise. Choosing \(n\) is itself a modelling decision.