how to find phase shift between two signals

robert bristow-johnson wrote:
> Andor wrote:

> > 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"?

The same way you project one vector on another in 2D or 3D
geometry. Sum the product of terms.

> how is it done in the DFT?

X(k) = sum _N x(n) * exp(jwn)

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

Normalization is taken care of by the polar form of the spectrum
coefficients X(k):

Cxy(k) = rx*ry *exp(j(phix - phiy))

No need to normalize.

> 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?

XCOR doesn't multiply the signals. It correlates them. This is
corresponds to multiplying the *spectra* of the signals.

Rune

Vladimir Vassilevsky wrote:
> 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

The comp.dsp regulars aren't more "clueless" than that both of you
find it worth your while to hang around, apparently...

Rune

Rune Allnor wrote:
> robert bristow-johnson wrote:
> > Andor wrote:
>
> > > 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"?
>
> The same way you project one vector on another in 2D or 3D
> geometry. Sum the product of terms.
>
> > how is it done in the DFT?
>
> X(k) = sum _N x(n) * exp(jwn)
>
> > (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.
>
> Normalization is taken care of by the polar form of the spectrum
> coefficients X(k):
>
> Cxy(k) = rx*ry *exp(j(phix - phiy))
>
> No need to normalize.

by having a quadrature component of the reference, the amplitude of
thereal part of that complex correlation is used to normalize the
amplitude of the imaginary part in the arctan(.) function to get the
angle.

but with that "single phase" method:

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

where you arccos(.) the mean of the above, you have to first normalize
with A1*A2/2 which you can get by computing the r.m.s. of each of x1(t)
and x2(t).  but you still have the problem that you don't know if x1 is
leading x2 or lagging by the same amount.  nevertheless, even though
this method requires some sort of filtering or averaging to remove the
cos(2*w*t + phi1+phi2) component, it uses no Hilbert transformer.

but using that more complete method where x1 goes through a Hilbert
transformer to get the quadrature component, the normalzation is taken
care of by the division in the arctan(.) function, so i agree.

r b-j

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
...
> > > Figure it out for yourself :)

that's persuasive.

> > > As far as patent is concerned,

why would anyone care?  this obviates your patent for sure (you haven't
demonstrated how it has relevance to this yet) and likely anyone elses
patent.  maybe Bell Labs patented this in the 50s, but i doubt it.
even so, that patent is so old that it wouldn't matter now.
...
> Using Matlab it should take less then one page of code, probably just
> 15-20 lines of some terse heavily vectorized matlab code.

how about 3 lines of MATLAB code that is trivially combined to 1 line
of MATLAB code.  it's below.  then i invite you to post to us your
solution that makes all other solutions obsolete.

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

i think, Dmitry, this is evidence of your cluelessness.  nobody else is
talking about looking for peaks in a cross-correlation.  we're talking
about cross-correlation at a lag of zero.  this is not a pitch
detection problem.

> Like I said, phase-shift measurement between two different signals is
> not covered by the patent claims. It's all yours :-)

how generous.



if exist('L') ~= 1,
    L = 16384
end;

if exist('f0') ~= 1,
    f0 = 500.2/L
end;

if exist('A1') ~= 1,
   A1 = 1.0e-5
end;

if exist('phi1') ~= 1,
    phi1 = 0.1
end;

if exist('A2') ~= 1,
   A2 = 1.0
end;

if exist('phi2') ~= 1,
    phi2 = 1.56
end;

t = linspace(-L/2, L/2-1, L);

if exist('noise_amp') ~= 1,
    noise_amp = 0.01
end;

if exist('noise_pole') ~= 1,
    noise_pole = 0.1
end;

% or design your own noise signal
noise = noise_amp*(rand([1 L]) - 0.5);
noise = filter([(1-noise_pole) 0], [1 -noise_pole], noise);

x1 = A1*cos(2*pi*f0*t + phi1);

x2 = A2*cos(2*pi*f0*t + phi2) + noise;

% the above generates data,  you can rerun this with different
% parameters: A1, A2, phi1, phi2, noise_amp, noise_pole
% the code below solves the problem

  x1_h = hilbert(x1);

  X_corr = sum( x1_h .* x2 .* hanning(length(x1))' );

  theta = -angle(X_corr);

% let's look at the difference
phi2 - phi1 - theta

(end of file)



this could be one line of MATLAB trivially:

theta = -angle(sum( hilbert(x1) .* x2 .* hanning(length(x1))' ));

now, Dmitry, it's time for you to wow us with your brilliant solution
to obsolete all other solutions that takes 15-20 lines of code.  i
invite you to use the same code to generate test data (or if your noise
is different, post the code that generates it, this is simple LPF white
noise with adjustable amplitude).

waiting with bated breath.

r b-j

decious wrote:
> 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 = sig_iq_1 .* conj( sig_iq_2 );
> delta_phase = angle( sig_xcorr );
>
> to prove it, look at the complex signals, and work throught he math, it
> falls out very easily.

yup just like this.  but the OP probably didn't have complex sinusoids,
but only real ones.  this is why you need to apply the Hilbert
transform to one of them.  and you don't need to compute the complex
analytic signal for both, only one:

delta_phase = -angle(sum( hilbert(x1) .* x2 .* hanning(length(x1))' ));

windowing down the edges is useful in case your arrays don't have an
integer number of cycles (or half cycles) in them.  just to lose any
edge effects.

ya know, i just figgered out that i can use the fact that the
hanning(.) is a column vector so the call to sum() might not be
necessary (i can dot product it).

delta_phase = -angle( (hilbert(x1) .* x2) * hanning(length(x1)) );

that is a one line MATLAB answer to the OP.

i still wanna see our two Russian friends demonstrate their brilliant
solutions.

r b-j

Vladimir Vassilevsky wrote:
> 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.

(birds of a feather?)

yeah, he just nailed us, V.

what's your brilliant and simple solution?  we know what Dmitry's is.
plug his state-space pitch detection algorithm as the be all and end
all of any measurement problem, even problems of little similarity to
the pitch detection problem.

have any more smoke to blow?

r b-j

fizteh89 wrote:
> 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)

hey Dmitry, i ran that number:

http://patft.uspto.gov/netacgi/nph-Parser?TERM1=20030088401&Sect1=PTO1&Sect2=HITOFF&d=PALL&p=1&u=%2Fnetahtml%2FPTO%2Fsrchnum.htm&r=0&f=S&l=50

and it says:

    Searching US Patents Collection...

    Results of Search in US Patents Collection db for:
    PN/20030088401: 0 patents.

    No patents have matched your query

can you give us a link that works?

i already read your paper and still never got any of your MATLAB code
to work.  but i still invite you to post a simple MATLAB program the
demos your pitch detection.  we have your audio file: "She had your
dark suit in greasy wash water all year."

i have lot's of other audio files of singing voices and monophonic
instruments (or polyphonic that are deliberately played monophonically,
like lead guitar).  i would like to test them on your pitch detection
algorithm.  i have MATLAB 5 on a PC also.  do you have a MATLAB program
that will work simply that i can test with other sound files?

r b-j

robert bristow-johnson wrote:
> fizteh89 wrote:
> > 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)
>
> hey Dmitry, i ran that number:
>
> http://patft.uspto.gov/netacgi/nph-Parser?TERM1=20030088401&Sect1=PTO1&Sect2=HITOFF&d=PALL&p=1&u=%2Fnetahtml%2FPTO%2Fsrchnum.htm&r=0&f=S&l=50
>
> and it says:
>
>     Searching US Patents Collection...
>
>     Results of Search in US Patents Collection db for:
>     PN/20030088401: 0 patents.
>
>     No patents have matched your query
>
> can you give us a link that works?

i've also tried the search with APN (Application Patent Number):
APN/20030088401

and still get zilch.

is there an actual link you can get us?

BTW, i know patents can sometimes take a couple of years to be finally
awarded (so it isn't still a Patent Application, but the real thing).
you've been at this for a few years.  isn't it curious that you are
still quoting us a US Patent Application Publication number (that i
can't get to work, anyway)?

r b-j

hi all,
       how about doing correlation on each and finding out phase (blind
method) and then doing the difference.

particle

robert bristow-johnson wrote:
> robert bristow-johnson wrote:
> > fizteh89 wrote:
> > > 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)
> >
> > hey Dmitry, i ran that number:
> >
> > http://patft.uspto.gov/netacgi/nph-Parser?TERM1=20030088401&Sect1=PTO1&Sect2=HITOFF&d=PALL&p=1&u=%2Fnetahtml%2FPTO%2Fsrchnum.htm&r=0&f=S&l=50
> >
> > and it says:
> >
> >     Searching US Patents Collection...
> >
> >     Results of Search in US Patents Collection db for:
> >     PN/20030088401: 0 patents.
> >
> >     No patents have matched your query
> >
> > can you give us a link that works?
>
> i've also tried the search with APN (Application Patent Number):
> APN/20030088401
>
> and still get zilch.
>
> is there an actual link you can get us?

okay i found it.

http://appft1.uspto.gov/netacgi/nph-Parser?Sect1=PTO1&Sect2=HITOFF&d=PG01&p=1&u=%2Fnetahtml%2FPTO%2Fsrchnum.html&r=1&f=G&l=50&s1=%2220030088401%22.PGNR.&OS=DN/20030088401&RS=DN/20030088401

typical inpenetrable patent drivel (or patent-speak), i've been there
before.   but your ICASSP paper is telling:

http://www.soundmathtech.com/pitch/paper.html

and posted these observations:

http://groups.google.com/group/comp.speech.research/browse_frm/thread/47c863a624aac888/1897d8af3742992a?lnk=st&q=&rnum=1&hl=en#1897d8af3742992a

this guy isn't very impressed either:

http://www2.cs.uregina.ca/~gerhard/publications/TRdbg-Pitch.pdf

i've still never been able to get your pdemo.zip to work on MATLAB 5
(Mac or PC).  can you post a simple MATLAB program that works and that
we can actually run our own soundfiles?


r b-j