# Time Delay Estimation Techniques - Part 1

### Introduction

In this article we will look at signal processing techniques for time delay estimation.

### Background

Time delay estimation has been a research topic of significant practical importance in many fields like radar, sonar, seismology, geophysics, ultrasonic’s, hands-free communications,Doppler positioning systems etc and is of fundamental importance in a variety of signal-processing applications.

The problem boils down to computing the location of a digitized “signature” waveform residing within a larger time-slice.

Let $s(n)$ be the sound source signal. Let us consider two spatially separated sensors and $x_{1}(t)$ and $x_{2}(t)$ be result of propagation of $s(t)$ through different paths to reach the respective sensors.

A simple propagation model is that the signals just encounter attenuation and delay and are corrupted by additive noise

$x_{i}(t) = a_{i} s(t- \tau_{i}) + b_{i}(t)$

where ,

$b_{i}(t)$ is the additive noise uncorrelated with the signal

We make a assumption that the distance from the sensors to the source is very large compared to the spatial separation between the sensors.Thus signal $s(t)$ received by both the receivers is the same,there is no significant change in the signal as it travels thorough different paths to reach the spatially separated sensors.

### Cross-correlation

One of the simplest method of time delay estimation is cross-correlation . The cross-correlation of two signals is a measure of similarity between the two sequences. The cross-correlation function is maximized when both the signals have significant overlap.

$R_{xy}(k) = \sum_{i} {x(i) y(i+k)} = x[n]*y[-n]$

Compute maximum absolute value of the correlation function to estimate the lag.

To test the results we create create two sequences,one a delayed version of another. We add white noise to the delayed sequence and use sample correlation to detect the lag. while performing correlation we normalize the signals,so that correlation measure is bounded between [0,1]

Now we increase the noise to try to get estimate of noise co-variance at which this technique fails. we run the experiment 100 times and estimate the number of times we get the proper answer.

#### Testing with Different Additive Noise Covariance

 *********** Information ************
time delay :  10
Noise : 0.1
SNR 15.4390474505
Correct  98  Incorrect  3
mean  9.90099009901 Std  1.0000490136


We see that we are always correct,when the noise levels are low or SNR is high

 *********** Information ************
time delay :  10
Noise : 0.3
SNR 5.89662235611
Correct  48  Incorrect  53
mean  9.79207920792 Std  1.2767287382


The SNR has dropped to about 6db and error in estimation is about 50%. The accuracy is linearly related with the SNR

 *********** Information ************
time delay :  10
Noise : 0.5
SNR 1.45964736379
Correct  36  Incorrect  65
mean  9.80198019802 Std  1.57969668714



At 2 db SNR we can see the variance of estimated time delay increasing with rising noise levels.

We can see in the auto-correlation plot the difference between the peak and its neighbors is not significant and depending on the random noise levels introduced,we may not always get the right answer.

As noise increases,we can see variance in the estimated time delay increases and error in estimation also increases.As the level of noise increases, the uncertainty in the time-delay estimate increases

We need a reliable time delay estimator in presence of additive noise.

If is always difficult to estimate the time delays for a base-band signal illustrated above.The due to the noise,we are not reliably able to estimate the delay corresponding to the maximal overlap between the noisy and ideal signal.

def delay(signal,N):
""" function introduces a delay of N samples

Parameters
-----------
signal : numpy-array,
The input signal

N      : integer
delay

Returns
--------
out : numpy-array
delayed signal

"""
d=numpy.zeros((1,N+1));
signal=numpy.append(d,signal)
return signal;

Parameters
-----------
s : numpy-array,
The input signal

N      : float
noise covariance

Returns
--------
out : numpy-array
noisy signal

"""
noise = np.random.normal(0,variance,len(s))
s=s+noise;
return s;

if __name__ == "__main__":

Fs=1000;

mode=0
if mode==0:
x = triang(20);
x2=x;

......

tdelay=10;
varnoise=0.001;
loop=1000

#delay the signal
dx=delay(x,tdelay);
result=numpy.zeros((1,loop));
for i in range(loop):
s=dx;

#normalize the signals
s=s/np.linalg.norm(s);
x1=x2/np.linalg.norm(x2);

unfiltered_signal=s;

r=numpy.correlate(s,x1,mode="full")
arg=np.argmax(r)
result[0,i]=abs(arg-len(x))

#plot the results
c=np.sum(result==tdelay)
ic=np.sum(result!=tdelay)
#10*np.log(np.sum(np.abs(x*x))/
print " *********** Information ************ "
print 'time delay : ',tdelay
print 'Noise :',varnoise
print "SNR",10*np.log(np.mean(abs(x*x))/(varnoise*varnoise))/np.log(10);
print "Correct ",str(c)," Incorrect ",str(ic)
print "mean ",np.mean(result),"Std ",np.std(result)
print result
plt.figure(1)
subplot(2,2,1)
plt.plot(range(len(x)),x)
xlabel('Time')
ylabel('Amplitude')

subplot(2,2,2)
plt.plot(range(len(unfiltered_signal)),unfiltered_signal)
xlabel('Time')
ylabel('Amplitude')

subplot(2,2,3)
plt.plot(range(len(r)),r)
xlabel('Time')
ylabel('Amplitude')


### Rectangular Pulse signals

Let us look at the results for a different signal in the form of a rectangular pulse . we will consider the wave of same duration 30.

 *********** Information ************
time delay :  10
Noise : 0.5
SNR 0.791812460476
Correct  86  Incorrect  15
mean  9.92079207921 Std  1.11411438271


we see that signal has a lower SNR,but accuracy and estimated time delay variance is lower.

Thus the type of pulse we use has a impact on the accuracy of auto-correlation function.

def rectangular(N,start,end):
""" fuction generates a rectangular pulse

Parameters
---------
N : integer
length of signal

start,end: integer
starting and ending index of pulse

"""
x = np.zeros((1,N))
x[:,start:end]=1;
x=x.flatten();
return x;


### Coded Pulses

Broadband techniques have a sequence of code pulses that increase the accuracy of time delay estimation in the presence of noise.

time delay :  10
Noise : 0.5
SNR 2.55272505103
Correct  95  Incorrect  6
mean  9.93069306931 Std  1.01725144864



x=rectangular(20,8,14)+rectangular(20,2,5)


But most of the time we do not have control of the source pulse.

In the remainder of the article we will assume that it is a rectangular pulse of duration 1000 with impulse of duration 50 starting at 100.

time delay :  10
Noise : 0.5
SNR -6.98970004336
Correct  856  Incorrect  144
mean  9.991 Std  0.561176442841


x=rectangular(1000,100,150)


Now given source signal,we need to see if we can do better in the presence of additive noise atleast.In real life situations there will be a host of other distortions and effects which will increase the estimation errors apart from noise.

For pulses like the ones observed above,noise is a dominant factor,signal energy is low compared to noise energy.

### Modulated Pulses

Often the pulses are modulated by sinusoidal waves for longer range transmissions. </pre> time delay : 10 Noise : 0.5 SNR -6.98970004336 Correct 988 Incorrect 12 mean 10.004 Std 0.109471457467 Actual time delay 10 </pre>

def sinepulse(N,start,end,f,Fs=1000,tau=0):
""" function generates modulated rectangular pulse

Parameters
---------
N : integer
length of signal

start,end: integer
starting and ending index of pulse

f   : integer
modulated carrier freuency

Fs  : integer
Sampling freuency

tau : integer
carrier phase delay in samples.

Returns
--------
out : numpy-array
modulated rectangular  signal

"""
x = np.zeros((1,N))
t=np.asarray(range(0,end-start));
x[:,start:end]=np.cos(2*np.pi*f*(t+tau)/Fs);

return x.flatten();


### Carrier Synchronization Issues

we can synchronize exactly with the carrier ,then like broadband techniques we can achieve a significantly enhanced SNR. However synchronizations are never possible.

we introduce a phase delay of 1 sample to check the effect of carrier phase errors

*********** Information ************
time delay :  10
Noise : 0.3
SNR -2.55272505103
Correct  2  Incorrect  998
mean  229.036 Std  244.973840857



In such cases another approach might be to use rectangular envelope but in the present case that also does not seem to help

 *********** Information ************
time delay :  10
Noise : 0.3
SNR -2.55272505103
Correct  0  Incorrect  1000
mean  325.336 Std  272.058543523



#### Envelope Detection

We perform envelope detection on the signal and then apply correlation. This improves the situation considerably

 *********** Information ************
time delay :  10
Noise : 0.3
SNR -2.55272505103
Correct  661  Incorrect  339
mean  9.746 Std  0.793400277288


            if mode==5 or mode ==6 or mode==7:
s=abs(scipy.signal.hilbert(s))


### Filtering

If we known the carrier frequency ,we can filter the noise outside the signal bandwidth before performing the correlation operation.

For the remaining examples we consider signal with carrier frequency of 1KHz and Sampling rate of 5Khz. The reason for that is explained below

Let us first look at the frequency characteristics of the rectangular pulse

The sharp transition at the edges of the rectangular pulse give rise to high frequency components. Thus while filtering we need to take these into account,else we may end up distorting the edge leading to errors in time delay estimation.

The width of the mainlobe in frequency domain is inversely proportional to width of rectangular pulse in the time domain.Thus a wider pulse in the time domain will provide higher spectral compression in the frequency domain.

we can see that main lobe of the rectangular pulse would have normalize freuency of 0.05 Few of the side lobes also contain significant information.


def plotFrequency(b,a=1):
""" the function plots the frequency and phase response """
w,h = signal.freqz(b,a)
h_dB = abs(h);#20 * np.log(abs(h))/np.log(10)
subplot(211)
plot(w/max(w),h_dB)
#plt.ylim(-150, 5)
ylabel('Magnitude (db)')
xlabel(r'Normalized Frequency (x$\pi$rad/sample)')
title(r'Frequency response')
subplot(212)
h_Phase = np.unwrap(np.arctan2(np.imag(h),np.real(h)))
plot(w/max(w),h_Phase)
xlabel(r'Normalized Frequency (x$\pi$rad/sample)')
title(r'Phase response')


We try with a normalized frequency components of 0.05 which contains around 3 adjacent side-lobes

we can see that the signal has been attenuated near the edges.Thus the distortion in the region where the pulse starts will lead to ambiguity of time delay computation.

Thus we need to consider the freuency components to retain so that sharp transition in time domain are not affected

The carrier freuency is 1Khz .The normalized carrier frequency at sampling rate of 5Khz is 0.4. Thus the spectrum of modulated rectangular pulse is centered at 0.4

Thus significant frequency band lies from 0.3 to 0.5.

If we bandpass filter frequencies in this band,we should get a relatively noiseless waveform,which is expected to improve the performance of correlation.

We apply band-pass butter worth filter of order 2 with normalized cutoff frequencies are 0.2 and 0.6.

def butter_bandpass(lowcut, highcut, fs, order=5):
""" function returns the bandpass butterworth filter coefficients

Parameters
-------------
lowcut,highcut : integer
lower and higher cutoff freuencies in Hz

Fs : Integer
Samping freuency in Hz

order : Integer
Order of butterworth filter

Returns
--------
b,a - numpy-array
filter coefficients

"""
nyq = 0.5 * fs
low = lowcut / nyq
high = highcut / nyq
b, a = butter(order, [low,high], btype='bandpass')
return b, a

def bandpass_filter(data, lowcut, highcut, fs, order,filter_type='butter_worth'):

""" the function performs bandpass filtering

Parameters
-------------
data : numpy-array
input signal

lowcut,highcut : integer
lower and higher cutoff freuencies in Hz

Fs : Integer
Samping freuency in Hz

order : Integer
Order of butterworth filter

Returns
--------
out : numpy-array
Filtered signal

"""
global once
if filter_type=='butter_worth':
b, a = butter_bandpass(lowcut, highcut, fs, order=order)
if once==0:
plt.figure(2)
plotFrequency(b,a)
once=1
y = filtfilt(b, a, data)
return y



The butter-worth band-pass filter is not a zero phase .It will introduce different phase delay depending of frequency. The.The phase plot of the band pass filter is not linear,

To see the phase delay effects of bandpass butter-worth filter,we consider the case of no noise . we can see that instead of 10,delay is introduced as 12 due to butter-worth filter. However we have a stable singular correlation peak being detected

 *********** Information ************
time delay :  10
Noise : 0.001
SNR 46.9897000434
Correct  0  Incorrect  1000
mean  12.0 Std  0.0


.....
#introduce delay
dx=delay(x,tdelay);
result=numpy.zeros((1,loop));
for i in range(loop):
s=dx;

if varnoise!=0:

#normalize signals
s=s/np.linalg.norm(s);
x1=x2/np.linalg.norm(x2);
if carrier!=None:
unfiltered_signal=s;
if mode==1:
s = bandpass_filter(s, carrier-500, carrier+500, Fs,order,filter_type)
if mode==2:
s = bandpass_filter(s, carrier-750, carrier+750, Fs,order,filter_type)

filtered_signal=s;
else:
unfiltered_signal=s;

#perform envelope detection
s=abs(scipy.signal.hilbert(s))
if mode==1 or mode==2:
r=numpy.correlate(s,x1,mode="full")

.......


To compensate for that we need to use zero phase filtering techniques like forward backward filtering.

To achieve zero phase the linear filter is applied twice, once forward and once backwards. The combined filter has zero phase.In general forward-backward filtering squares the amplitude response and zeros the phase response if phase of filter is linear.

we see that there are few errors inspire of applying the forward backward algorithm due to distortions introduced by bandpass filtering of high frequency components.Though the mean is closer to actual delay of 10.

 *********** Information ************
time delay :  10
Noise : 0
Correct  0  Incorrect  1000
mean  9.0 Std  0.0


This tells us that we cannot eliminate the effects introduced due to filtering.

we now introduce noise and compare the performance against envelope detection.

 *********** Information ************
time delay :  10
Noise : 0.3
SNR -2.55272505103
Correct  442  Incorrect  558
mean  9.49 Std  0.678159273327

 *********** Envelope Detection ************
time delay :  10
Noise : 0.3
SNR -2.55272505103
Correct  661  Incorrect  339
mean  9.746 Std  0.793400277288


Compared to just envelope detection,we can see that standard deviation has reduced.Thus the neighborhood of estimated TDOA values are reduced,though the mean is shifted away from 10 may be due to the phase delay effects of bandpass filter.

Though band-pass filtering did not lead to a significant improvement due no the noise within the frequency bandwidth .It is essential to keep out unwanted frequency components due to environmental noise and other factors.

### Code

The code used in the article can be found in github repository “Github Link”

Files

• TDOA1.py
• TODA2.py - Band Pass Filtering

All the plots included in the article can generated by changing the mode variable in the files