v_paper_3_1.m

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

%两个信号，13个矢量传感器


close all;
clear all;
clc;
warning off;

j_2p=j*2*pi;
degrad=pi/180;
snap_num=800;
num_sig=2;
N=13;
%入射信号的参数信息
a(:,1)=v_a(58.07,44.05,45,90);
a(:,2)=v_a(59.1,36.47,45,-90);


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


%构造信号源
r_1=(1+j)*randn(1,snap_num)/sqrt(2);
r_2=(1+j)*randn(1,snap_num)/sqrt(2);

s(1,:)=r_1.*exp(j*5/25*(1:snap_num)+pi/3);
s(2,:)=r_2.*exp(j*7/25*(1:snap_num)+pi/6);

%矢量天线空间相位矢量
p_s=[0 0 0;0.5 0 0;0 0.5 0;1.35 0 0;-1.35 0 0;0 0.5 0;0 -0.5 0;0 1.35 0;0 -1.35 0;2 2 0;2 -2 0;-2 -2 0;-2 2 0];
for i=1:N
    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


%构建阵列流型
A=[kron(a(:,1),q(:,1)),kron(a(:,2),q(:,2))];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%获得快拍后的数据
for SNR=-10:5:40
for monte_loop=1:500
z=A*s+(10^(-(SNR)/20))*v_noise(6*N,snap_num)/sqrt(2);
Rx=1/snap_num*z*z';
[u,v]=eigs(Rx,size(Rx,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];  %%N,5*num_sig
    v2=[v2 v1];  %%num_sig,num_sig*5
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=[angle_esti,pola_esti];
if sin(ang(1,2)*degrad)<sin(ang(2,2)*degrad)
    temp=ang(1,:);
    ang(1,:)=ang(2,:);
    ang(2,:)=temp;
end
re_ang(:,:,monte_loop)=ang;
p_es(1,1,monte_loop)=sin(ang(1,1)*degrad)*cos(ang(1,2)*degrad);
p_es(2,1,monte_loop)=sin(ang(1,1)*degrad)*sin(ang(1,2)*degrad);

p_es(1,2,monte_loop)=sin(ang(2,1)*degrad)*cos(ang(2,2)*degrad);
p_es(2,2,monte_loop)=sin(ang(2,1)*degrad)*sin(ang(2,2)*degrad);
end

bias_p_1(1,(SNR+15)/5)=abs(max(p_es(1,1,:)-p(1,1)));
bias_p_1(2,(SNR+15)/5)=abs(max(p_es(2,1,:)-p(2,1)));

bias_p_2(1,(SNR+15)/5)=abs(max(p_es(1,2,:)-p(1,2)));
bias_p_2(2,(SNR+15)/5)=abs(max(p_es(2,2,:)-p(2,2)));
%%%%std
std_ang_1(1,(SNR+15)/5)=v_std(re_ang(1,1,:),monte_loop);
std_ang_1(2,(SNR+15)/5)=v_std(re_ang(1,2,:),monte_loop);
std_ang_1(3,(SNR+15)/5)=v_std(re_ang(1,3,:),monte_loop);
std_ang_1(4,(SNR+15)/5)=v_std(re_ang(1,4,:),monte_loop);

std_ang_2(1,(SNR+15)/5)=v_std(re_ang(2,1,:),monte_loop);
std_ang_2(2,(SNR+15)/5)=v_std(re_ang(2,2,:),monte_loop);
std_ang_2(3,(SNR+15)/5)=v_std(re_ang(2,3,:),monte_loop);
std_ang_2(4,(SNR+15)/5)=v_std(re_ang(2,4,:),monte_loop);
end

% figure
SNR=-10:5:40
figure
semilogy(SNR,bias_p_1(1,:),'-*');
hold on;
semilogy(SNR,bias_p_1(2,:),'-s');
hold on;
semilogy(SNR,bias_p_2(1,:),'-diamond');
hold on;
semilogy(SNR,bias_p_2(2,:),'->');
xlabel('信噪比(dB)');
ylabel('最大偏差');
title('入射信号U，V估计值的最大偏差与信噪比关系');
hold on;

figure
semilogy(SNR,std_ang_1(1,:),'-*');
hold on;
semilogy(SNR,std_ang_1(2,:),'-s');
hold on;
semilogy(SNR,std_ang_2(1,:),'-diamond');
hold on;
semilogy(SNR,std_ang_2(2,:),'->');
xlabel('信噪比(dB)');
ylabel('标准偏差(deg)');
title('入射信号空间到达角估计值的标准偏差与信噪比关系');
hold on;

figure
semilogy(SNR,std_ang_1(3,:),'-*');
hold on;
semilogy(SNR,std_ang_1(4,:),'-s');
hold on;
semilogy(SNR,std_ang_2(3,:),'-diamond');
hold on;
semilogy(SNR,std_ang_2(4,:),'->');
xlabel('信噪比(dB)');
ylabel('标准偏差(deg)');
title('入射信号极化角估计值的标准偏差与信噪比关系');
hold on;