# How can I extract a sinusoid from a very noisy signal for computing the phase shift?

I have done some measurements and I got two very noisy signals. I want to extract a certain sinusoid from the signals and I want to compute the phase shift between them.

I made some simulations and I saw that the noise modifies a lot the value of the phase shift.

How can I get a trustful value for the phase shift?

You are correct. Estimating the instantaneous phase of a noisy sinusoid is NOT easy. I suggest you design a narrow bandpass filter such that your sinusoid-of-interest is in the filter's passband. (The better the filter the more noise that will be eliminated.) Pass your two signals through the bandpass filter to generate filtered signals $x_1[n]$ and $x_2[n]$. Next, pass your $x_1[n]$ and $x_2[n]$ signals through a Hilbert transformer to generate $\hat{x}_1[n]$ and $\hat{x}_2[n]$. Create two analytic (complex) signals as:

$$z_1[n] = x_1[n] + j \, \hat{x}_1[n],$$

and $$z_2[n] = x_2[n] + j \, \hat{x}_2[n]$$

where

\begin{align} \hat{x}[n] & = \mathcal{H}\{ x[n] \} \\ & = \sum\limits_{m=-\infty}^{+\infty} \frac{1 - (-1)^{m}}{\pi \, m} x[n-m] \\ \end{align}

Next, compute two instantaneous phase sequences:

$$\phi_1[n] = \arg\{z_1[n]\}$$

and $$\phi_2[n] = \arg\{z_2[n]\}.$$

Finally, compare the instantaneous phase difference between the $\phi_1[n]$ and $\phi_2[n]$ sequences.

• After the Hilbert Transform, I would pass the signal through a low pass filter, in order to further reduce the presence of residual noise. – fpe Aug 28 '15 at 10:37
• Are you asking a question? Hopefully the bandpass filter will eliminate unwanted noise from your original real-valued signal. If your original real-valued signal is low in frequency, then yes, lowpass filtering the x1(n), x2(n), y1(n), and y2(n) signals may further reduce unwanted noise. – Richard Lyons Aug 30 '15 at 0:21
• i did some expanding of Rick's answer. of course, you can't do the summation as shown (it's gotta be made finite and, for real-time, you'll have to delay both $\hat{x}[n]$ and $x[n]$. – robert bristow-johnson Sep 13 '15 at 3:43
• Just so readers know, the above infinite-limits summation was not in my original answer. I wish it hadn’t been inserted into my answer. And yes, the x1[n] and x2[n] sequences in the z1[n] and z2[n] equations are the original x1[n] and x2[n] sequences delayed-in-time by the time delay of the Hilbert transformer. – Richard Lyons Sep 14 '15 at 13:43

we've been over this before, somewhere (maybe comp.dsp).

if you think that $x_1[n]$ and $x_2[n]$ have the virtually the same frequency and you want to measure phase angle between them, the simplest noise-immune method is to beat one against the other and the Hilbert transform of the other. so pick $x_1[n]$ to be the reference sinusoid.

using the same notation of Rick's above:

$$u[n] + j \, v[n] = x_2[n] \left( x_1[n] + j \, \hat{x}_1[n] \right)$$

$$\phi_2[n] - \phi_1[n] = \arg \{ u[n] + j \, v[n] \}$$

to denoise it, you gotta low-pass filter the phase difference $\phi_2[n] - \phi_1[n]$ to the extent that you need and dare. (too much LPFing will slow down or cover a fluctuation of phase that you might want to see. too little LPF will leave the phase difference noisy.)

float phis=0;
void phdelim(short* dwData, short* dwData2, DWORD dwLength, float fphi,
float fepsilon, BOOL* bOk, DWORD* dwLength2){
if(bDebug){
MessageBox(cvar.get_hwnd(),
"compare the phase between dwData-dwData2 and then\
computes the difference between them.\
if the difference crosses zero or it is smaller than\
fepsilon bOK=TRUE. BOk=FALSE if the difference > fepsilon\
in every instance.\n", "titlu phdelim", 0);
};

DWORD dwi=0;
DWORD dwk=0;
DWORD dwj=0;
BOOL bOKintern;

do{
if (dwData2[dwi]!=0)
if((1-fepsilon<dwData[dwi]/(float)dwData2[dwi])&&\
(dwData[dwi]/(float)dwData2[dwi]<=1+fepsilon)){
bOKintern=TRUE;
dwk++;
};

if (bOKintern==TRUE)
if (dwData2[dwi+1]!=0)
if((1-fepsilon<dwData[dwi+1]/(float)dwData2[dwi+1])&&\
(dwData[dwi+1]/(float)dwData2[dwi+1]<=1+fepsilon)){
//                  bOk[0]=TRUE;
if (dwk==dwi)
dwj=dwk;
};

bOKintern=FALSE;
dwi++;
}while (dwi<dwLength&&bOKintern==FALSE);

if (bOKintern==TRUE&&dwk==dwi){
bOk[0]=TRUE;
dwLength2[0]=dwLength;
phis=fphi;
};

if ((float)dwj/(float)dwi<fepsilon){
bOk[0]=TRUE;
dwLength2[0]=dwj;
}else{
dwLength2[0]=dwi*fepsilon;
};

}

• Would you care to explain the code? – JRE Sep 2 '15 at 13:54
• dwData and dwData2 is the feed of input values (of dwLength in size, each), dwphi is the phase difference between them, and bOk is the return value if the phases are not shifted by fepsilon ratio. If the dwDatas values are not in sync with a ratio of fepsilon (there exists a further modulation) then bOk is still TRUE but deLength2 is equal to the size of dwData/dwData2 that are in close phase. thx – mypaul Sep 2 '15 at 14:01
• hope it works :) i will try it myself soon. – mypaul Sep 2 '15 at 14:02
• the function is called "phase delimiter" – mypaul Sep 2 '15 at 14:03