Time varying phase correction

Hello all,

I have a bunch of datasets that I'm trying to phase align so I can
average them.  These data sets are interferograms, consisting of a
carrier frequency at 50 kHz.  They're sampled at ~1 MHz and are ~1.5
million samples long.

My method currently is to start with the first data set (D1) and phase
align the second data set (D2) to it by getting the phase delay from
the auto-correlation.  I correct this phase offset by something like
real(ifft(fft(D2).*exp(2*pi*f*delay/Fs)))  where Fs is the sampling
frequency.   I then average D1 and D2_phaseshifted and move on to the
next data set and phase align it to this average and average them
cumulatively.

This works great, however it only corrects the average phase
difference.  It seems my phase difference drifts by +/- 1 sample over
the whole of the data set.  This means the data sets are usually in
phase in the center of the interferogram but can be out of phase by 20
degrees near the ends.  This leads to a distorted average (think cigar
shaped, fat in the middle tapering towards the ends).

This phase variation is to be expected because the interferometer will
drift from acoustic/thermal vibrations.

I've been looking for a method that would allow me to correct this
time varying phase difference.  So far I've come up with a method that
bandpasses the interferograms about the carrier, checks for zero
crossings, calculates the number of samples between zero crossings of
two data sets then does a moving average of that.  Unfortunately this
is slow and kludgy.

This strikes me as something that would be commonplace, so does anyone
know of a more elegant or at the very least faster way of aligning
this kind of data and avoiding the distortion I am talking about.

Ideally, I'd like to know for each sample of the interferogram the
number of samples to phase shift it to align it with the reference
interferogram.  And I'd like to do this quickly and with as little
babysitting as possible.

(I do realize that taking the first data set as a reference isn't a
good idea either since it has as much noise as the other data sets, so
now I am thinking I should create an ideal phase reference based on
the average of the first two sets or something...)

Any insight would be greatly appreciated.

-Frank

Frank McDermott wrote:

> Hello all,
> 
> I have a bunch of datasets that I'm trying to phase align so I can
> average them.  These data sets are interferograms, consisting of a
> carrier frequency at 50 kHz.  They're sampled at ~1 MHz and are ~1.5
> million samples long.
> 
> My method currently is to start with the first data set (D1) and phase
> align the second data set (D2) to it by getting the phase delay from
> the auto-correlation.  I correct this phase offset by something like
> real(ifft(fft(D2).*exp(2*pi*f*delay/Fs)))  where Fs is the sampling
> frequency.   I then average D1 and D2_phaseshifted and move on to the
> next data set and phase align it to this average and average them
> cumulatively.
> 
> This works great, however it only corrects the average phase
> difference.  It seems my phase difference drifts by +/- 1 sample over
> the whole of the data set.  This means the data sets are usually in
> phase in the center of the interferogram but can be out of phase by 20
> degrees near the ends.  This leads to a distorted average (think cigar
> shaped, fat in the middle tapering towards the ends).
> 
> This phase variation is to be expected because the interferometer will
> drift from acoustic/thermal vibrations.
> 
> I've been looking for a method that would allow me to correct this
> time varying phase difference.  So far I've come up with a method that
> bandpasses the interferograms about the carrier, checks for zero
> crossings, calculates the number of samples between zero crossings of
> two data sets then does a moving average of that.  Unfortunately this
> is slow and kludgy.
> 
> This strikes me as something that would be commonplace, so does anyone
> know of a more elegant or at the very least faster way of aligning
> this kind of data and avoiding the distortion I am talking about.
> 
> Ideally, I'd like to know for each sample of the interferogram the
> number of samples to phase shift it to align it with the reference
> interferogram.  And I'd like to do this quickly and with as little
> babysitting as possible.
> 
> (I do realize that taking the first data set as a reference isn't a
> good idea either since it has as much noise as the other data sets, so
> now I am thinking I should create an ideal phase reference based on
> the average of the first two sets or something...)

Break your signal in two halves, correlate each half and find the timing 
offsets for both, find the difference between the offsets, calculate the 
clock drift from there, resample the signal so to compensate for the 
clock drift and the timing offset.

> Any insight would be greatly appreciated.

How much is the great appreciation?



Vladimir Vassilevsky
DSP and Mixed Signal Design Consultant
http://www.abvolt.com

Hello,

If resampling on a uniform grid with +/- n samples will do (I doubt it):
- fft
- zero-pad or remove n samples from the center bins
- ifft. Depending on n, this may be quick or slow.
- you can zero-pad the input signal to increase resolution. For example,
10Msamples should be still manageable.

A more general solution: 
1) estimate the delay at multiple points. 
2) use (linear) interpolation between the delay values, and get a
resampling time instant for each individual sample.
3) resample the your original signal (FFT) for the time instant of each
sample. 

For 1):
apply some windowing function to cut a piece out of signal and reference,
short enough that drift is negligible during its length.
[1] estimates the time offset.

For 2): matlab interp1()

For 3) quick-and-dirty: http://www.dsprelated.com/showcode/206.php
this will work for small examples, but not with the full dataset size. 
Try Lagrange interpolation instead, for example
http://www.dsprelated.com/showcode/208.php
"outTime" would be the target sampling instants.

If it needs to be more accurate and / or it has strong high frequency
content, a polyphase FIR filter should be the solution.
If performance matters, the whole thing could be combined into some kind of
overlay-add algorithm
- extract window
- fft
- estimate local delay
- use fft time shift
- merge into output.

And / or, maybe your zero crossings give already enough information if you
lowpass-filter the result to a very narrow bandwidth, let's say 10 cycles
over the duration of your data set. Then use this as delay estimate,
construct the resampling grid and interpolate.

Vladimir Vassilevsky <nospam@nowhere.com> wrote:

>> I have a bunch of datasets that I'm trying to phase align so I can
>> average them.  These data sets are interferograms, consisting of a
>> carrier frequency at 50 kHz.  They're sampled at ~1 MHz and are ~1.5
>> million samples long.

(snip)
>> This works great, however it only corrects the average phase
>> difference.  It seems my phase difference drifts by +/- 1 sample over
>> the whole of the data set.  

(snip)
> Break your signal in two halves, correlate each half and find the timing 
> offsets for both, find the difference between the offsets, calculate the 
> clock drift from there, resample the signal so to compensate for the 
> clock drift and the timing offset.

I was reading the OP and figuring that the correction should be
a constant term plus a linear (drift) term, but hadn't thought
how you determine the drift term.  Then I hit next.

That does sound about right, though.  The average of the two gives
the constant term, the difference the linear term.  But if that
isn't enough, do you go to a quadratic term?  Break into thirds,
then find the quadratic that goes through the three offsets?

-- glen

this link was missing in my earlier post
[1] http://www.dsprelated.com/showcode/207.php
accurately estimates the delay between two signals.

On Wed, 24 Aug 2011 20:44:12 -0500, Vladimir Vassilevsky wrote:

> Frank McDermott wrote:
> 
>> Hello all,
>> 
>> I have a bunch of datasets that I'm trying to phase align so I can
>> average them.  These data sets are interferograms, consisting of a
>> carrier frequency at 50 kHz.  They're sampled at ~1 MHz and are ~1.5
>> million samples long.
>> 
>> My method currently is to start with the first data set (D1) and phase
>> align the second data set (D2) to it by getting the phase delay from
>> the auto-correlation.  I correct this phase offset by something like
>> real(ifft(fft(D2).*exp(2*pi*f*delay/Fs)))  where Fs is the sampling
>> frequency.   I then average D1 and D2_phaseshifted and move on to the
>> next data set and phase align it to this average and average them
>> cumulatively.
>> 
>> This works great, however it only corrects the average phase
>> difference.  It seems my phase difference drifts by +/- 1 sample over
>> the whole of the data set.  This means the data sets are usually in
>> phase in the center of the interferogram but can be out of phase by 20
>> degrees near the ends.  This leads to a distorted average (think cigar
>> shaped, fat in the middle tapering towards the ends).
>> 
>> This phase variation is to be expected because the interferometer will
>> drift from acoustic/thermal vibrations.
>> 
>> I've been looking for a method that would allow me to correct this time
>> varying phase difference.  So far I've come up with a method that
>> bandpasses the interferograms about the carrier, checks for zero
>> crossings, calculates the number of samples between zero crossings of
>> two data sets then does a moving average of that.  Unfortunately this
>> is slow and kludgy.
>> 
>> This strikes me as something that would be commonplace, so does anyone
>> know of a more elegant or at the very least faster way of aligning this
>> kind of data and avoiding the distortion I am talking about.
>> 
>> Ideally, I'd like to know for each sample of the interferogram the
>> number of samples to phase shift it to align it with the reference
>> interferogram.  And I'd like to do this quickly and with as little
>> babysitting as possible.
>> 
>> (I do realize that taking the first data set as a reference isn't a
>> good idea either since it has as much noise as the other data sets, so
>> now I am thinking I should create an ideal phase reference based on the
>> average of the first two sets or something...)
> 
> Break your signal in two halves, correlate each half and find the timing
> offsets for both, find the difference between the offsets, calculate the
> clock drift from there, resample the signal so to compensate for the
> clock drift and the timing offset.

Or, if the carrier signal is strong enough, do all of them to a reference 
rather than doing n-1 of them to the first one you happen to pull out of 
a bag.

I'm not sure if it'll make much of a difference, other than perhaps being 
just a hair less noisy in an already fairly noise-proof set of data.

-- 
Tim Wescott
Control system and signal processing consulting
www.wescottdesign.com

Tim <tim@seemywebsite.please> wrote:

(snip)
>>> I have a bunch of datasets that I'm trying to phase align so I can
>>> average them.  These data sets are interferograms, consisting of a
>>> carrier frequency at 50 kHz.  They're sampled at ~1 MHz and are ~1.5
>>> million samples long.

(snip)
> Or, if the carrier signal is strong enough, do all of them to a reference 
> rather than doing n-1 of them to the first one you happen to pull out of 
> a bag.

Well, you could do all pairs, or some number of pairs, and then
combine them such that the effect isn't so far off.

Maybe if you have constant and linear shift terms from all
relative to the first, then you can combine them in some way
to reduce noise and other effects from the choice of a first one.

Or do consecutive pairs (1,2), (2,3), (3,4), ...  (n,1)
that removes any special 'first' data set.

> I'm not sure if it'll make much of a difference, other than perhaps being 
> just a hair less noisy in an already fairly noise-proof set of data.

-- glen

Thanks for the suggestions guys.

I feel pretty dumb for not thinking of Vladimir's suggestion. I probably thought doing multiple cross-correlations would slow things down immensely since it seems to be the bottleneck with regards to computation currently.  

Glen is right that I think I'd have to do more than just split the data set since it seems for some worst case scenarios the data can shift +1, -1, and then back to +1 over the course of a data set.  Out of 100 datasets, the worst case scenario was a 5 sample drift.  However if you overlay all datasets you can see the drift is about a definitive mean which makes me think Tim's idea of referencing to an ideal carrier is a good one (did you suggest this idea in another post? I think I read it somewhere on comp.dsp and thought it sounded like a good idea).

I am wondering about re-sampling now.  This touches on the questions of accuracy that mnentwig brought up.  The datasets are actually sampled equally with respect to optical frequency (the clock signal that drives the ADC is generated from a reference interferometer).  The data sets are FFT'd, window'ed about the carrier, and then iFFT'd to extract certain parameters.  Would all resampling have to occur at integer multiples of the sampling frequency?  What sort of noise addition would I expect from the interpolation? 

-Frank

Hi,

I had a look, and unfortunately my earlier proposal to estimate the timing
error does not work. It needs to be done differently, but shouldn't be
impossible. With the time-varying estimate, resampling should be easy and
fast. 

Now getting the time from the zero crossings doesn't sound like a bad idea.
It should be quite accurate, if you fit a curve to it. And fast, too.

If the frequencies are otherwise stable (50k and 1M = 20 * 50k), this would
fit a least-squares sine- and cosine wave across a zero crossing location:

    sBase = sin((0:19)/20 * 2 * pi);
    cBase = cos((0:19)/20 * 2 * pi);
    
    p = pinv([sBase.' cBase .'])
    c = p * sigNoisy(ix:ix+19).'

Then get the phase from the elements in c via atan2().
If it's too sensitive to frequency error, some polynomial fit to a few
nearby samples might be an alternative.

.. or multiply directly with sine / cosine ...