DOA - direction of arrival

**colonel_vostrikoff** · 15/05/2014, 18h37

Bonjour à tous,

J'ai un projet de sonar passif, je me suis fait ma carte DSP et ma carte d'acquisition.
L'algo DOA permet de déterminer de quel angle arrive un signal, ce qui est ce que je recherche.
Pour choisir l'algo à implémenter (et déterminer les caractéristiques physiques du sonar) j'ai fait un programme en c pour simuler et faire varier les différents paramètres.
Mais voila : j'avais toujours "deux" signaux !
puis j'ai eu quelques doutes : faut-il des signaux complexes pour calculer ? mes calculs / algos de vecteurs propres sont-ils corrects ?

j'ai découvert Scilab lorsque je cherchais à simuler tout cela pour voir d'ou venait le problème : il me fallait bien reconstruire la partie imaginaire du signal.
ci-dessous mon "programme" scilab pour calculer la direction d'arrivée d'un signal.

les caractéristiques sont les suivantes : 4 micros, dans l'eau, à 44Khz d'échantillonnage, un seul signal qui arrive à un angle de + 25°
L'algo en question est le MUSIC (multiple signal classification) c'est une méthode de sous espace signal qui requiert le calcul de la matrice de covariance.

//----------------- version nombres réels : 2 direction vues alors qu'il y en a qu'une , car les signaux d'entrée ne sont pas complexes ... ---------------------------------------

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
 
clear;
k=0:500;
Fsignal=3000;
Fech = 44000;
Tech = 1/Fech;
Angle = +25;
ECART_MIC_M = 0.1;
VitSon=1500;
Lambda = VitSon / Fsignal
retard = sin(3.141592*Angle/180) * ECART_MIC_M/VitSon;
Mic1= real ( exp( - 2 * %i * 3.141592 * Fsignal *( k(:)*Tech-0*retard)   )  );
Mic2= real ( exp( - 2 * %i * 3.141592 * Fsignal *( k(:)*Tech-1*retard)   )  );
Mic3= real ( exp( - 2 * %i * 3.141592 * Fsignal *( k(:)*Tech-2*retard)   )  );
Mic4= real ( exp( - 2 * %i * 3.141592 * Fsignal *( k(:)*Tech-3*retard)   )  );
 
X=[Mic1, Mic2, Mic3, Mic4]';
R=X*conj(X)';
 
[vect, val]=spec (R);
vp1=[vect(:,1)];
vp2=[vect(:,2)];
vp3=[vect(:,3)];
vp4=[vect(:,4)];
 
//verification
//R*vp1
//1024*vp1
 
 
ValX = 0;
ValY = 0;
Index = 1;
for angle_degre=-90:90
 SinAngle = sin ( angle_degre *3.141592/180 );
 arg = - 2.0 * 3.141592654 * ECART_MIC_M * SinAngle  / Lambda;
 A = [ exp(0); exp(%i*arg); exp(%i*2*arg); exp(%i*3*arg)];
 //A = real ( A );
 denom = conj(A)'*vp1*conj(vp1)'*A   +  conj(A)'*vp2*conj(vp2)'*A ; //  +    conj(A)'*vp4*conj(vp4)'*A;
 disp(denom);
 ValY(Index) = 1 / (0.000000001 + (denom)) ;
 //ValY(Index) = 1 / (0.000000001 + real(denom)) ;) ;
 ValX(Index) = angle_degre;
 Index = Index + 1;
end
plot (ValX, ValY);

//----------------- version nombres complexes : cette fois ci c'est bon ---------------------------------------

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
clear;
k=0:500;
Fsignal=3000;
Fech = 44000;
Tech = 1/Fech;
Angle = +25;
ECART_MIC_M = 0.1;
VitSon=1500;
Lambda = VitSon / Fsignal
retard = sin(3.141592*Angle/180) * ECART_MIC_M/VitSon;
Mic1=exp( - 2 * %i * 3.141592 * Fsignal *( k(:)*Tech-0*retard)   );
Mic2=exp( - 2 * %i * 3.141592 * Fsignal *( k(:)*Tech-1*retard)   );
Mic3=exp( - 2 * %i * 3.141592 * Fsignal *( k(:)*Tech-2*retard)   );
Mic4=exp( - 2 * %i * 3.141592 * Fsignal *( k(:)*Tech-3*retard)   );
 
X=[Mic1, Mic2, Mic3, Mic4]';
R=X*conj(X)';
 
[vect, val]=spec (R);
vp1=[vect(:,1)];
vp2=[vect(:,2)];
vp3=[vect(:,3)];
vp4=[vect(:,4)];
 
//verification
//R*vp1
//1024*vp1
 
 
ValX = 0;
ValY = 0;
Index = 1;
for angle_degre=-90:90
 SinAngle = sin ( angle_degre *3.141592/180 );
 arg = - 2.0 * 3.141592654 * ECART_MIC_M * SinAngle  / Lambda;
 A = [ exp(0); exp(%i*arg); exp(%i*2*arg); exp(%i*3*arg)];
 denom = conj(A)'*vp2*conj(vp2)'*A   +  conj(A)'*vp3*conj(vp3)'*A   +    conj(A)'*vp4*conj(vp4)'*A;
 
 ValY(Index) = 1 / abs (denom);
 ValX(Index) = angle_degre;
 Index = Index + 1;
end
plot (ValX, ValY);

//-----------------------------

voilà
a noter :

- il ne faut pas trop de bruit donc filtrage !!
- sur les signaux réels il faut chercher les directions pour chaque fréquence (d'ou le filtrage et la TF que je n'ai pas faite ici )
- c'est donc bien une méthode de bande étroite : narrow band signal
- la direction d'arrivée avec une simple correlation est une catastrophe car résolution vraiment affreuse et seulement un signal décelable... à éviter !

j'espère que ça pourra servir à quelqu'un ou attiser la curiosité !

**colonel_vostrikoff** · 26/07/2014, 13h40

Version double sources complexe, avec normalisation du résultat

c'est un peu brouillon mais ça me permet de valider l'algo !

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
 
clear;
 
NbEch=100;
MargeEch=50
k=0:NbEch+MargeEch;
Fsignal=1000;
Fsignal2=1500;
Fech = 96000;
Tech = 1/Fech;
Angle = -40;
Angle2 = -35;
ECART_MIC_M = 0.1;
VitSon=1500;
Lambda = VitSon / Fsignal;
Lambda2 = VitSon / Fsignal2;
retard = sin(3.141592*Angle/180) * ECART_MIC_M/VitSon;
retard2 = sin(3.141592*Angle2/180) * ECART_MIC_M/VitSon;
 
Mic1Cplx = exp( - 2 * %i * 3.141592 * Fsignal *( k(:)*Tech-0*retard)    );
 
Mic1 =   exp( - 2 * %i * 3.141592 * Fsignal *( k(:)*Tech-0*retard)   )   +  exp( - 2 * %i * 3.141592 * Fsignal2 *( k(:)*Tech-0*retard2)   ) ;
Mic2 =   exp( - 2 * %i * 3.141592 * Fsignal *( k(:)*Tech-1*retard)   )   +  exp( - 2 * %i * 3.141592 * Fsignal2 *( k(:)*Tech-1*retard2)   ) ;
Mic3 =   exp( - 2 * %i * 3.141592 * Fsignal *( k(:)*Tech-2*retard)   )   +  exp( - 2 * %i * 3.141592 * Fsignal2 *( k(:)*Tech-2*retard2)   ) ;
Mic4 =   exp( - 2 * %i * 3.141592 * Fsignal *( k(:)*Tech-3*retard)   )   +  exp( - 2 * %i * 3.141592 * Fsignal2 *( k(:)*Tech-3*retard2)   ) ;
 
Mic1 =   real ( Mic1 );
Mic2 =   real ( Mic2 );
Mic3 =   real ( Mic3 );
Mic4 =   real ( Mic4 );
 
QuartPeriodeK = (Fech/Fsignal) / 4;
for i=1 : (NbEch + MargeEch - QuartPeriodeK)
  Mic1(i) = Mic1(i) + %i * Mic1(i+QuartPeriodeK);
  Mic2(i) = Mic2(i) + %i * Mic2(i+QuartPeriodeK);
  Mic3(i) = Mic3(i) + %i * Mic3(i+QuartPeriodeK);
  Mic4(i) = Mic4(i) + %i * Mic4(i+QuartPeriodeK);
end
 
QuartPeriodeK = (Fech/Fsignal2) / 4;
for i=1 : (NbEch + MargeEch - QuartPeriodeK)
  Mic1(i) = Mic1(i) + %i * Mic1(i+QuartPeriodeK);
  Mic2(i) = Mic2(i) + %i * Mic2(i+QuartPeriodeK);
  Mic3(i) = Mic3(i) + %i * Mic3(i+QuartPeriodeK);
  Mic4(i) = Mic4(i) + %i * Mic4(i+QuartPeriodeK);
end
 
Mic1 = Mic1(1:NbEch);
Mic2 = Mic2(1:NbEch);
Mic3 = Mic3(1:NbEch);
Mic4 = Mic4(1:NbEch);
 
X=[Mic1, Mic2, Mic3, Mic4]';
R=X*conj(X)';
 
[vect, val]=spec (R);
 
vp1=[vect(:,1)];
vp2=[vect(:,2)];
vp3=[vect(:,3)];
vp4=[vect(:,4)];
 
ValX = 0;
ValY1 = 0;
Index = 1;
maxDOA1 = -100000;
for angle_degre=-90:90
  SinAngle = sin ( angle_degre *3.141592/180 );
  arg = - 2.0 * 3.141592654 * ECART_MIC_M * SinAngle  / Lambda;
  A = [ exp(0); exp(%i*arg); exp(%i*2*arg); exp(%i*3*arg)];
  //PowMUSIC = conj(A)'*vp2*conj(vp2)'*A   +conj(A)'*vp3*conj(vp3)'*A     +conj(A)'*vp4*conj(vp4)'*A  ;
  PowMUSIC = conj(A)'*vp3*conj(vp3)'*A     +conj(A)'*vp4*conj(vp4)'*A  ;
  ValY1(Index) = 1 / abs (PowMUSIC);
  if ( ValY1(Index) > maxDOA1 )
	maxDOA1 = ValY1(Index);
  end
  ValX(Index) = angle_degre;
  Index = Index + 1;
end
Index = 1;
for angle_degre=-90:90
  ValY1(Index) = 65536 * ValY1(Index) / maxDOA1;
  Index = Index + 1;
end
//plot (ValX, ValY1);
 
 
ValX = 0;
ValY2 = 0;
Index = 1;
maxDOA2 = -100000;
for angle_degre=-90:90
  SinAngle = sin ( angle_degre *3.141592/180 );
  arg = - 2.0 * 3.141592654 * ECART_MIC_M * SinAngle  / Lambda2;
  A = [ exp(0); exp(%i*arg); exp(%i*2*arg); exp(%i*3*arg)];
  //PowMUSIC = conj(A)'*vp2*conj(vp2)'*A   +conj(A)'*vp3*conj(vp3)'*A     +conj(A)'*vp4*conj(vp4)'*A  ;
  PowMUSIC = conj(A)'*vp3*conj(vp3)'*A     +conj(A)'*vp4*conj(vp4)'*A  ;
  ValY2(Index) = 1 / abs (PowMUSIC);
  if ( ValY2(Index) > maxDOA2 )
	maxDOA2 = ValY2(Index);
  end
  ValX(Index) = angle_degre;
  Index = Index + 1;
end
Index = 1;
for angle_degre=-90:90
  ValY2(Index) = 65536 * ValY2(Index) / maxDOA2;
  Index = Index + 1;
end
//plot (ValX, ValY2);
 
 
Index = 1;
ValY_tot = 0;
for angle_degre=-90:90
  ValY_tot(Index) = ValY1(Index) + ValY2(Index);
  Index = Index + 1;
end
plot (ValX, ValY_tot);
 
abs(val)

DOA - direction of arrival

Scilab

Discussions similaires

Partager

Partager