Técnica ESPRIT para la estimación de DOA en Matlab

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