how to find phase shift between two signals

"fizteh89" <dt@soundmathtech.com> wrote in news:1156265930.842902.75290
@m73g2000cwd.googlegroups.com:

> Cross-correlation is not the method of choice anymore.
> 
> State-space embedding is...
> 
> Read US Patent Application Pub. 20030088401 at
> http://www.uspto.gov/patft
> (Already allowed with all of the original and some new claims)
> 
> Also, ICASSP 2002 paper and Matlab demo can be found at:
> http://www.soundmathtech.com/pitch
> 
> Works equally well for measuring the period of any periodic signal or
> just measuring the phase shift between two signals
> Figure it out for yourself :)
> 
> As far as patent is concerned, you don't need to worry as long as you
> stay away from measuring the actual period/frequency of your sine wave
> and just measure the phase shift between two different sine waves.
> 
> You can safely disregard all comments by comp.dsp "experts": some of
> them are clueless, others are jealous :)
> 
> Best Regards
> 

You are kidding, right?  The guy is describing a single-frequency 
sinusoid in, and the same frequency out, at steady state.  He wants a 
phase identification algorithm that is quick and easy to implement, 
probably for someone without a signal processing background.  None of the 
obstacles to traditional methods that the algorithm was designed to 
overcome apply.

-- 
Scott
Reverse name to reply

Scott Seidman wrote:
> "fizteh89" <dt@soundmathtech.com> wrote in news:1156265930.842902.75290
> @m73g2000cwd.googlegroups.com:
>
> > Cross-correlation is not the method of choice anymore.

when the only tool you have is a hammer, every problem looks like a
nail.

> You are kidding, right?

i dunno if that can be determined.  i dunno how a pitch-detection
algorithm (even the "Mother of all Pitch-Detection Algorithms" or "The
Pitch-Detection Algorithm to end all pitch-detection algorithms") is
used to measure phase shift between input and output of a supposed LTI
system, but Dmitry evidently knows how.

>  The guy is describing a single-frequency
> sinusoid in, and the same frequency out, at steady state.  He wants a
> phase identification algorithm that is quick and easy to implement,
> probably for someone without a signal processing background.

and maybe even robust and noise-immune.

looks like i need to append my other post.

r b-j

Scott Seidman wrote:
> "fizteh89" <dt@soundmathtech.com> wrote in news:1156265930.842902.75290
> @m73g2000cwd.googlegroups.com:
>
> > Cross-correlation is not the method of choice anymore.
> >
> > State-space embedding is...
> >
> > Read US Patent Application Pub. 20030088401 at
> > http://www.uspto.gov/patft
> > (Already allowed with all of the original and some new claims)
> >
> > Also, ICASSP 2002 paper and Matlab demo can be found at:
> > http://www.soundmathtech.com/pitch
> >
> > Works equally well for measuring the period of any periodic signal or
> > just measuring the phase shift between two signals
> > Figure it out for yourself :)
> >
> > As far as patent is concerned, you don't need to worry as long as you
> > stay away from measuring the actual period/frequency of your sine wave
> > and just measure the phase shift between two different sine waves.
> >
> > You can safely disregard all comments by comp.dsp "experts": some of
> > them are clueless, others are jealous :)
> >
> > Best Regards
> >
>
> You are kidding, right?  The guy is describing a single-frequency
> sinusoid in, and the same frequency out, at steady state.  He wants a
> phase identification algorithm that is quick and easy to implement,
> probably for someone without a signal processing background.  None of the
> obstacles to traditional methods that the algorithm was designed to
> overcome apply.
>
> --
> Scott
> Reverse name to reply

No, I am not kidding...

It's about as quick and easy as implementing cross-correlation (OK, a
little harder,
but this is not the Poincar=E9 conjecture after all :-)

Using Matlab it should take less then one page of code, probably just
15-20 lines of some terse heavily vectorized matlab code.

But, of course, you are welcome to use cross-correlation if you prefer
to locate some very wide peaks as opposed to locating very narrow
peaks (as narrow as you can possibly get them with your given sampling
rate) produced by the new method...
Like I said, phase-shift measurement between two different signals is
not covered by the patent claims. It's all yours :-)

Andor wrote:
> robert bristow-johnson wrote:
> >
> > let:
> >
> >     x1(t) = A1 * cos(w*t + phi1)
> > and
> >     x2(t) = A2 * cos(w*t + phi2)
> >
> >     x1(t)*x2(t) = A1*A2 * cos(w*t + phi1)*cos(w*t + phi2)
> >
> >        = A1*A2/2 * ( cos(2*w*t + phi1+phi2) + cos(phi1-phi2) )
> >
> > if you low pass filter (average) the result the high frequency term
> > goes away.  another way to compute the mean is to integrate over one
> > period and divide by the length of that period.
> >
> >     A1*A2/2 * cos(phi1-phi2) = mean{ x1(t)*x2(t) }
> >
> >     cos(phi1-phi2) = 2 * mean{ x1(t)*x2(t) } / A1*A2
> >
> >     |phi1-phi2| = arccos( 2 * mean{ x1(t)*x2(t) } / A1*A2 )
> >
> > so you'll still have to find a way to determine the amplitudes A1 and
> > A2, but i'll bet you can figure that out.
> >
> > one last problem, this doesn't tell you if phi1-phi2 is positive or
> > negative (because cos(phi1-phi2) = cos(phi2-phi1)), if x1(t) is leading
> > or lagging x2(t).  the really fancy mathematical method to determine
> > who is leading whom is to compute the Hilbert transform of either x1(t)
> > or x2(t) (whichever you think is the reference signal) and do this same
> > correlation.  the math is very similar.
>
> That sounds awfully complicated. I would just project the two sine
> waves (input and output) onto a complex phasor at the given frequency,
> sort of like computing a single bin of a DFT, except that the frequency
> is not a multiple of the sampling frequency (make sure that enough
> periods are covered in the projection to reduce edge effects).

how do you *do* that "projection"?  how is it done in the DFT?  (you
*multiply*.)  if you use, as your reference frequency a copy of the
input x1(t), how do you normalize it?  if you don't normalize it what
will happen in your math.

Andor, what i illustrated above, *is* essentially the same thing you
are suggesting, but with some of the holes filled in.

Scott Seidman wrote:
>
> In my experience, if no other info is needed, the xcorr is probably the
> quickest to implement.

yes, that's how you multiply the two signals together.  but what do you
do about the amplitudes of the input and output?

> > the really fancy mathematical method to determine
> >who is leading whom is to compute the Hilbert transform of either x1(t)
> >or x2(t) (whichever you think is the reference signal) and do this same
> >correlation.  the math is very similar.
>
> Assuming you're talking about using the Hilbert Transform of the product of
> the two sinusoids, you can't do this when the two sinusoids are of the same
> frequency.  More accurately, you can, of course, take the Hilbert
> Transform, but the phase you get out won't be what you want it to be.
> Aliasing, you know.

no, no, no, no, no...  there was (a lazy) reason i presented it in
continuous time.  i said (and meant) applying the Hilbert Transform to
either x1(t) OR to x2(t), whichever one is considered the reference
sinusoid.  lte's say it's x1(t), then

    x1'(t) = Hilbert{ x1(t) } =  A1 * sin(w*t + phi1)

so

    x1(t)*x2(t) = A1*A2/2*( cos(2*w*t + phi1+phi2) + cos(phi1-phi2) )

and

    x1'(t)*x2(t) = A1*A2/2*( sin(2*w*t + phi1+phi2) + sin(phi1-phi2) )

so, if you LPF (or average, or however you'll do it) out the high
frequency terms of both you have:

    mean{ x1(t)*x2(t) } = A1*A2/2 * cos(phi2-phi1)

    mean{ x1'(t)*x2(t) } = -A1*A2/2 * sin(phi2-phi1)

and now, unambiguously,

    ph2-phi1 = arg{ mean{ x1(t)*x2(t) } - j*mean{ x1'(t)*x2(t) } }

or

    phi2-ph1 =  sgn( Im )*pi/2 - arctan(Re/Im)

where Re = mean{ x1(t)*x2(t) }
   and Im = -mean{ x1'(t)*x2(t) }

that's what i meant.  it's robust for any phase difference from -pi to
+pi, it's pretty noise and DC immune, and you don't need to know the
frequency or look for zero-crossings or the like.

if *you* are the one generating the reference x1(t) (say using a NCO or
DDS or whatever we call a table-lookup sine wave generator), you can
generate x1'(t) trivially as another output (may just keep this signal
in the mind of the DSP if it is also inputting the result, x2(t)), so
you need not compute a Hilbert transform of x1(t), just synthesize it.

this is NOT exotic.  it's the normal, simple, and robust way of doing
this sort of thing.

r b-j

robert bristow-johnson wrote:

> Andor wrote:
> > robert bristow-johnson wrote:
> > >
> > > let:
> > >
> > >     x1(t) = A1 * cos(w*t + phi1)
> > > and
> > >     x2(t) = A2 * cos(w*t + phi2)
> > >
> > >     x1(t)*x2(t) = A1*A2 * cos(w*t + phi1)*cos(w*t + phi2)
> > >
> > >        = A1*A2/2 * ( cos(2*w*t + phi1+phi2) + cos(phi1-phi2) )
> > >
> > > if you low pass filter (average) the result the high frequency term
> > > goes away.  another way to compute the mean is to integrate over one
> > > period and divide by the length of that period.
> > >
> > >     A1*A2/2 * cos(phi1-phi2) = mean{ x1(t)*x2(t) }
> > >
> > >     cos(phi1-phi2) = 2 * mean{ x1(t)*x2(t) } / A1*A2
> > >
> > >     |phi1-phi2| = arccos( 2 * mean{ x1(t)*x2(t) } / A1*A2 )
> > >
> > > so you'll still have to find a way to determine the amplitudes A1 and
> > > A2, but i'll bet you can figure that out.
> > >
> > > one last problem, this doesn't tell you if phi1-phi2 is positive or
> > > negative (because cos(phi1-phi2) = cos(phi2-phi1)), if x1(t) is leading
> > > or lagging x2(t).  the really fancy mathematical method to determine
> > > who is leading whom is to compute the Hilbert transform of either x1(t)
> > > or x2(t) (whichever you think is the reference signal) and do this same
> > > correlation.  the math is very similar.
> >
> > That sounds awfully complicated. I would just project the two sine
> > waves (input and output) onto a complex phasor at the given frequency,
> > sort of like computing a single bin of a DFT, except that the frequency
> > is not a multiple of the sampling frequency (make sure that enough
> > periods are covered in the projection to reduce edge effects).
>
> how do you *do* that "projection"?  how is it done in the DFT?  (you
> *multiply*.)

Standard dot-product. Mutliply and accumulate (integrate). Upon reading
your post again, this is probably what you meant with

> if you low pass filter (average) the result the high frequency term
> goes away

I just skimmed over your post without trying to grasp the details.

> if you use, as your reference frequency a copy of the
> input x1(t), how do you normalize it?

I was suggesting to project both signals onto a normalized phasor (this
way, you get relative phase shift and amplitude). Alternatively, one
can indeed project one signal onto the other - the autocorrelation
approach. For that, one needs to know at least one amplitude but not
the frequency (which one needs to know for my approach).

> if you don't normalize it what
> will happen in your math.
>
> Andor, what i illustrated above, *is* essentially the same thing you
> are suggesting, but with some of the holes filled in.

Your derivation just seemed to arrive from a strange angle. I like to
the geometric approach to signal processing, it's very intuitive for
me.

Regards,
Andor

> "The Pitch-Detection Algorithm to end all pitch-detection algorithms"

:-D

> looks like i need to append my other post.

Maybe you should patent it.

"robert bristow-johnson" <rbj@audioimagination.com> wrote in 
news:1156274201.907780.75130@i42g2000cwa.googlegroups.com:

> yes, that's how you multiply the two signals together.  but what do you
> do about the amplitudes of the input and output?
> 

You can pull out the magnitude-- you can cross correlate the input and 
output  each with the sine and the cosine, pulling out the real and 
imaginary parts of both.

More importantly, though, the original poster did not seem to care about 
the magnitude.  I wouldn't want to create requirements that don't exist! --
though, likely, I'm sure the OP will find he wants the amplitude as well.  
Personally, I'd go with the fft.

-- 
Scott
Reverse name to reply

this seems pretty straight forward,
take IQ data for Signal A, multiple it by the conjugate of signal B,
take the angle of the output.

so for the MATLAB types
sig_xcorr =3D sig_iq_1 .* conj( sig_iq_2 );
delta_phase =3D angle( sig_xcorr );

to prove it, look at the complex signals, and work throught he math, it
falls out very easily.


fizteh89 wrote:
> Scott Seidman wrote:
> > "fizteh89" <dt@soundmathtech.com> wrote in news:1156265930.842902.75290
> > @m73g2000cwd.googlegroups.com:
> >
> > > Cross-correlation is not the method of choice anymore.
> > >
> > > State-space embedding is...
> > >
> > > Read US Patent Application Pub. 20030088401 at
> > > http://www.uspto.gov/patft
> > > (Already allowed with all of the original and some new claims)
> > >
> > > Also, ICASSP 2002 paper and Matlab demo can be found at:
> > > http://www.soundmathtech.com/pitch
> > >
> > > Works equally well for measuring the period of any periodic signal or
> > > just measuring the phase shift between two signals
> > > Figure it out for yourself :)
> > >
> > > As far as patent is concerned, you don't need to worry as long as you
> > > stay away from measuring the actual period/frequency of your sine wave
> > > and just measure the phase shift between two different sine waves.
> > >
> > > You can safely disregard all comments by comp.dsp "experts": some of
> > > them are clueless, others are jealous :)
> > >
> > > Best Regards
> > >
> >
> > You are kidding, right?  The guy is describing a single-frequency
> > sinusoid in, and the same frequency out, at steady state.  He wants a
> > phase identification algorithm that is quick and easy to implement,
> > probably for someone without a signal processing background.  None of t=
he
> > obstacles to traditional methods that the algorithm was designed to
> > overcome apply.
> >
> > --
> > Scott
> > Reverse name to reply
>
> No, I am not kidding...
>
> It's about as quick and easy as implementing cross-correlation (OK, a
> little harder,
> but this is not the Poincar=E9 conjecture after all :-)
>
> Using Matlab it should take less then one page of code, probably just
> 15-20 lines of some terse heavily vectorized matlab code.
>
> But, of course, you are welcome to use cross-correlation if you prefer
> to locate some very wide peaks as opposed to locating very narrow
> peaks (as narrow as you can possibly get them with your given sampling
> rate) produced by the new method...
> Like I said, phase-shift measurement between two different signals is
> not covered by the patent claims. It's all yours :-)

"robert bristow-johnson" <rbj@audioimagination.com> wrote in 
news:1156274201.907780.75130@i42g2000cwa.googlegroups.com:

> 
> this is NOT exotic.  it's the normal, simple, and robust way of doing
> this sort of thing.
> 


That IS pretty nifty.

-- 
Scott
Reverse name to reply

fizteh89 wrote:


> 
> You can safely disregard all comments by comp.dsp "experts": some of
> them are clueless, others are jealous :)
> 

You nailed them! 33 wizards solwing a great problem of 2 x 2 = 4
Thank you for the excellent comment.

VLV