Estimating phase of Noisy Sinusoidal Signal

I am developing ultrasound based imaging system for research purposes.
I am looking for a very accurate algoritham/method for estimating a phase
of a sinusoidal signals corrupted by noise. Frequency of the sinusoid is a
prior information (fairly accurate), I do not care about the magnitude of a
sinusoid. Sinusiod is excited by mechanical resonator, it will have a
transition cycles before it reaches a steady state.  I end up with several
sinusoids which are delayed with respect to each other. I am interested in
estimating a phase shift between them. Sinusoids are in the frequency range
of 25-200 Hz. Sampling frequency is lower then 5 kHz.  I was applying next
algorithm:

Choose region of interest with respect to most delayed sinusoid (to get
rid of the transient effect). I would get complex envelope (make IQ data
using Hilbert transform), and then use R(1) (lag one autocorrelation
estimator) to estimate phase shift. This is a stanard Doppler processing
(radars/ultrasound) for phae shift estimation. This works but is not that
accurate. Main reason is that I am trying to preserve resolution in phase
estimation, thus I am only using ensembles of 2 signals for estimating a
phase shift( larger the ensemble and higher the SNR lower the variance of
estimation).

I would appreciate any interesting suggestion/discussion.


Thanks,

ultrasonics

On Thu, 05 Feb 2009 21:21:30 -0600, ultrasonics wrote:

> I am developing ultrasound based imaging system for research purposes. I
> am looking for a very accurate algoritham/method for estimating a phase
> of a sinusoidal signals corrupted by noise. Frequency of the sinusoid is
> a prior information (fairly accurate), I do not care about the magnitude
> of a sinusoid. Sinusiod is excited by mechanical resonator, it will have
> a transition cycles before it reaches a steady state.  I end up with
> several sinusoids which are delayed with respect to each other. I am
> interested in estimating a phase shift between them. Sinusoids are in
> the frequency range of 25-200 Hz. Sampling frequency is lower then 5
> kHz.  I was applying next algorithm:
> 
> Choose region of interest with respect to most delayed sinusoid (to get
> rid of the transient effect). I would get complex envelope (make IQ data
> using Hilbert transform), and then use R(1) (lag one autocorrelation
> estimator) to estimate phase shift. This is a stanard Doppler processing
> (radars/ultrasound) for phae shift estimation. This works but is not
> that accurate. Main reason is that I am trying to preserve resolution in
> phase estimation, thus I am only using ensembles of 2 signals for
> estimating a phase shift( larger the ensemble and higher the SNR lower
> the variance of estimation).
> 
> I would appreciate any interesting suggestion/discussion.
> 
> 
> Thanks,
> 
> ultrasonics

Since you already know the frequency, why not correlate the sine wave 
with its I & Q prototypes over some integer number of cycles, then find 
the polar angle of the resulting complex number?

I do this sort of thing all the time, and it is easy to make quite 
accurate; if you know the frequency exactly it's very likely the optimal 
phase estimator in the presence of additive white noise.

-- 
http://www.wescottdesign.com

On Feb 5, 11:09&#4294967295;pm, Tim Wescott <t...@seemywebsite.com> wrote:
> On Thu, 05 Feb 2009 21:21:30 -0600, ultrasonics wrote:
> > I am developing ultrasound based imaging system for research purposes. I
> > am looking for a very accurate algoritham/method for estimating a phase
> > of a sinusoidal signals corrupted by noise. Frequency of the sinusoid is
> > a prior information (fairly accurate), I do not care about the magnitude
> > of a sinusoid. Sinusiod is excited by mechanical resonator, it will have
> > a transition cycles before it reaches a steady state. &#4294967295;I end up with
> > several sinusoids which are delayed with respect to each other.

do you know anything about the range of phase shift between any pair?
if it is between 0 and pi, then you won't need to H.T. the reference
sinusoid to generate a quadrature reference sine.  if you have no idea
what the phase relationship is, you need both the I and the Q.

> > I am
> > interested in estimating a phase shift between them. Sinusoids are in
> > the frequency range of 25-200 Hz. Sampling frequency is lower than 5
> > kHz. &#4294967295;I was applying next algorithm:
>
> > Choose region of interest with respect to most delayed sinusoid (to get
> > rid of the transient effect).

in other words, your region of interest is the most current data.
this is what you get because the LP filters used for the correlation
Tim refers to tend to forget about input that happened longer ago.

> > I would get complex envelope (make IQ data
> > using Hilbert transform), and then use R(1) (lag one autocorrelation
> > estimator) to estimate phase shift. This is a standard Doppler processing
> > (radars/ultrasound) for phase shift estimation. This works but is not
> > that accurate.

we were talking about that in the other thread also.

> > Main reason is that I am trying to preserve resolution in
> > phase estimation, thus I am only using ensembles of 2 signals for
> > estimating a phase shift( larger the ensemble and higher the SNR lower
> > the variance of estimation).
...
> Since you already know the frequency, why not correlate the sine wave
> with its I & Q prototypes over some integer number of cycles, then find
> the polar angle of the resulting complex number?

isn't this what we were just talking about in this thread?:

 http://groups.google.com/group/comp.dsp/msg/0008e4899b8d85ee

heck, i even put a "noise" term in it.

depending on what kinda math library the OP has to work with, there
may be some effort needed to perform the arctan() function.

r b-j

robert bristow-johnson  <rbj@audioimagination.com> wrote:

>do you know anything about the range of phase shift between any pair?
>if it is between 0 and pi, then you won't need to H.T. the reference
>sinusoid to generate a quadrature reference sine.  if you have no idea
>what the phase relationship is, you need both the I and the Q.

I'm not seeing why you would ever need to Hilbert Transform
any signal for this problem.

Steve

On Feb 6, 12:24�am, spop...@speedymail.org (Steve Pope) wrote:
> robert bristow-johnson &#4294967295;<r...@audioimagination.com> wrote:
>
> >do you know anything about the range of phase shift between any pair?
> >if it is between 0 and pi, then you won't need to H.T. the reference
> >sinusoid to generate a quadrature reference sine. &#4294967295;if you have no idea
> >what the phase relationship is, you need both the I and the Q.
>
> I'm not seeing why you would ever need to Hilbert Transform
> any signal for this problem.

because of two reasons.  (if either was not the case, he wouldn't have
to.)

first, the OP is not generating the reference sinusoid.  if he was, he
could just as well generate a synced sine that was 90 degrees outa
phase (like in the other current thread about this same topic).

second, because of quadrant ambiguity (if he already knows he's in
quadrants I or II, it's no problem, he does an arccos() instead of an
arctan()), he needs to correlate the sinusoid of unknown phase against
both the reference sinusoid and the quadrature component of the
reference sinusoid.  if you can settle the quadrant ambiguity more
easily (say, by looking at zero crossings, but remember this is
noisy), then fine - no need for HT.

anyway, that's my spin on it.

r b-j

robert bristow-johnson wrote:
(snip)

> first, the OP is not generating the reference sinusoid.  if he was, he
> could just as well generate a synced sine that was 90 degrees outa
> phase (like in the other current thread about this same topic).

> second, because of quadrant ambiguity (if he already knows he's in
> quadrants I or II, it's no problem, he does an arccos() instead of an
> arctan()), he needs to correlate the sinusoid of unknown phase against
> both the reference sinusoid and the quadrature component of the
> reference sinusoid.  if you can settle the quadrant ambiguity more
> easily (say, by looking at zero crossings, but remember this is
> noisy), then fine - no need for HT.

The arctan is convenient in that it is independent of the amplitude
of the signal.

Though the place this is done most often is in the lock-in amplifier:

http://en.wikipedia.org/wiki/Lock-in_amplifier

but they usually use a square wave into the demodulator.  (It is
easier to multiply by a square wave than a sine, and it is usually
good enough.)

Many lock-in amplifiers will generate both in phase and quadrature,
some will compute the arctan for you.

They might be available for reasonable prices on eBay.

-- glen

>On Feb 6, 12:24=A0am, spop...@speedymail.org (Steve Pope) wrote:
>> robert bristow-johnson =A0<r...@audioimagination.com> wrote:
>>
>> >do you know anything about the range of phase shift between any pair?
>> >if it is between 0 and pi, then you won't need to H.T. the reference
>> >sinusoid to generate a quadrature reference sine. =A0if you have no
idea
>> >what the phase relationship is, you need both the I and the Q.
>>
>> I'm not seeing why you would ever need to Hilbert Transform
>> any signal for this problem.
>
>because of two reasons.  (if either was not the case, he wouldn't have
>to.)
>
>first, the OP is not generating the reference sinusoid.  if he was, he
>could just as well generate a synced sine that was 90 degrees outa
>phase (like in the other current thread about this same topic).
>
>second, because of quadrant ambiguity (if he already knows he's in
>quadrants I or II, it's no problem, he does an arccos() instead of an
>arctan()), he needs to correlate the sinusoid of unknown phase against
>both the reference sinusoid and the quadrature component of the
>reference sinusoid.  if you can settle the quadrant ambiguity more
>easily (say, by looking at zero crossings, but remember this is
>noisy), then fine - no need for HT.
>
>anyway, that's my spin on it.

Under noisy conditions, the arc cos approach has an accuracy that varies
hugely with the angle. You mentioned the noise issue in relation to
identifying zero crossings, but the effect of noise at particular angles is
a bigger problem. The arc tan approach offers comparable accuracy all the
way from 0 to 360 degrees.

Steve

On Feb 5, 10:21 pm, "ultrasonics" <ore...@gmail.com> wrote:
> I am developing ultrasound based imaging system for research purposes.
> I am looking for a very accurate algoritham/method for estimating a phase
> of a sinusoidal signals corrupted by noise. Frequency of the sinusoid is a
> prior information (fairly accurate), I do not care about the magnitude of a
> sinusoid. Sinusiod is excited by mechanical resonator, it will have a
> transition cycles before it reaches a steady state.  I end up with several
> sinusoids which are delayed with respect to each other. I am interested in
> estimating a phase shift between them. Sinusoids are in the frequency range
> of 25-200 Hz. Sampling frequency is lower then 5 kHz.  I was applying next
> algorithm:
>
> Choose region of interest with respect to most delayed sinusoid (to get
> rid of the transient effect). I would get complex envelope (make IQ data
> using Hilbert transform), and then use R(1) (lag one autocorrelation
> estimator) to estimate phase shift. This is a stanard Doppler processing
> (radars/ultrasound) for phae shift estimation. This works but is not that
> accurate. Main reason is that I am trying to preserve resolution in phase
> estimation, thus I am only using ensembles of 2 signals for estimating a
> phase shift( larger the ensemble and higher the SNR lower the variance of
> estimation).
>
> I would appreciate any interesting suggestion/discussion.
>
> Thanks,
>
> ultrasonics

Do you filter before the hilbert transform? If you use a simple FIR
very narrow filter around the frequency, you will get significant
noise suppression. Also, once you have phase from the hilbert, you can
generate a running average (filter(ones(1,32),1,phase_vals)/32). That
might improve the accuracy of the phase difference. Also, if you were
to perform FFT of the signals that you receive and pick the bin for
that particular frequency, calculate the phase and take the
difference. You will have to take rather short ffts if you dont have
enough points.

"parth_vakil@yahoo.com" <vakilp@gmail.com> writes:

> On Feb 5, 10:21 pm, "ultrasonics" <ore...@gmail.com> wrote:
>> I am developing ultrasound based imaging system for research purposes.
>> I am looking for a very accurate algoritham/method for estimating a phase
>> of a sinusoidal signals corrupted by noise. Frequency of the sinusoid is a
>> prior information (fairly accurate), I do not care about the magnitude of a
>> sinusoid. Sinusiod is excited by mechanical resonator, it will have a
>> transition cycles before it reaches a steady state.  I end up with several
>> sinusoids which are delayed with respect to each other. I am interested in
>> estimating a phase shift between them. Sinusoids are in the frequency range
>> of 25-200 Hz. Sampling frequency is lower then 5 kHz.  I was applying next
>> algorithm:
>>
>> Choose region of interest with respect to most delayed sinusoid (to get
>> rid of the transient effect). I would get complex envelope (make IQ data
>> using Hilbert transform), and then use R(1) (lag one autocorrelation
>> estimator) to estimate phase shift. This is a stanard Doppler processing
>> (radars/ultrasound) for phae shift estimation. This works but is not that
>> accurate. Main reason is that I am trying to preserve resolution in phase
>> estimation, thus I am only using ensembles of 2 signals for estimating a
>> phase shift( larger the ensemble and higher the SNR lower the variance of
>> estimation).
>>
>> I would appreciate any interesting suggestion/discussion.
>>
>> Thanks,
>>
>> ultrasonics
>
> Do you filter before the hilbert transform? If you use a simple FIR
> very narrow filter around the frequency, you will get significant
> noise suppression. 

That is an excellent point that I've not heard anyone mention
previously. 

If you're going to have a signal in noise, a pure sinusoid of known
frequency is the best one to have.
-- 
%  Randy Yates                  % "Midnight, on the water... 
%% Fuquay-Varina, NC            %  I saw...  the ocean's daughter." 
%%% 919-577-9882                % 'Can't Get It Out Of My Head' 
%%%% <yates@ieee.org>           % *El Dorado*, Electric Light Orchestra
http://www.digitalsignallabs.com

On Thu, 05 Feb 2009 21:17:49 -0800, robert bristow-johnson wrote:

> On Feb 5, 11:09&nbsp;pm, Tim Wescott <t...@seemywebsite.com> wrote:
>> On Thu, 05 Feb 2009 21:21:30 -0600, ultrasonics wrote:
>> > I am developing ultrasound based imaging system for research
>> > purposes. I am looking for a very accurate algoritham/method for
>> > estimating a phase of a sinusoidal signals corrupted by noise.
>> > Frequency of the sinusoid is a prior information (fairly accurate), I
>> > do not care about the magnitude of a sinusoid. Sinusiod is excited by
>> > mechanical resonator, it will have a transition cycles before it
>> > reaches a steady state. &nbsp;I end up with several sinusoids which are
>> > delayed with respect to each other.
> 
> do you know anything about the range of phase shift between any pair? if
> it is between 0 and pi, then you won't need to H.T. the reference
> sinusoid to generate a quadrature reference sine.  if you have no idea
> what the phase relationship is, you need both the I and the Q.
> 
>> > I am
>> > interested in estimating a phase shift between them. Sinusoids are in
>> > the frequency range of 25-200 Hz. Sampling frequency is lower than 5
>> > kHz. &nbsp;I was applying next algorithm:
>>
>> > Choose region of interest with respect to most delayed sinusoid (to
>> > get rid of the transient effect).
> 
> in other words, your region of interest is the most current data. this
> is what you get because the LP filters used for the correlation Tim
> refers to tend to forget about input that happened longer ago.
> 
>> > I would get complex envelope (make IQ data using Hilbert transform),
>> > and then use R(1) (lag one autocorrelation estimator) to estimate
>> > phase shift. This is a standard Doppler processing
>> > (radars/ultrasound) for phase shift estimation. This works but is not
>> > that accurate.
> 
> we were talking about that in the other thread also.
> 
>> > Main reason is that I am trying to preserve resolution in phase
>> > estimation, thus I am only using ensembles of 2 signals for
>> > estimating a phase shift( larger the ensemble and higher the SNR
>> > lower the variance of estimation).
> ...
>> Since you already know the frequency, why not correlate the sine wave
>> with its I & Q prototypes over some integer number of cycles, then find
>> the polar angle of the resulting complex number?
> 
> isn't this what we were just talking about in this thread?:
> 
>  http://groups.google.com/group/comp.dsp/msg/0008e4899b8d85ee
> 
> heck, i even put a "noise" term in it.
> 
> depending on what kinda math library the OP has to work with, there may
> be some effort needed to perform the arctan() function.
> 
> r b-j

So it is.

Yes, the arctan function can be a barrier, but there are many ways to 
implement them.

-- 
http://www.wescottdesign.com