retrieve phase of an extrapolated periodic 1D signal

Question

I have a periodic signal on a given grid, say (in matlab format):

t = 1:30;
omega = 2*pi/18.431;
phi = -pi+2*pi*rand(1);   % a random phase [-pi,pi]
x = sin(omega*t+phi);     % the signal
x = x+0.5*rand(1,length(x)); % add some noise

Now I want to retrieve the phase phi. There are a few approaches to that, one is to use fft. The problem is that my grid is not good enough to pick that frequency, hence the phase. (another way is to fit to a sin, but that takes way too long if I need to do it 1e6 times)
I tried to implement the "Autoregressive all-pole model parameters — Burg's method" that is mention in an answer to this question. I manage to pick up the frequency quite well, but fail to pick up the phase.

I was trying to use angle(fft_vec(idx)), where fft_vec is the spectrum obtained after extrapolating the original signal (using arburg), and idx is the index where the frequency omega is picked up. What am I doing wrong? and how can I manage to get that phase by other means?

Answer 1

if you know, in advance, the frequency of the sinusoid, so only the phase (and perhaps the amplitude) is unknown, the simplest noise-immune method i know of to get amplitude and phase is (i'm gonna do this in continuous time and let you translate it to discrete time and eventually to MATLAB):

$$ x(t) = \sin(\omega t + \phi) + n(t) $$

compute

$$ \begin{align} a & = \int\limits_0^{\frac{2 \pi}{\omega}} x(t) \sin(\omega t) dt \\ & = \int\limits_0^{\frac{2 \pi}{\omega}} \left(\sin(\omega t + \phi) + n(t) \right) \sin(\omega t) dt \\ & = \int\limits_0^{\frac{2 \pi}{\omega}} \sin(\omega t + \phi) \sin(\omega t) dt + \int\limits_0^{\frac{2 \pi}{\omega}} n(t) \sin(\omega t) dt \\ & = \frac{1}{2} \int\limits_0^{\frac{2 \pi}{\omega}} \cos(\phi) dt - \frac{1}{2} \int\limits_0^{\frac{2 \pi}{\omega}} \cos(2\omega t + \phi) dt + \int\limits_0^{\frac{2 \pi}{\omega}} n(t) \sin(\omega t) dt \\ & \approx \frac{1}{2} \int\limits_0^{\frac{2 \pi}{\omega}} \cos(\phi) dt \\ & = \frac{\pi}{\omega} \cos(\phi) \\ \end{align} $$

and

$$ \begin{align} b & = \int\limits_0^{\frac{2 \pi}{\omega}} x(t) \cos(\omega t) dt \\ & = \int\limits_0^{\frac{2 \pi}{\omega}} \left(\sin(\omega t + \phi) + n(t) \right) \cos(\omega t) dt \\ & = \int\limits_0^{\frac{2 \pi}{\omega}} \sin(\omega t + \phi) \cos(\omega t) dt + \int\limits_0^{\frac{2 \pi}{\omega}} n(t) \cos(\omega t) dt \\ & = \frac{1}{2} \int\limits_0^{\frac{2 \pi}{\omega}} \sin(\phi) dt + \frac{1}{2} \int\limits_0^{\frac{2 \pi}{\omega}} \sin(2\omega t + \phi) dt + \int\limits_0^{\frac{2 \pi}{\omega}} n(t) \cos(\omega t) dt \\ & \approx \frac{1}{2} \int\limits_0^{\frac{2 \pi}{\omega}} \sin(\phi) dt \\ & = \frac{\pi}{\omega} \sin(\phi) \\ \end{align} $$

in both expressions (on the bottom line) we know the middle integral is zero and expect the integral on the right to tend toward zero because $n(t)$ is bipolar and ain't correlated with anything.

i think, then, that you'll find the phase $\phi$ to be:

$$ \phi = \arg\left\{a + jb \right\} $$

which is, if $a>0$

$$ \phi = \arctan\left(\frac{b}{a} \right) $$

otherwise, you might have to add or subtract $\pi$ to that expression.