# DOA estimation 2 sources time delay estimation uniform array

I am new to signal processing, so sorry if the question is superbanal. I have been spending more days than I would admit on this and can't find the error.

I have two triangle waves that are emitted from loudspeakers and M sensor in a linear uniform array. I need to estimate the delays of the single wave in the single sensor.

I went through the math and understand I have to use the correlation, but when I implement it in Matlab I clearly do some mistakes and get a correlation that has only one peak, in the center (in time) of the signal.

I don't think I can add pictures because I have no reputation on this SO, I hope the question is clear.

EDIT: Apparently I can, so this is the delay per sensor .

In Matlab I basicly do xcorr(X(:,sensor),temp), where X(:,sensor) is the signal received in a sensor, that is the signal received from the two sources plus a small gaussian error and temp is the triangle signal emitted by both sources.

Why does the crosscorrelation have only one peak and not two?

EDIT: full code

clear all
clc
% close all
%dati
v=350;%metri/sec
xa=1; %metri
ya=8;
xb=9;
yb=8;
spazio=0.25;%spazio intra sensori
spt=12.5*10^-3;%secondi
M=41;

%ipotesi
sig=1/10;

%distanze trasmettitore/microfoni
indic_min_dista=xa/spazio+1;%5
indic_min_distb=xb/spazio+1;%37

da=zeros(M,1); %vettore delle distanze
for i=(indic_min_dista):M
da(i)=sqrt(ya^2+(spazio*(i-indic_min_dista))^2);
end

for i=1:(indic_min_dista -1)
da(i)=sqrt(ya^2+(spazio*(indic_min_dista-i))^2);
end

db=zeros(M,1); %vettore delle distanze
for i=(indic_min_distb):M
db(i)=sqrt(yb^2+(spazio*(i-indic_min_distb))^2);
end

for i=1:(indic_min_distb -1)
db(i)=sqrt(yb^2+(spazio*(indic_min_distb-i))^2);
end
%tempo trasmettitore/microfoni

ta=zeros(M,1); %vettore dei tempi
for i=1:M
ta(i)=da(i)/v;
end

tb=zeros(M,1); %vettore dei tempi
for i=1:M
tb(i)=db(i)/v;
end
% ta-(fliplr(tb'))'==0 i tempi sono = se letti dagli estremi opposti
%vettore dei tempi di campionamento

win=10^-3;%secondi
N=ceil((max(ta)+2*spt)/win)+1;
% N=ceil((max(tb)+2*spt)/win)+1; sono = la situazione è del ttt simmetrica

samples=zeros(N,1);
for i=0:(N-1)
samples(i+1)=-spt + win*(i);
end
rit=+spt+samples;

%segnale inviato
S=zeros(N,M);
temp=zeros(N,1);
for i =0:(N-1)
for j=1:M

if((samples(i+1)>ta(j)-spt) && (samples(i+1)<=ta(j)))
S(i+1,j)=S(i+1,j)+1-abs(samples(i+1)-ta(j))/spt;
if(j==1)
temp(i+1)=+1-abs(samples(i+1)-ta(j))/spt;
end
end

if((samples(i+1)>tb(j)-spt) && (samples(i+1)<=tb(j)))
S(i+1,j)=S(i+1,j)+1-abs(samples(i+1)-tb(j))/spt;
end

if((samples(i+1)>ta(j)) && (samples(i+1)<=ta(j)+spt))
S(i+1,j)=S(i+1,j)+1-abs(samples(i+1)-ta(j))/spt;
if(j==1)
temp(i+1)=+1-abs(samples(i+1)-ta(j))/spt;
end
end

if((samples(i+1)>tb(j)) && (samples(i+1)<=tb(j)+spt))
S(i+1,j)=S(i+1,j)+1-abs(samples(i+1)-tb(j))/spt;
end
end
end

figure
imagesc(0:41,rit,S)
xlabel('Sensori');%\sigma^2_\omega
ylabel('Ritardo [s]');

W=zeros(N,M);

%costruisco W
for i =1:N
for j=1:M
W(i,j)=normrnd(0,sig);
end
end

X=W+S;

%vettore che uso per calcolo correlazioni
vv=find(temp~=0);
temp=temp(min(vv): max(vv));

%%
figure
imagesc(0:41,rit,X)
xlabel('Sensori');%\sigma^2_\omega
ylabel('Ritardo [s]');
figure
surf(S)
set(gca,'xlim',[0 41])
xlabel('Sensori'), ylabel('Ritardo [s]');

% ticklabels = get(gca,'YTickLabel');
newTicklabels = rit;%[ticklabels(2:end,:);'60'];
set(gca,'YTickLabel',newTicklabels);

figure
surf(X)
xlabel('Sensori'), ylabel('Ritardo [s]');

% ticklabels = get(gca,'YTickLabel');
newTicklabels = rit;%[ticklabels(2:end,:);'60'];
set(gca,'YTickLabel',newTicklabels);
set(gca,'xlim',[0 41])
set(gca,'zlim',[0 2])

%%
clc
% close all
% Stima ritardo per ogni sensore, in funzione di sigma^2
% valuta MSE cfr CRB almeno m=10 e y=5m,

sig=1/100;%[1/10 1/5 1/2];
Nrun=100;
%segnale vero S(:,41)
%segnale registrato X(:,41)
MSEA=zeros(M,length(sig));
MSEB=zeros(M,length(sig));
mediaA=zeros(M,length(sig));
mediaB=zeros(M,length(sig));

W=zeros(N,M);
psi=(-N:N)*win;%mi serve per valutare phi nella differenza tra i tau

for nrun=1:Nrun
ia =zeros(length(sig),M); %indici del massimo della funzione di correlazione
ib =zeros(length(sig),M); %indici del massimo della funzione di correlazione
for k=1:length(sig)
%costruisco W
for i =1:N
for j=1:M
W(i,j)=normrnd(0,sig(k));
end
end
X=W+S;
for sensor=41:41
phisx=zeros(N+1,1);
phiss=zeros(length(psi),1);%N*2+1
corrtemp=xcorr(X(:,sensor),temp);
figure
plot(corrtemp)
end
end
end


EDIT comparison of the received signal in two sensor (one on the side the blue one) and the central one (green line), the red line corresponds to the original triangle wave.

EDIT Correlation between triangle wave and received signal

In response to the question, this is the hardware positioning. Please remember I am not an art major :D

• Prior to taking the cross-correlation and just plotting the sensors either on a single plot or on a series of subplots, do you see the triangle waves? You should be able to validate your results using visual inspection first to make sure the data is as you expect it. – porten May 5 '14 at 1:18
• yes the wave is triangular everything visually looks ok. I checked that, I really do not see what can be wrong. – Irene May 5 '14 at 7:14
• I'm still a little confused about what your problem is. Your simulation almost looks fine. The only problem I see is that towards the end of the array you have flat top triangles. I think you should have two distinct triangle waves as you go away from the center of the array, with the size of the valley in between them increasing. Also, learn how to vectorize your code and use matrix operations. That should clean and speed it up. – porten May 5 '14 at 12:49
• Sorry, I am not sure I am following you. My problem is that the correlation function should have in theory two peaks, corresponding to the delay from the two sources in every sensor, am I wrong? Instead the correlation (see my edit) presents only one. Am I wrong theoretically or is the code wrong? – Irene May 5 '14 at 15:47
• @porten i had misunderstood your comment. I don't think the simulated signal is wrong because the transmitters are on a line parallel to the array and are closer to the edges of the array, symmetrically in respect to the axis of simmetry of the array. I'll put a picture of it. What do you think? – Irene May 5 '14 at 16:12

The correlation you are seeing looks right - the correlation you are seeing is the superposition of two 'peaks' that are too wide to resolve. Take a look at this example:

%pulses
fs = 10e3;
t1 = -0.1:1/fs:0.1;
t2 = t1 + 100/fs;
w = 40e-3;

%generate two triangular pulses
x = tripuls([t1', t2'],w);

% compute cross correlation
y = [xcorr(x(:,1),x(:,1)) xcorr(x(:,2),x(:,1))];

% plot cross correlation
figure; subplot(3,1,1);
plot(y);
subplot(3,1,2);
plot(sum(y,2));
subplot(3,1,3);
plot(xcorr(sum(x,2),x(:,1)));


Use the convolution instead of the correlation, you will get the two peaks that correspond to the ML estimator of delay time.