Aquí se presenta
el código para el método ESPRIT en Matlab.
ESPRIT es un
método utilizado para estimar la dirección de llegada de una o más señales, mediante
el uso de un array de sensores compuesto por duplas (doublets). Las señales que
son captadas por los sensores cuentan con ruido. Para poder separar este ruido
de la señal, se busca maximizar la proyección de estas señales de entrada en
una base ortogonal, la cual expanda un espacio de igual dimensión al número de
fuentes. Este método funciona
para señales con una banda frecuencial (frequency bandwidth) pequeña, limitación prinicipal respecto a otros métodos como es MUSIC. Para información mas detallada ver la referencia al final de la entrada.
Se presentan dos
códigos. El primero, en el que se ejecuta una realización para distintas
estructuras de arrays, muestra gráficamente la estimación de dos fuentes en el
caso 2D:
clearvars -except RXX_2
clc;
close all;
rng shuffle
twpi = 2.0*pi;
derad = pi / 180.0;
radeg = 180.0 / pi;
c=3*10^8;
% the speed of light
%% 1. Scenarion parameters
% 1.1 Source Definition
w1=2e6;
% Frequency source 1
w2=1e6;
% Frequency source 2
w0=1e9;
% Carrier Frequency
amp=1; % Signal amplitude
angles=[20,-60];
% Sources Angles (degrees)
theta = derad*angles;
lambda0=c*twpi/w0;
SNR=5;
var_n=amp/(10^(SNR/10));
sigma_n=sqrt(var_n);
% 1.2 Array Geometry
N = 9; %
(-) elements number array
delta=[1 0]*lambda0/4;
% distance between doublets
d=lambda0/2;
% Separation between doublets
% delta=[2 0]; %
distance between doublets
% d=4;%lambda0/2;
delta_mod=((delta(1)^2+delta(2)^2)^0.5);
type_array='random';
if strcmp(type_array,'ULA')
% 1.1.1 Linear array
Pos=zeros(2,N);
Pos(1,:)=0:d:(N-1)*d;
Pos(2,:)=zeros(1,N);
elseif strcmp(type_array,'Y')
% 1.1.2 Y-array
angleY=[30 150 270]*derad;
for j=1:3
Pos(1,((j-1)*(N/3)+1):(j*(N/3)))=(d:d:(N/3)*d)*cos(angleY(j));
Pos(2,((j-1)*(N/3)+1):(j*(N/3)))=(d:d:(N/3)*d)*sin(angleY(j));
end;
% Y-array (120º)
% \ /
% \ /
% \ /
% \ /
% +
% |
% |
% |
% |
elseif strcmp(type_array,'square')
% 1.1.3 Square-array
Row=3;
Col=2;
N=Row*Col;
Pos=zeros(2,N);
for j=1:Row
for k=1:Col
Pos(1,(j-1)*Col+k)=(k-1)*d;
Pos(2,(j-1)*Col+k)=(j-1)*d;
end;
end;
% Square-array (5x4)
% +---+---+---+
% | |
| |
% +---+---+---+
% | |
| |
% +---+---+---+
% | |
| |
% +---+---+---+
% | |
| |
% +---+---+---+
elseif strcmp(type_array,'random')
% 1.1.4 random-array
N=5;
for j=1:N
Pos(1,j)=d*rand;
Pos(2,j)=d*rand;
end;
elseif strcmp(type_array,'paper')
% 1.1.5 Linear array and scenario Paper (pag. 993)
N = 5;
Pos=zeros(2,N);
Pos(1,:)=[0,1,3,5.5,7]*(lambda0/2);
delta=[1 0]*(lambda0/4);
delta_mod=((delta(1)^2+delta(2)^2)^0.5);
angles=[24 28];
theta = derad*angles;
SNR=0;
var_n=amp/(10^(SNR/10));
sigma_n=sqrt(var_n);
% Take into account Rayleigh of 3dBs
end
% 1.3 Simulation parameters
n=100;
% sample number
N_cov=1;
% Number of meassures for calculate
the cov matrix
% 1.4. Generate Input Sampled Signal
Ts=1e6;
% 1MHz
Fs=1/Ts;
t=linspace(0,Ts,n);
phase=0.5;
% The sources don't arrive at the
same moment
s1=amp*exp(1i*w1*t);
s2=amp*exp(1i*w2*t);
s=[s1.*exp(w0*1i*t);s2.*exp(w0*1i*t+phase)];
theta_delta=atan(delta(2)/delta(1));
% We use the relative angle between the direction of arrival and the
delta vector
gamma=(w0*delta_mod/c)*[sin(theta(1)-theta_delta)
sin(theta(2)-theta_delta)];
gamma1=gamma;
tau=zeros(2,N);
% Genereal case
for i=1:N
tau(1,i)=(Pos(:,1)-Pos(:,i))'*[-sin(theta(1)); cos(theta(1))]; % /c
tau(2,i)=(Pos(:,1)-Pos(:,i))'*[-sin(theta(2));
cos(theta(2))]; % /c
end;
i=sqrt(-1);
a1 = exp(2i*pi*tau(1,:)');
a2 = exp(-2i*pi*tau(2,:)');
A=[a1 a2];
phi=diag(exp(1i*gamma));
%% 2. Calculate the cov matrix of z=[x;y]
Z=zeros(2*N,n,N_cov);
RZZ=zeros(2*N,2*N);
x = zeros(N,n,N_cov);
y=x;
for j=1:N_cov
x(:,:,j)=A*s+sigma_n*randn(N,n);
y(:,:,j)=A*phi*s+sigma_n*randn(N,n);
Z=[x(:,:,j);y(:,:,j)];
RZZ=RZZ+Z*Z';
end
RZZ=RZZ/(N*N_cov);
[U,D,V] = svd(RZZ);
%% 3. Estimate number of sources
% We keep at least the 95% of the energy
D2=diag(D*conj(D));
Energy=sum(D2);
threshold=0.98;
D2_cumu=0;
j=0;
while (D2_cumu<(threshold*Energy))
j=j+1;
D2_cumu=D2_cumu+D2(j);
end;
num_sources=j;
% num_sources=2;
%% 3. Determine Sources
ES=U(:,1:num_sources);
EX=ES(1:N,:);
EY=ES(N+1:2*N,:);
psi=pinv(EX)*EY; %??
phi=eigs(psi);
gamma_est=angle(phi);
angles_est=asin(gamma_est*c/(w0*delta_mod))*radeg;
%% 4. Result Display
% -------------------------Graphic configuration---------------------------
% Propiedades de las figuras (cm)
% Size
width = 18; % 18 para DIN A4
height = 17;
width_line=2;
width_dot=3;
% Margins
margin_left = 1.5; %0.95;
margin_right = 0.25; %0.25;
margin_top = 3;
margin_bottom = 1.2;
% Font
font_name = 'Times'; % Times,
Helvetica
font_size = 12.0; % Tamaño de fuente (ejes)
legend_font_size = 10.0; % Tamaño de
fuente (leyenda)
title_font_size = 16.0; % Tamaño de
fuente del título (leyenda)
marker_size = 6; % Tamaño de los
marcadores
color=[[1 0 0];[0 1 0];[0 0 1];[0 1 1];[1 0 1];[1 1 0]];
% -------------------end graphic
configuration-----------------------------
r=0:1:10;
% 4.1 Array geometry
figure(1)
clf;
position = [1 2 width height];
set(gcf,'Units','centimeters',...
'Position',position,...
'PaperUnits','centimeters',...
'PaperPosition',position,...
'PaperPositionMode','auto');
% Crear ejes
axesposition = [margin_left, margin_bottom,
position(3)-margin_left-margin_right position(4)-margin_bottom-margin_top];
axes('FontSize',font_size,...
'Units','centimeters',...
'FontName',font_name,...
'Position',axesposition);
hold on; grid on;
plot(Pos(1,:),Pos(2,:),'+','LineWidth',width_dot,'Color','black');
hold on;
plot(Pos(1,:)+delta(1),Pos(2,:)+delta(2),'x','LineWidth',width_dot,'Color','black');
% 4.2 Sources
mid_point=[(min(Pos(1,:))+max(Pos(1,:)))/2
(min(Pos(2,:))+max(Pos(2,:)))/2];
for j=1:num_sources
plot(mid_point(1)-r*sin(angles(j)*derad),mid_point(2)+r*cos(angles(j)*derad),'-.r','LineWidth',width_line);
end;
% 4.3 Sources Estimation
for j=1:num_sources
plot(mid_point(1)-r*sin(angles_est(j)*derad),mid_point(2)+r*cos(angles_est(j)*derad),'LineWidth',width_line);
end;
axis([min(Pos(1,:))-2
max(Pos(1,:))+2 min(Pos(2,:))-2 max(Pos(2,:))+2]);
title(sprintf('Direction of
Arrival Estimation (SNR %3.1f)\n(Source 1: %3.1f Source Est. 1: %3.1f)
\n(Source 2: %3.1f Source Est. 1: %3.1f)\n',...
SNR,angles(1), angles_est(1), angles(2),
angles_est(2)),...
'interpreter','latex',...
'HorizontalAlignment','center',...
'FontSize',title_font_size);
my_legend = {'x doublet','y doublet','Source 1','Source 2','Est. Source
1','Est. Source 2'};
aux=legend(my_legend,...
'Interpreter','latex',...
'FontSize',12,...
'Orientation','Vertical',...
'Units','centimeters',...
'Location','southeast');
El segundo,
muestra un análisis de Monte-Carlo para distintas configuraciones de ULAs en un
escenario de SNR=5. Se hace variar el tamaño del array viendo cómo varían los resultados.
%% Montecarlo Analysis for ESPRIT estimation
clearvars;
clc;
close all;
rng shuffle; % Initializate the seed of random numbers
num_simu=3000;
% Simulations number
twpi = 2.0*pi;
derad = pi / 180.0;
radeg = 180.0 / pi;
c=3*10^8;
% the speed of light
%% 1. Scenarion parameters
% 1.1 Source Definition
w1=2e6;
% Frequency source 1
w2=1e6;
% Frequency source 2
w0=1e9;
% Carrier Frequency
amp=1;
% Signal amplitude
angles=[24,28];
% Sources Angles (degrees)
theta = derad*angles;
lambda0=c*twpi/w0;
SNR=10;
% We set the SNR
var_n=(amp^2)/(10^(SNR/10));
% We calculate the variance
sigma_n=sqrt(var_n);
% 1.2 Array Geometry
delta=[1 0]*lambda0/4;
% distance between doublets
d=lambda0/2;
% Separation between doublets
delta_mod=((delta(1)^2+delta(2)^2)^0.5);
%% We compare the three architectures
num_architectures=4;
angles_est=zeros(2,num_simu,num_architectures);
names={'ULA (2 Doublets)','ULA2 (4
Doublets)','ULA3
(6 Doublets)','ULA4
(8 Doublets)'};
N=0;
for comp=1:num_architectures
type_array=names(comp)
N=N+2;
% 1.1.1 Linear array 1
d=5;
Pos=zeros(2,N);
Pos(1,:)=0:d:(N-1)*d;
Pos(2,:)=zeros(1,N);
% 1.3 Simulation parameters
n=100; % sample number
N_cov=1; % Number of meassures for calculate the cov matrix
% 1.4. Generate Input Sampled Signal
Ts=1e6; % 1MHz
Fs=1/Ts;
t=linspace(0,Ts,n);
phase=0.5; % The sources don't arrive at the same moment
s1=amp*exp(1i*w1*t);
s2=amp*exp(1i*w2*t);
s=[s1.*exp(w0*1i*t);s2.*exp(w0*1i*t+phase*1i)];
theta_delta=atan(delta(2)/delta(1));
% We use the relative angle between the direction of arrival and the
delta vector
gamma=(w0*delta_mod/c)*[sin(theta(1)-theta_delta)
sin(theta(2)-theta_delta)];
gamma1=gamma;
tau=zeros(2,N); % Genereal case
for i=1:N
tau(1,i)=(Pos(:,1)-Pos(:,i))'*[-sin(theta(1));
cos(theta(1))]; % /c
tau(2,i)=(Pos(:,1)-Pos(:,i))'*[-sin(theta(2)); cos(theta(2))]; % /c
end;
i=sqrt(-1);
a1 = exp(2i*pi*tau(1,:)');
a2 = exp(-2i*pi*tau(2,:)');
A=[a1 a2];
phi=diag(exp(1i*gamma));
%% Simulation Star
for simu=1:num_simu
%% 2. Calculate the cov
matrix of z=[x;y]
Z=zeros(2*N,n,N_cov);
RZZ=zeros(2*N,2*N);
x = zeros(N,n,N_cov);
y=x;
for j=1:N_cov
x(:,:,j)=A*s+sigma_n*randn(N,n);
y(:,:,j)=A*phi*s+sigma_n*randn(N,n);
Z=[x(:,:,j);y(:,:,j)];
RZZ=RZZ+Z*Z';
end
RZZ=RZZ/(N*N_cov);
[U,D,V] = svd(RZZ);
%% 3. Estimate number of
sources
% We dont estimate here the number sources
num_sources=2;
%% 4. Determine Sources
ES=U(:,1:num_sources);
EX=ES(1:N,:);
EY=ES(N+1:2*N,:);
psi=pinv(EX)*EY;
phi_est=eigs(psi);
gamma_est=angle(phi_est);
angles_est(:,simu,comp)=asin(gamma_est*c/(w0*delta_mod))*radeg;
end;
end;
%% Plot histogram
% -------------------------Graphic
configuration---------------------------
% Propiedades de las figuras (cm)
% Size
width = 18; % 18 para DIN A4
height = 17;
% Margins
margin_left = 1.5; %0.95;
margin_right = 0.25; %0.25;
margin_top = 2.6;
margin_bottom = 1.2;
% Font
font_name = 'Times'; % Times,
Helvetica
font_size = 12.0; % Tamaño de fuente (ejes)
legend_font_size = 10.0; % Tamaño de
fuente (leyenda)
title_font_size = 16.0; % Tamaño de
fuente del título (leyenda)
marker_size = 6; % Tamaño de los
marcadores
color=[[1 0 0];[0 1 0];[0 0 1];[0 1 1];[1 0 1];[1 1 0]];
% -------------------end graphic
configuration-----------------------------
num_div=100;
figure(1)
clf;
position = [1 2 width height];
set(gcf,'Units','centimeters',...
'Position',position,...
'PaperUnits','centimeters',...
'PaperPosition',position,...
'PaperPositionMode','auto');
% Crear ejes
axesposition = [margin_left, margin_bottom,
position(3)-margin_left-margin_right position(4)-margin_bottom-margin_top];
axes('FontSize',font_size,...
'Units','centimeters',...
'FontName',font_name,...
'Position',axesposition);
hold on; grid on;
for j=1:num_architectures
[counts,centers] =
hist(reshape(squeeze(angles_est(:,:,j)),[1,2*num_simu]),num_div);
plot(centers,counts,'Color',color(j,:),'LineWidth',2);
hold on;
end;
aux=title(sprintf('Histogram
Estimation DOA (%i Simulations)\n(Source 1: %3.1f Source 2: %3.1f)\n(SNR %3.1f and %i
snapshots)',...
num_simu,angles(1),
angles(2),SNR,n),...
'interpreter','latex',...
'HorizontalAlignment','center',...
'FontSize',title_font_size);
my_legend = names;
aux=legend(my_legend,...
'Interpreter','latex',...
'FontSize',legend_font_size,...
'Orientation','Vertical',...
'Units','centimeters',...
'Location','northwest');
aux = xlabel('$^\circ$ DOA',...,
'Interpreter','latex',...
'FontSize',font_size,...
'Units','centimeters',...
'VerticalAlignment','bottom');
set(aux,'Position',get(aux,'Position') - [0 .3 0])
aux = ylabel('N$^\circ$ of
simulations',...,
'Interpreter','latex',...
'FontSize',font_size,...
'Units','centimeters',...
'VerticalAlignment','bottom');
(Nota: Corregir saltos de línea del código)
Referencias: ESPRIT-Estimation of Signal Parameters Via Rotational
Invariance Techniques RICHARD ROY AND THOMAS KAILATH, FELLOWI, IEEE
No hay comentarios:
Publicar un comentario