Delay estimation by FFT

Posted by Markus Nentwig on Sep 22 2007 under Matlab | Basics | Tips and Tricks   

Summary: The delay between a signal and a reference appears in the frequency domain as linear phase shift. It can be estimated with subsample accuracy using FFT methods. This article provides a Matlab / Octave example program and discussion.

Hello,

there were several recent posts on comp.dsp that involved the following problem:

"The transmitted signal is a pulse, for example five cycles taken from a sine wave. How can I find the pulse location from the received signal?"

It's a textbook problem (matched filter, correlation receiver) and I do hope that it's not a homework problem...

Still, there was quite some discussion. An interesting idea was brought up:
"The delay between transmitted and received pulse appears as linear phase shift in the frequency domain. Therefore, I could take the FFT of received signal and original pulse, calculate the phase difference over frequency, and obtain the delay from the slope of the phase."
One advantage is that the delay can be resolved with subsample accuracy without need for excessive amounts of oversampling.
The problem is that the phase wraps around and is not that straightforward to unroll. But of course it can be done. The given example program is more robust than standard "phase unwrapping" and may work in cases, where the latter fails.

The approach may not be the first choice for hardware implementation (and neither for "matched filtering" homework problems). But it can be quite useful when it comes to for example offline processing of measured data.

Note that FFT processing as used in the example treats signals as cyclic. Usually this is not a restriction, simply use as much zero-padding as appropriate.

The picture below shows an example: the red curve is the transmitted pulse. Blue is the received signal, delayed by 10.45 samples. The example uses a rather short delay for clarity, but it works just as well with longer delay, of course at the expense of more memory and CPU cycles for the FFTs.

The "ringing" before and after the time-shifted pulse remains invisible in the original signal, because it is sampled exactly on the zero-crossings of the waveform (Nyquist's pulse shaping rule for an ideal lowpass filter).

The crosscorrelation is obtained by taking:

ifft(fft(signal) .* conj(fft(reference)))

Below the result for the windowed sin()-pulse:

The windowed sin()-pulse has rather "bad" autocorrelation properties. It is "too similar too itself" at some offsets, resulting in several crosscorrelation peaks. Nonetheless, the highest pulse corresponds to the delay between signal and reference.

For comparison, here is the same result using a different pulse consisting entirely of random numbers: When it is time shifted, it has very little similarity to itself, and there is only one clearly defined peak:

Ideally, it it wouldn't be necessary to search for the peak, because the delay can be determined from the phase of the cross spectral density fft(signal) .* conj(fft(reference)) alone.

The absolute value of above expression is related to the energy (amplitude squared) per bin, of both reference and received signal.
The phase is "the phase of the received signal, rotated back by the phase of the reference", in other words, the phase shift at each frequency caused by the delay.

The picture below shows said phase when using the windowed sin()-pulse: The leftmost bin corresponds to 0 Hz, the center bin to half the sampling frequency. The right half of the picture represents negative frequencies:

By eye, it is easy to see that the phase goes through about ten-and-a-half full turns (exactly: 10.45).
But there are many "spikes", caused by nulls in the spectrum of the pulse. Consequently, it is difficult to automatically unwrap the phase in a general way that is robust in the presence of noise.
Now the delay is already known from the location of the crosscorrelation peak, at least with an accuracy of +/- 0.5 samples.
Therefore, the next step is to rotate the phase backwards to compensate for the known delay (with one sample resolution), leaving only the phase corresponding to a subsample delay. The remaining phase is shown below:

The slope can be determined by using a least-squares fit to the phase of a unity delay, as shown in red.

This works, but is not very accurate. The reason is that most of the bins in the cross-spectrum contain little or no energy. In particular, the windowed sin()-pulse has a rather narrow spectrum.
The solution is to weight the phase values on "both sides" of the least-mean-squares-solution with the energy contained in the respective bin. For the windowed sin-pulse, the result looks as follows:

The location of the peaks correspond to the frequency of the windowed sinewave (positive and negative frequencies). Said pulse is a "difficult" case, because the energy is concentrated in a narrow frequency region.

As it turns out, the weighting improves the accuracy by about eight orders of magnitude! For a noise-free signal, a typical error is below 10^-10 the sampling time interval, which could not be easily achieved by simple oversampling.

The program can be found here:
Download (commented Matlab / Octave program)
The relevant code excluding comments and test signal generation consists of the 20 lines below:

% ref: transmitted signal % sig: received signal n=length(sig); binFreq=(mod(((0:n-1)+floor(n/2)), n)-floor(n/2))/n; sig_FD=fft(sig); ref_FD=fft(ref, n); u=sig_FD .* conj(ref_FD); if mod(n, 2) == 0 u(length(u)/2+1)=0; end Xcor=abs(ifft(u)); integerDelay=find(Xcor==max(Xcor)); integerDelay=integerDelay(1)-1; rotN=exp(2i*pi*integerDelay .* binFreq); uDelayPhase=-2*pi*binFreq; u=u .* rotN; weight=abs(u); uDelayPhase=uDelayPhase.* weight; ang=arg(u).* weight; fractionalDelay=transpose(uDelayPhase)\transpose(ang); result=integerDelay+fractionalDelay;

Comments, additional references, bug reports are of course welcome!

Cheers

Markus

Update 16.8.2011:
An extended ready-to-use implementation for Matlab can be found in the code snippet.

References:

Similar approach, but uses FFT and oversampling to locate the crosscorrelation peak: http://dspace.dsto.defence.gov.au/dspace/bitstream/1947/4441/1/DSTO-TR-1705%20PR.pdf

Methods that could be used to locate the peak without oversampling, when it falls between bins:
http://www.dspguru.com/howto/tech/peakfft.htm

A similar approach is discussed in:
Bendat & Piersol: "Random Data", Wiley, 20090.
Thanks to Rune Allnor for this reference!
 

 

Rate this article:
posted by Markus Nentwig
Markus received his Dipl. Ing. degree in electrical engineering / communications in 1999. Work interests include RF transceiver system design, implementation and modeling. He works as senior architect for Renesas Mobile Europe in Finland.

Previous post by Markus Nentwig: Polyphase filter / Farrows interpolation
Next post by Markus Nentwig: A brief look at multipath radio channels
all articles by Markus Nentwig

Would you like to be notified by email when Markus Nentwig publishes a new blog?

  

Hello Markus,
i was wondering if something similar to your program could be used in an ultrasonic anemometer to resolve the wave-length ambiguity in determining the time between the emitted ultrasound wave and the received signal.
By measuring the phase shift of the received signal from the zero crossing samples, I get quite accurate values within one wave length but there is no obvious way to determine the number of pulses of the delay. Due to the inertia of the mechanical transceiver, there is no sharp rise in amplitude...

Hello,
yes, that should work (depending on the pulse and its autocorrelation properties).
Can you provide some sample data for transmitted and received signal? It might be a nice example.

Hello Markus

the autocorrelation is where the difficulty lies.
Whatever modulation one tries to impress on the sender signal (phase reversal or starting and stopping after a short pulse train), the reaction in the receiver is always delayed by the mechanical inertia. If there is more cross-wind, then less energy would be transferred and the delay is longer.
Here are two pictures of the waves seen on an oscilloscope:
http://www.flashgenie.net/anemometer50uS.gif

http://www.flashgenie.net/41kHzPulse.gif
(the upper wave is the receiver in each picture, the lower one the sender, which has some ringing.)
The best method would appear to be to renounce directly determining the number of pulses that have elapsed and instead use your algorithm to make a more accurate measurement of the exact phase shift of the received signal, and to repeat the measurement at a slightly different frequency. The difference between the two phase shift values would then allow one to resolve the ambiguity of the number of pulses.
If one used 39 kHz and 41 kHz, for example, one could resolve the ambiguity within 20 pulses, which would be enough for an anemometer where the transducers are not too far apart.
In the present experimental circuit I was only measuring the zero-crossing times of the received signal, because the microprocessor can do that to within 10 ns, but there is too much jitter in the received signal to obtain sufficiently accurate calculations.
If one used an AD-converter one could measure ca 25 analog values on each pulse and use your calculation method to find the best fitting phase-shift, right ?

Hello,

sorry, it took me a while, busy week...

I think it would work just fine, but
it's maybe overkill for the given problem.

If your hardware is able to transmit and receive a 39 kHz / 41 kHz two-tone signal, the common cycle length is 1 ms.
For example, I could multiply the received signal with two complex-valued NCOs at -39 / - 41 kHz, and integrate both (complex) over 1 ms. Afterwards, it should be possible to map between the phase angles of both complex sums and the delay.
How to implement that mapping is another question, but I'm sure it will require much less CPU cycles than the FFTs.
I have a C example program on my disk that could be used to put together a prototype with a soundcard and microphone (of course at audio frequencies). Please let me know, if interested.

Cheers

Markus

Hello Markus

yes, i would definitely be very grateful to try your idea to get me out of the present dead-lock!
My mathematical knowledge in this area is very limited... To be able to follow you, I would need to understand a little better what this is about -- could you point me to any literature ?
Right now I am using a pwm signal to alternately send 39 kHz and 41 kHz pulse trains. To reduce ringing in the rest of the circuit, I send the signal through a simple sine-filter. The ultrasonic transducers have a rather narrow resonant frequency, so it is possible there is some beat between the pwm frequency and their own resonant frequency.

Hi,
here is a program that you could use as basis for a prototype with a soundcard, uses portaudio library:

http://www.elisanet.fi/mnentwig/webroot/looprecord/index.html

I found also another program, I wrote this a while ago to measure the acoustic path length with an arbitrary test signal and FFT-based crosscorrelation (requires fftw3).

It uses simple oversampling, not the more complicated (and accurate) method in the blog article.

http://www.elisanet.fi/mnentwig/dspblog/fft_align/portaudio_example.c

many belated thanks!

(it would be nice if there were some way to have an option of receiving an email update when there is a new entry in the log)

Hello Markus,
I want to get a frequency independent phase shift for my baseband signal. This phase shift should be something like 20,30,40 and so on degrees. I was thinking of using a FFT on the baseband and then altering the phase information somehow without changing the magnitude and then taking the IFFT to get a phase shifted signal.
The goal being that each frequnecy component of the input should be shifted by the same amount of degrees.
Do you know how can I change the phase information of the FFT signal , so that when I take an IFFT it results in the correct ouput.
Or if you have some other suggestions please do let me know

Thanks

>> I want to get a frequency independent phase shift for my baseband signal

Right, you can take FFT and multiply each frequency component individually with exp(i angle). At least it's going to work as long as there is a full cycle length of the (periodic!) signal inside the FFT.

An alternative solution in time domain for nonperiodic signals could use a Hilbert transform filter, which gives 90 degree phase shift for bandpass signals.
As the lower corner frequency approaches 0 Hz, the filter length grows to infinity, though.

It might be a good idea to ask the question on comp.dsp.

hi markus .

I read two waveforms (acoustic) picked up by mics as a .wav file in matlab using wavread command. One of these is reference (ref in ur code) and other is signal..however running ur code with this shows the following error :-

"??? Error using ==> run
Error using ==> times
Matrix dimensions must agree."

I cant decipher the error as of now..

please guide..

regards

dks

hi markus ,

the problem seems to be coming with
the following line of code in ur code :-

" u=u .* rotN;"
I am still trying to figure out.

dks

Hello,

you could compare size(u) and size(rotN).
Probably wavread gives a column vector, but it should be a row vector. If so, simply transpose the data and it should be fine.
Also, check whether the file is one-channel audio..

hi markus..

I did transpose my input data and the system as of now works ..i am yet to try it on field data to check the accuracy of delay estimation..The field data through a software has already estimated time delay between various received waveforms..As of now , I have to manually clip the waveform from the software ( Adobe Audition ) and then read it in matlab..do u know of any means to directly read it from the software.. Ur guidance has been of great help in reaching the final goal..thanks always ..n regards..dks

hi markus

i have sent 2 mails at ur website mail address and have been anticipating ur reply..I have sent these mails from my mail id : rashdk@iitk.ac.in.
In trying to compile ur above program (the advanced version from ur website) , i get an error in compilation for as :
Error using ==> vertcat
CAT arguments dimensions are not consistent.

Error in ==> delay at 78
r=transpose([constRotPhase; uDelayPhase])\transpose(ang);

i have tried resizing but to no avail..
could u suggest ways to oversome this ..

rashdk

hello markus
i tried ur progrm in my set up of 2 mics (placed 3 m each other ) and recorded a gunshot at distance of 60 m. The cross corr peak and hence the time delay that i get is in samples ?? so if i divide that by sampling freq , shouldn't i get the time delay of arrival at farther mic ...the progrma as of now yields this in seconds !!!! can u suggest me the errors..

Hi,

I want to do the delay shift in the time domain. But obviously I cannot shift 0.5 samples in time domain. So I thougth to do oversampling. But to perform the fft I need to down sample again and the shift is either 0 or 1.

Does anyone have an idea?

Thanks!
oliht

Hi,

please have a look at this other page: http://www.dsprelated.com/showarticle/22.php

One straighforward solution is to design an FIR filter that has a delay offset of 1/2 sample. The filter will do some "damage" to the signal, that can't be avoided (but reduced by using a longer filter).

Compared to the "polyphase" filter in the above page, you're only interested in a "single" phase. So you'd use only the "red" (or blue or green, those that match your 1/2 sample delay offset) filter coefficients.

Hello Markus,

Why in the program fft_align.m for definition integerDelay the maximum of correlation function is used, instead of a phase slope in frequency domain?

Thanks!

Hi Markus
If there is a Frequency shift,and the sequence is a Zadoffchu_sequence ,Where the peak will be located ,And i think the peak position is dependent on (U - logicalrootsequence) which is used to generate the zadoff_chu sequence
please help in undestanding this
this is the psuedocode

U = 140;
for n= 0:838

conv_result(n+1) = exp((-j*pi*U*n*n+1)/839));
shift_result(n+1)= conv_result(n+1)*exp((j*2*pi*6*n)/839);
end
[corr lag] = xcorr(shift_result,conv_result);

Is ther any thing iam missing here iam not able to mathematically locate the peakposition .

thanks

Hi Markus
Did u know about reconstruction of image using Fourier transform?in matlab..