v_paper_4_2.m

极化阵列信号处理DOA及极化参数估计； 阵列同何结构不作要求

% three vector sensors 
% the forth chapter
clear all;
close all;
clc;
warning off;

j_2p=j*2*pi;
degrad=pi/180;
snapshot=800;
m=3;
num_sig=2;
N=3;


a(:,1)=v_a(30.93,37.09,45,90);
a(:,2)=v_a(50.08,39.71,45,-90);

s(1,:)=exp(j_2p*5/25*[0:snapshot-1]);
s(2,:)=exp(j_2p*7/25*[0:snapshot-1]);

%Poynting vectors
p(:,1)=cross(a(1:3,1),conj(a(4:6,1)));
p(:,2)=cross(a(1:3,2),conj(a(4:6,2)));

% 
for d=1:1:4
%d=0.5
% position and phase delay
p_s=[0 0 0;d 0 0;0 d 0];
for i=1:m
    q(i,1)=exp(j_2p*p_s(i,1)*p(1,1))*exp(j_2p*p_s(i,2)*p(2,1))*exp(j_2p*p_s(i,3)*p(3,1));
    q(i,2)=exp(j_2p*p_s(i,1)*p(1,2))*exp(j_2p*p_s(i,2)*p(2,2))*exp(j_2p*p_s(i,3)*p(3,2));
end

% receive signals
%###########
for SNR=0:5:30
%SNR=10
    for monte_loop=1:50
A=[kron(q(:,1),a(:,1)),kron(q(:,2),a(:,2))];
%@@@@@@@@
A_kt=[kron(a(:,1),q(:,1)),kron(a(:,2),q(:,2))];

z_kt=A_kt*s(1:2,:)+(10^(-(SNR)/20))*v_noise(6*m,snapshot)/sqrt(2);

R=1/snapshot*z_kt*z_kt';

%获得信号子空间的6个子空间
[u,v]=eigs(R,size(R,1));
Es=u(:,1:num_sig);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u2=[];
v2=[];
for loop=1:5
    %Es1=Es(N*(loop-1)+1:N*loop,:);
    Es1=Es(1:N,:);
    Es2=Es(loop*N+1:N*(1+loop),:);
    %vary_Es3=inv(Es1'*Es1)*Es1'*Es2;
    %[T1 v1]=eigs(vary_Es3,num_sig);
    %u1=1/2*(Es1*inv(T1)+Es2*inv(T1)*inv(v1));
    [u1,v1]=eigs(Es2*inv(Es1'*Es1)*Es1',size(Es2,1));
    u1=u1(:,1:num_sig);
    v1=v1(1:num_sig,1:num_sig);
%%%%%%%%%%%%%%%%%%%%%solve the permutation,pairing%%%%%%%%     
   if (loop~=1)         
    varya=u2(:,1:num_sig);
    varyb=u1;
    varyc=abs(varya'*varyb);
    [value1 index1]=max(varyc.');
    varyd=diag(v1);
    varye=varyd(index1);
    v1=diag(varye);
    u1=u1(:,index1);
   end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    u2=[u2 u1];  
    v2=[v2 v1];  
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%angle estimation%%%%%%%%%%%%%%%
for loop=1:num_sig
    varyf=u2(:,loop:num_sig:size(u2,2));
    varyg=v2(loop,loop:num_sig:size(v2,2));
    Ex=[1 varyg(1:2)];
    Hx=conj(varyg(3:5));
    direction_cos=[Ex(2)*Hx(3)-Ex(3)*Hx(2) Ex(3)*Hx(1)-Ex(1)*Hx(3) Ex(1)*Hx(2)-Ex(2)*Hx(1)]/(norm(Ex)*norm(Hx));
    angle_esti(loop,:)=[asin(norm(direction_cos(1:2))) atan(real(direction_cos(2)/direction_cos(1)))]/pi*180;
end                      %%elevation                            azimuth
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%polarization estimation
for loop=1:num_sig
    elvation=angle_esti(loop,1)/180*pi;
    azimuth=angle_esti(loop,2)/180*pi;
    electro_esti=[cos(azimuth)*cos(elvation)      -sin(azimuth);...
                  sin(azimuth)*cos(elvation)      cos(azimuth);...
                  -sin(elvation)                     0                 ;...
                  -sin(azimuth)                    -cos(azimuth)*cos(elvation);...
                  cos(azimuth)                     -sin(azimuth)*cos(elvation);...
                  0                                 sin(elvation)           ];
     
    varyg=v2(loop,loop:num_sig:size(v2,2));
    a_esti=[1 varyg].';
    g_esti=inv(electro_esti'*electro_esti)*electro_esti'*a_esti;
    %gg=[gg g]
    pola_esti(loop,:)=[atan(abs(g_esti(1)/g_esti(2))) ANGLE(g_esti(1))-ANGLE(g_esti(2))]/pi*180;
end
ang_kt(:,:,monte_loop)=[angle_esti,pola_esti];
% sort
if sin(ang_kt(1,1)*degrad)>sin(ang_kt(2,1)*degrad)
    temp=ang_kt(1,:);
    ang_kt(1,:)=ang_kt(2,:);
    ang_kt(2,:)=temp;
end
p_kt(:,1)=v_p(ang_kt(1,1),ang_kt(1,2));
p_kt(:,2)=v_p(ang_kt(2,1),ang_kt(2,2));





%% my ideal
%$$$$$$$$$$$$$$$$$$$$
z=A*s+(10^(-(SNR)/20))*v_noise(6*m,snapshot)/sqrt(2);
R=1/snapshot*z*z';
[u v w]=svd(R);

eo=u(1:6,1:2);
ex=u(7:12,1:2);
ey=u(13:18,1:2);
u=[eo;ex;ey];

wxo=inv(eo'*eo)*eo'*ex;
wyo=inv(eo'*eo)*eo'*ey;

[uu uv]=eig(wxo);
[vu vv]=eig(wyo);

p_e(1,:)=[angle(uv(1,1)),angle(uv(2,2))]./(2*pi*d);
p_e(2,:)=[angle(vv(1,1)),angle(vv(2,2))]./(2*pi*d);

%%% search solve spatio 
%$$$$$$$$$$$$$$$$$$$$$
if d>0.5
for ii=1:2
    d_te=floor(d*(-1-p_e(1,ii)));
    e_te=ceil(d*(1-p_e(1,ii)));
    b_u=d_te:e_te;
    u_temp=p_e(1,ii)+b_u/d;
    uu_te=u_temp-p_kt(1,ii);
    [a_te b_te]=min(abs(uu_te));
    p_e(1,ii)=u_temp(b_te);
    
    f_te=floor(d*(-1-p_e(2,ii)));
    g_te=ceil(d*(1-p_e(2,ii)));
    b_v=f_te:g_te;
    v_temp=p_e(2,ii)+b_v/d;
    vv_te=v_temp-p_kt(2,ii);
    [a_te b_te]=min(abs(vv_te));
    p_e(2,ii)=v_temp(b_te);
end
if p_e(1,1)-p_kt(1,1)>0.22
    temp=p_e(:,1);
    p_e(:,1)=p_e(:,2);
    p_e(:,2)=temp;
    for ii=1:2
    d_te=floor(d*(-1-p_e(1,ii)));
    e_te=ceil(d*(1-p_e(1,ii)));
    b_u=d_te:e_te;
    u_temp=p_e(1,ii)+b_u/d;
    uu_te=u_temp-p_kt(1,ii);
    [a_te b_te]=min(abs(uu_te));
    p_e(1,ii)=u_temp(b_te);
    
    f_te=floor(d*(-1-p_e(2,ii)));
    g_te=ceil(d*(1-p_e(2,ii)));
    b_v=f_te:g_te;
    v_temp=p_e(2,ii)+b_v/d;
    vv_te=v_temp-p_kt(2,ii);
    [a_te b_te]=min(abs(vv_te));
    p_e(2,ii)=v_temp(b_te);
end
end
end
%$$$$$$$$$$$$$$$$$$$$$

ang_half(1,1)=asin(norm(p_e(:,1)))*180/pi;
ang_half(1,2)=atan(p_e(2,1)/p_e(1,1))*180/pi;
ang_half(2,1)=asin(norm(p_e(:,2)))*180/pi;
ang_half(2,2)=atan(p_e(2,2)/p_e(1,2))*180/pi;



h1=[cos(ang_half(1,2)*pi/180)*cos(ang_half(1,1)*pi/180)    -sin(ang_half(1,2)*pi/180)
    sin(ang_half(1,2)*pi/180)*cos(ang_half(1,1)*pi/180)     cos(ang_half(1,2)*pi/180)
   -sin(ang_half(1,1)*pi/180)                                   0
   -sin(ang_half(1,2)*pi/180)                              -cos(ang_half(1,2)*pi/180)*cos(ang_half(1,1)*pi/180)
    cos(ang_half(1,2)*pi/180)                              -sin(ang_half(1,2)*pi/180)*cos(ang_half(1,1)*pi/180)
    0                                                       sin(ang_half(1,1)*pi/180)];

h2=[cos(ang_half(2,2)*pi/180)*cos(ang_half(2,1)*pi/180)    -sin(ang_half(2,2)*pi/180)
    sin(ang_half(2,2)*pi/180)*cos(ang_half(2,1)*pi/180)     cos(ang_half(2,2)*pi/180)
   -sin(ang_half(2,1)*pi/180)                                   0
   -sin(ang_half(2,2)*pi/180)                              -cos(ang_half(2,2)*pi/180)*cos(ang_half(2,1)*pi/180)
    cos(ang_half(2,2)*pi/180)                              -sin(ang_half(2,2)*pi/180)*cos(ang_half(2,1)*pi/180)
    0                                                       sin(ang_half(2,1)*pi/180)];



ang_v(:,:,monte_loop)=ang_half;

end
% end of monte loop

% kt ideal   mean

% v ideal
re_v_1_mean((SNR+5)/5,1,(d))=mean(ang_v(1,1,:));
re_v_1_mean((SNR+5)/5,2,(d))=mean(ang_v(1,2,:));

re_v_2_mean((SNR+5)/5,1,(d))=mean(ang_v(2,1,:));
re_v_2_mean((SNR+5)/5,2,(d))=mean(ang_v(2,2,:));

%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%*********************************
% kt bias
re_kt_1_bias((SNR+5)/5,1,(d))=v_std(ang_kt(1,1,:),monte_loop);
re_kt_1_bias((SNR+5)/5,2,(d))=v_std(ang_kt(1,2,:),monte_loop);

re_kt_2_bias((SNR+5)/5,1,(d))=v_std(ang_kt(2,1,:),monte_loop);
re_kt_2_bias((SNR+5)/5,2,(d))=v_std(ang_kt(2,2,:),monte_loop);


% v ideal bias
re_v_1_bias((SNR+5)/5,1,(d))=v_std(ang_v(1,1,:),monte_loop);
re_v_1_bias((SNR+5)/5,2,(d))=v_std(ang_v(1,2,:),monte_loop);


re_v_2_bias((SNR+5)/5,1,(d))=v_std(ang_v(2,1,:),monte_loop);
re_v_2_bias((SNR+5)/5,2,(d))=v_std(ang_v(2,2,:),monte_loop);

end 
%  end SNR
end % end d

% figure
SNR=0:5:30
figure
semilogy(SNR,re_kt_1_bias(:,1,1)','-*');
hold on;
semilogy(SNR,re_kt_1_bias(:,2,1)','-square');
hold on;

semilogy(SNR,re_v_1_bias(:,1,1),'-diamond');
hold on;
semilogy(SNR,re_v_1_bias(:,2,1),'->');
xlabel('信噪比(dB)');
ylabel('标准偏差(deg)');
title('俯仰角估计值的标准偏差与信噪比关系');
%%%%
figure
semilogy(SNR,re_v_1_bias(:,1,1),'-*');
hold on;
semilogy(SNR,re_v_1_bias(:,1,2),'-square');
hold on;

figure
semilogy(SNR,re_v_1_bias(:,2,1),'-*');
hold on;
semilogy(SNR,re_v_1_bias(:,2,2),'-square');
xlabel('信噪比(dB)');
ylabel('标准偏差(deg)');
title('方位角估计值的标准偏差与信噪比关系');
hold on;