Estimation of the Amplitude and Phase of Sinusoids in Noise

 

by Barry Van Veen

In the previous post we developed a method for estimating the amplitude and phase of $K$ sinusoids from $N\ge 2K$ data samples assuming the frequencies are known.  This method applies to sinusoids with arbitrary frequencies.  The approach we used to develop the estimator was to convert the amplitude and phase of each sinusoid into a pair of amplitude parameters that are inserted in a vector ${\bf f}$. We then defined an $N$-by-$2K$ matrix ${\bf H}$ whose columns are  cosine and sine vectors corresponding to the frequencies of interest. This allows the $N$-by-1 vector of data samples ${\bf x}$ to be expressed as

${\bf x}={\bf Hf}$

The matrix ${\bf H}$ is known since the frequencies are known.  The unknown amplitudes and phases are represented by the vector ${\bf f}$.

We showed how to obtain ${\bf f}$ from ${\bf x}$ using the equation

${\bf f}=\left({\bf H}^T{\bf H}\right)^{-1}{\bf H}^T{\bf x}$                                  (1)

All the quantities on the right-hand side of Eqn. (1) are known.  This implies that we can recover the amplitude and phase exactly for arbitrary frequencies.

In practice noise will introduce errors into the estimation.  The effect of noise is the topic of this post.

An Example of Estimating the Amplitude and Phase of Sinusoids

First we illustrate this procedure with an example.  The MATLAB code used to generate the example is available at the end of this post.  We assumed $K=10$ sinusoids and $N=100$ samples.  The frequencies of the sinusoids are chosen as prime numbers 101, 103, 107, 109, 113, 127, 137 149, 157, and 167 Hz. Thus, the only harmonic common to all of them is one Hz. The data was generated assuming integer amplitudes of 3, 2, 1, 4, 1, 3, 2, 1, 4, 1 and phases 0, 30, 45, 60, 90, 0, 30, 45, 60, 90 degrees, respectively. The result of applying Eqn. (1) to this data is shown in Fig. 1.  All the amplitudes and phases are recovered exactly.

Estimation of the Amplitude and Phase of Sinusoids in Noise

There are two effects that make exact estimation impossible in practice: noise and errors in the assumed frequencies.  Here we examine the effect of noise since noise is always present - even if only due to quantization errors when the data is sampled.

Suppose we measure $y[n]=x[n]+w[n]$ where $w[n]$ is noise. In the analysis here we assume the noise is white - that is, uncorrelated from sample to sample and has constant power. This implies that the power of the noise is distributed equally across all frequencies.  We write our noisy data in vector form as ${\bf y}={\bf x+w}$ where ${\bf y}$ and ${\bf w}$ are defined analogously to  ${\bf x}$.  Now apply the estimation procedure in Eqn. (1) to ${\bf y}$ to obtain

$\hat{\bf f}=\left({\bf H}^T{\bf H}\right)^{-1}{\bf H}^T{\bf y}$

$=\left({\bf H}^T{\bf H}\right)^{-1}{\bf H}^T{\bf x} +\left({\bf H}^T{\bf H}\right)^{-1}{\bf H}^T{\bf w}$

$={\bf f}+{\bf e}$

where ${\bf e}=\left({\bf H}^T{\bf H}\right)^{-1}{\bf H}^T{\bf w}$. The hat over ${\bf f}$ explicitly denotes that this is an estimate and no longer the true value.  We see that $\hat{\bf f}$ is a noise corrupted version of the true ${\bf f}$.

The impact of the noise is typically assessed by considering the statistics of the error term ${\bf e}$. If the noise is zero mean, then ${\bf e}$ is also zero mean.  Thus $\hat{\bf f}$ is an unbiased estimate of ${\bf x}$. The covariance of the error is

E$\{{\bf e}{\bf e}^T\}=\left({\bf H}^T{\bf H}\right)^{-1}{\bf H}^T$E$\{{\bf w}{\bf w}^T\}{\bf H}\left({\bf H}^T{\bf H}\right)^{-1}$

Since we are assuming the noise is white, E$\{{\bf w}{\bf w}^T\}=\sigma^2{\bf I}$ where $\sigma^2$ is the noise variance and ${\bf I}$ is the identity matrix. Substituting for the noise covariance gives the estimation error covariance as a relatively simple function of ${\bf H}$

E$\{{\bf e}{\bf e}^T\}=\sigma^2\left({\bf H}^T{\bf H}\right)^{-1}$                                                  (2)

The diagonal elements of E$\{{\bf e}{\bf e}^T\}$ tell us the error variance in the corresponding element of $\hat{\bf f}$ while the off-diagonal elements tell us the correlation between the errors in different elements of $\hat{\bf f}$.

The impact of the noise depends solely on the properties of matrix ${\bf H}$ which in turn depend on the frequencies of interest. It turns out that if the discrete-time frequencies are separated by integer multiples of $2\pi/N$ rads (that is, integer multiples of $f_s/N$ Hz where $f_s$ is the sampling frequency), then the columns of ${\bf H}$ are orthogonal and $\left({\bf H}^T{\bf H}\right)=(N/2){\bf I}$. Thus, E$\{{\bf e}{\bf e}^T\}=\sigma^2(2/N){\bf I}$.  This means that the errors in all components of $\hat{\bf f}$ are uncorrelated and of equal variance, $\sigma^2(2/N)$.  This is the lowest possible error variance. The error variance increases as the frequencies of the sinusoids become closer than $2\pi/N$ rads (or $f_s/N$ Hz). In any case, however, increasing $N$ decreases the error variance as we illustrate in the following example.

An Example of Estimation Error Variance

Here we consider the impact of noise on the estimation of the amplitude and phase of the ten sinusoids we considered in the previous example.  Recall the frequencies are 101, 103, 107, 109, 113, 127, 137 149, 157, and 167 Hz and the sampling frequency is $f_s=1000$ Hz. We consider three different values for $N$: 100, 250, and 500.  Note that for $N=500$, $f_s/N=2$ Hz, which is the minimum spacing between the sinusoid frequencies. Hence, we expect the variance of each term to be equal in this case and to have a value of $\sigma^2/250$.

The variance of the odd-indexed elements of $\hat{\bf f}$ (the $b_k$ in Eqn. (2) of the previous post) are shown in Fig. 2 assuming $\sigma^2=1$.  These represent the variances associated with estimating the amplitudes of the cosine terms.  The variance for the even-indexed elements of $\hat{\bf f}$ (associated with the sine terms) has very similar characteristics and is not shown. Note that the vertical axis is displayed in dB to capture the wide range of variances.  The lower bound on the estimatin error variance is $1/250$ or -24 dB.

Several trends are evident. First, the sinusoids that are closely spaced (101, 103 and 107, 109 Hz) have greater estimation noise variance for $N=100$ and $N=250$.  Second, as $N$ increases, the noise variance decreases significantly.

The key factor that determines whether it is hard (i.e., large error variance) or easy (i.e., small error variance) to accurately estimate the parameters of a sinusoid is how distinct its frequency is relative to nearby sinusoids. If the frequency separation exceeds $f_s/N$ Hz, then it is relatively easy. Adjacent sinusoids that are closer than $f_s/N$ Hz are much more difficult to estimate accurately.  For example, the sinusoids at 107 and 109 Hz have the greatest error variance.  The 2 Hz separation here is significantly less than $f_s/N=10$ Hz when $N=100$. In contrast, all the sinusoids above 120 Hz are separated by 10 Hz or more and have much lower estimation error variance.  When $N=250$ we have $f_s/N=4$ Hz and the sinusoids separated by 2 Hz are estimated much more accurately, although still less than those above 120 Hz.

Is It Possible to Estimate Closely Spaced Sinusoid Amplitudes and Phases?

The estimation error variance for sinusoids spaced by 2 Hz is more than 30 dB greater when $N=100$ compared to $N=250$.  The usefulness of the estimates is not  dependent solely on the error variance, but rather on the amplitude of the quantity of interest relative to the error variance.  For example, if our error variance is 1000, we still may obtain very useful results if the parameter we are trying to estimate has amplitude $10^6$.  On the other hand, if the parameter we are trying to estimate has amplitude 10, and the estimation error variance is 1000, then the estimation error is dominant and the estimate is likely useless.

Thus, it is possible to reliably estimate amplitudes and phases for closely spaced sinusoids provided the amplitudes of the sinusoids are sufficiently large relative to the estimation error variance.  Figure 2 illustrates that as the frequency separation becomes less than $f_s/N$ Hz, the sinusoid amplitude must increase to obtain reliable estimates due to significant increases in the error variance.

This property is intuitively sensible.  If the frequencies are very similar, the two sinusoids are also very similar and separating them relies on relatively subtle differences.  The amplitudes need to be large in order for these subtle differences to stand out above the background noise.

What About the DFT?

The DFT assumes the sinusoids are located at integer multiples of $f_s/N$ Hz.  This implies that estimation of amplitude and phase parameters in our example requires $f_s/N=1$ Hz or $N=1000$ as all the frequencies of interest are odd integers.

In contrast, the procedure described in this and the previous post applies for arbitrary $N$ and arbitrary frequency spacings.

MATLAB Script

m-file: Amplitude and Phase of Sinusoids in Noise