Algorithm for finding best matching sin wave to input signal

Question

Question

Suppose I take a signal $y(t) = \cos(2\pi \cdot 4.1t + 2)$ and I sample it uniformly up to 2 seconds.

Given only these time samples, how would I design an algorithm that finds the pair of values (freq=4.1 Hz, phase=2 Rad)?

Notes

I do not want the DFT. The DFT assumes the signal is circular at the boundaries and works by using many values in many frequency bins. What I want is to find the best singular phase and frequency that minimise mean squared error between the input data and the reconstructed signal.

For reference, the DFT power spectrum at 400 samples gives this:

My Attempt

I attempted a grid search approach. I choose 1000 phases $\phi$ in $[0, 2\pi)$ and 1000 frequencies $f$ up to the nyquist frequency ($100$ Hz)

Once I had my grid of cos waves, I evaluated the signal at each pair ($\phi$, $f$)

$$y_\text{guess}[n] = \cos(2\pi \cdot f \cdot t[n] + \phi)$$

And calculated the mean square error between the guess and the original signal $$MSE = \frac{\sum_n(y[n] - y_\text{guess}[n])^2}{N}$$ where $N$ is the number of samples.

Running this algorithm, I found the pair with the lowest MSE was ($\phi$ = 2.0357..., $f$=4.0938...) with an MSE of 0.000247.

Those are pretty good values. Problem is, it's really slow. I can't help but feel like there's probably a better, more elegant way to do what I want to do.

I will also mention that while it isn't in the scope of this question, once I have figured this out this algorithm I would like an approach where I try to reconstruct a signal with N sin waves (say, N=5 would find the parameters for 5 sin waves that best minimise error).

Answer 1

Well, in this noiseless case, you can just estimate the frequency directly: take the $\arccos(y(t))$ of each sample, find the difference between consecutive samples: that's how much the phase changes per sample, i.e. your $f$.

You might sometimes need to "unwrap" the phase: if the absolute of the difference between two consecutive samples is $\left\lvert\arccos(y(t+1))-\arccos(y(t))\right\rvert > \pi$, then you just need to add/subtract $2\pi$ to keep things consistent.

Your phase $\varphi$ is just the arctan of the first sample, since $t=0$ at that point.

In the noisy case, you can average after unwrapping. Or, if this is for a continous signal, and up to medium SNRs, use a PLL to lock onto the signal, and read the frequency of the signal from the frequency that the PLL oscillator gets regulated to.

Under much noise, or for a sum of more than one frequency component, use a parametric spectral estimator. An FFT / DFT isn't your tool of choice there – it'll only be exact on its own for frequencies that land exactly on the grid defined by the sampling rate divided by the length of the DFT, and interpolative methods are sensitive to noise.

would like an approach where I try to reconstruct a signal with N sin waves (say, N=5 would find the parameters for 5 sin waves that best minimise error).

ESPRIT sounds like the right tool here in reality: you know it's pure tones, and you know the number of tones. The only thing missing is noise - completely without noise, this algorithm becomes unstable. It's a projective method that maps the signal into a signal and a noise subspace - and by properties of orthogonality, that minimizes the squared error.

Answer 2

When starting with a DFT or FFT, one can use Sinc interpolation of the complex components with binary search successive approximation to estimate a narrow-band or single frequency spectral peak between DFT result bins. Or, if computational efficiency isn't required, you can zero-pad a FFT by a factor of 1000X+ to get that much finer a bin resolution.

Answer 3

Use the inverse solutions to the sine DFT solution:

• Frequency
• Amplitude & Phase

If noiseless, recovery is exact. If noisy, it's competitively robust (maybe state of the art in compute-limited settings) and improved with greater $N$. Frequency is most robust, phase least.

Examples (see linked posts for more + explanations):

The DFT assumes the signal is circular at the boundaries

No. Just another case of "DFT periodicity" confusion.