1
$\begingroup$

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?

Improve this question
$\endgroup$

3 Answers 3

Reset to default
1
$\begingroup$

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.

Improve this answer
$\endgroup$
4
  • $\begingroup$ After the Hilbert Transform, I would pass the signal through a low pass filter, in order to further reduce the presence of residual noise. $\endgroup$
    – fpe
    Aug 28, 2015 at 10:37
  • $\begingroup$ 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. $\endgroup$
    – Richard Lyons
    Aug 30, 2015 at 0:21
  • $\begingroup$ 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]$. $\endgroup$
    – robert bristow-johnson
    Sep 13, 2015 at 3:43
  • $\begingroup$ 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. $\endgroup$
    – Richard Lyons
    Sep 14, 2015 at 13:43
1
$\begingroup$

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.)

Improve this answer
$\endgroup$
0
$\begingroup$ 
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;
    };

}
Improve this answer
$\endgroup$
4
  • $\begingroup$ Would you care to explain the code? $\endgroup$
    – JRE
    Sep 2, 2015 at 13:54
  • $\begingroup$ 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 $\endgroup$
    – mypaul
    Sep 2, 2015 at 14:01
  • $\begingroup$ hope it works :) i will try it myself soon. $\endgroup$
    – mypaul
    Sep 2, 2015 at 14:02
  • $\begingroup$ the function is called "phase delimiter" $\endgroup$
    – mypaul
    Sep 2, 2015 at 14:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.