Finding the Amplitude of a Sinusoid in Noise

I'm trying to solve for the amplitude and frequency of a sinusoid embedded in zero-mean gaussian white noise.

I am supplied a sample file of a 40000 element array. I first took the autocorrelation function in matlab and then took the fft of this to find the power spectral density. From there I was able to find the frequency of the sinusoid due to the peaks. I am unsure how I could go about finding the amplitude, as I am unsure how I can go about normalizing the PSD or Autocorrelation function.

ps this is all done in matlab, thanks!

• It is frequent question on dsp.stackexchange. For example, you can see useful answers in question link – SergV Apr 15 '16 at 4:45

If you know you have a single sinusoid in Gaussian noise, then a least squares type parametric estimator might be one method of estimating the amplitude. See DSP Related, specifically:

Use a flat top window and then perform an fft to obtain the peak and its gain, which will be the maximum of the two bins covering the signal.

You could try an approach that is inspired by a Lock-In-amplifier:

Assumptions

You have the correct frequency of the sine wave. I assume that you do also have the initial phase; therefore I will assume that without loss of generality this initial phase is 0.

Calculation of the amplitude, version 1

Multiply your signal by $\sin(2\pi\,ft)$ for a reasonable time, e.g. one period $T$ of your sine. The frequency $f = \frac{1}{T}$. Both is therefore known from your first analysis step.

Since your signal is $A\sin(2\pi\,ft) + N$, where N is the noise, you'll end up with

$\int_0^T (A \sin(2 \pi \, f t)^2 + N \sin(2 \pi\, f t)) \textrm{d}t = \frac{4 \pi A f T - A \sin(4 \pi \, f T) + 8 N \sin(\pi \, f T)^2}{8 \pi f}$

The sines in the result will be $0$, since $\sin(n\pi) = 0$ for integer values of $n$, and $fT = 1$. Hence,

$\int_0^T (A \sin(2 \pi \, f t)^2 + N \sin(2 \pi\, f t)) \textrm{d}t = \frac{4 \pi A f T}{8 \pi f} = \frac{AT}{2}$

which can be solved for $A$.

If you use multiple periods, i.e. an integration time $KT$ with an integer K, you would end up with $\frac{AKT}{2}$. Since each half-period constitutes a kind of averaging in this approach, and you do not integrate, but sum up the values in your discrete case, your results should become better when you take multiple periods into account.

Calculation of the amplitude, version 2

It is basically the same idea, but in another form. $X = \sin(2\pi\, ft)$ is a matrix with only one column, that holds the pure sine with the known frequency. $\vec{y}$ is your noisy measurement as a column vector. With the linear approach $XA + \vec{N} = \vec{y}$, you could estimate the amplitude $A$ in a linear least squares sense:

$A = \frac{1}{X^TX}X^T\vec{y}$.

Since $X$ is a one-column matrix, it essentially behaves like a vector. Hence, this equation simply is $A = \frac{<X,\vec{y}>}{<X,X>}$, where $<\cdot,\cdot>$ denotes the scalar product.

There is an IEEE Standard for Digitizing Waveform Recorders (IEEE Std 1057). It can be found under terms like three-, four- (or seven-)parameter sine fit, for instance: $$y = a + b \sin(c+d x)$$ MatlabCentral offers a well-received tool for that: Four-Parameter Sinefit. You can use Fourier or PSD estimates to intialize such algorithms.

Several additional works offer to estimates the bias and variance, for instance: Amplitude estimation using IEEE-STD-1057 three-parameter sine wave fit: Statistical distribution, bias and variance.

It can be intersting as well to consider the effect of signal quantization $Q$, and estimate thee related uncertainty, see for instance Uncertainty of the Estimates of Sine Wave Fitting of Digital Data in the Presence of Additive Noise.