📄 v_paper_3_2.m

📁 极化阵列信号处理DOA及极化参数估计；阵列同何结构不作要求
💻 M
字号:
% three vector sensors ,two signals,the third character
close all;
clear all;
clc;
warning off;

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


%入射信号的参数信息
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);

%矢量天线空间相位矢量
for d=1:2:7
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

for SNR=-10:5:30
    for monte_loop=1:500
%构建阵列流型
A=[kron(a(:,1),q(:,1)),kron(a(:,2),q(:,2))];
A_new=[kron(q(:,1),a(:,1)),kron(q(:,2),a(:,2))];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the first signal
z1=A(:,1)*s(1,:)+(10^(-(SNR)/20))*v_noise(6*m,snap_num)/sqrt(2);
R=1/snap_num*z1*z1';
[es bb cc]=svd(R);

Es=es(:,1);
Es_1=v_j(1,m)*Es;
Es_2=v_j(2,m)*Es;
Es_3=v_j(3,m)*Es;
Es_4=v_j(4,m)*Es;
Es_5=v_j(5,m)*Es;
Es_6=v_j(6,m)*Es;

%获得5个旋转不变量
w1_2=inv(Es_1'*Es_1)*(Es_1'*Es_2);
w2_3=inv(Es_2'*Es_2)*(Es_2'*Es_3);
w3_4=inv(Es_3'*Es_3)*(Es_3'*Es_4);
w4_5=inv(Es_4'*Es_4)*(Es_4'*Es_5);
w5_6=inv(Es_5'*Es_5)*(Es_5'*Es_6);

%获得阵列导向矢量的五个不变量
[u x1_2]=eig(w1_2);
[u x2_3]=eig(w2_3); 
[u x3_4]=eig(w3_4); 
[u x4_5]=eig(w4_5); 
[u x5_6]=eig(w5_6); 

%获得信号Poynting矢量的估计值
tmp_1_1=[1;x1_2(1,1);x1_2(1,1)*x2_3(1,1)];
tmp_1_2=[x1_2(1,1)*x2_3(1,1)*x3_4(1,1);x1_2(1,1)*x2_3(1,1)*x3_4(1,1)*x4_5(1,1);x1_2(1,1)*x2_3(1,1)*x3_4(1,1)*x4_5(1,1)*x5_6(1,1)];
p_e1=cross(tmp_1_1,conj(tmp_1_2))/sqrt(tmp_1_1'*tmp_1_1)/sqrt(tmp_1_2'*tmp_1_2);
%仰俯角和方位角,1为仰俯角，2为方位角
ang(1,1)=(acos(p_e1(3))*180/pi);
ang(1,2)=atan((p_e1(2)/p_e1(1)))*180/pi;
h1=[cos(ang(1,2)*pi/180)*cos(ang(1,1)*pi/180) -sin(ang(1,2)*pi/180)
    sin(ang(1,2)*pi/180)*cos(ang(1,1)*pi/180)  cos(ang(1,2)*pi/180)
   -sin(ang(1,1)*pi/180)                        0
   -sin(ang(1,2)*pi/180)                       -cos(ang(1,2)*pi/180)*cos(ang(1,1)*pi/180)
    cos(ang(1,2)*pi/180)                       -sin(ang(1,2)*pi/180)*cos(ang(1,1)*pi/180)
    0                                           sin(ang(1,1)*pi/180)];
a1=1/sqrt((tmp_1_1'*tmp_1_1)*(tmp_1_2'*tmp_1_2))*[1;x1_2(1,1);x1_2(1,1)*x2_3(1,1);x1_2(1,1)*x2_3(1,1)*x3_4(1,1);x1_2(1,1)*x2_3(1,1)*x3_4(1,1)*x4_5(1,1);x1_2(1,1)*x2_3(1,1)*x3_4(1,1)*x4_5(1,1)*x5_6(1,1)];
G1=inv(h1'*h1)*h1'*a1;
%计算两信号的参数，，3为辅助极化角，4为相位差
ang(1,3)=(atan(abs(G1(1)/G1(2)))*180/pi);
ang(1,4)=((angle(G1(1))-angle(G1(2)))*180/pi);
ang_1(monte_loop,:)=ang;
% new method
z1_new=A_new(:,1)*s(1,:)+(10^(-(SNR)/20))*v_noise(6*m,snap_num)/sqrt(2);
R_new=1/snap_num*z1_new*z1_new';
[aa bb cc]=svd(R_new);
Eo=aa(1:6,1);
Ex=aa(7:12,1);
Ey=aa(13:18,1);

wu=inv(Eo'*Eo)*Eo'*Ex;
wv=inv(Eo'*Eo)*Eo'*Ey;

[qq u]=eig(wu);
[ww v]=eig(wv);

u=angle(u)/(2*pi*d);
v=angle(v)/(2*pi*d);

if d>0.5
    bu=ceil(d*(-1-u)):floor(d*(1-u));
    bv=ceil(d*(-1-v)):floor(d*(1-v));
    temp_u=u+bu/d;
    temp_v=v+bv/d;
    [qu ii_u]=min(abs(temp_u-p_e1(1)));
    [qv ii_v]=min(abs(temp_v-p_e1(2)));
    u=temp_u(ii_u);
    v=temp_v(ii_v);
end
ang_1_new(monte_loop,:)=[asin(sqrt(u^2+v^2))*180/pi,atan(v/u)*180/pi];
        
% the second signal
z2=A(:,2)*s(2,:)+(10^(-(SNR)/20))*v_noise(6*m,snap_num)/sqrt(2);
R=1/snap_num*z2*z2';
[es bb cc]=svd(R);

Es=es(:,1);
Es_1=v_j(1,m)*Es;
Es_2=v_j(2,m)*Es;
Es_3=v_j(3,m)*Es;
Es_4=v_j(4,m)*Es;
Es_5=v_j(5,m)*Es;
Es_6=v_j(6,m)*Es;

%获得5个旋转不变量
w1_2=inv(Es_1'*Es_1)*(Es_1'*Es_2);
w2_3=inv(Es_2'*Es_2)*(Es_2'*Es_3);
w3_4=inv(Es_3'*Es_3)*(Es_3'*Es_4);
w4_5=inv(Es_4'*Es_4)*(Es_4'*Es_5);
w5_6=inv(Es_5'*Es_5)*(Es_5'*Es_6);

%获得阵列导向矢量的五个不变量
[u x1_2]=eig(w1_2);
[u x2_3]=eig(w2_3); 
[u x3_4]=eig(w3_4); 
[u x4_5]=eig(w4_5); 
[u x5_6]=eig(w5_6); 

%获得信号Poynting矢量的估计值
tmp_1_1=[1;x1_2(1,1);x1_2(1,1)*x2_3(1,1)];
tmp_1_2=[x1_2(1,1)*x2_3(1,1)*x3_4(1,1);x1_2(1,1)*x2_3(1,1)*x3_4(1,1)*x4_5(1,1);x1_2(1,1)*x2_3(1,1)*x3_4(1,1)*x4_5(1,1)*x5_6(1,1)];
p_e1=cross(tmp_1_1,conj(tmp_1_2))/sqrt(tmp_1_1'*tmp_1_1)/sqrt(tmp_1_2'*tmp_1_2);
%仰俯角和方位角,1为仰俯角，2为方位角
ang(1,1)=(acos(p_e1(3))*180/pi);
ang(1,2)=atan((p_e1(2)/p_e1(1)))*180/pi;
h1=[cos(ang(1,2)*pi/180)*cos(ang(1,1)*pi/180) -sin(ang(1,2)*pi/180)
    sin(ang(1,2)*pi/180)*cos(ang(1,1)*pi/180)  cos(ang(1,2)*pi/180)
   -sin(ang(1,1)*pi/180)                        0
   -sin(ang(1,2)*pi/180)                       -cos(ang(1,2)*pi/180)*cos(ang(1,1)*pi/180)
    cos(ang(1,2)*pi/180)                       -sin(ang(1,2)*pi/180)*cos(ang(1,1)*pi/180)
    0                                           sin(ang(1,1)*pi/180)];
a1=1/sqrt((tmp_1_1'*tmp_1_1)*(tmp_1_2'*tmp_1_2))*[1;x1_2(1,1);x1_2(1,1)*x2_3(1,1);x1_2(1,1)*x2_3(1,1)*x3_4(1,1);x1_2(1,1)*x2_3(1,1)*x3_4(1,1)*x4_5(1,1);x1_2(1,1)*x2_3(1,1)*x3_4(1,1)*x4_5(1,1)*x5_6(1,1)];
G1=inv(h1'*h1)*h1'*a1;
%计算两信号的参数，，3为辅助极化角，4为相位差
ang(1,3)=(atan(abs(G1(1)/G1(2)))*180/pi);
ang(1,4)=((angle(G1(1))-angle(G1(2)))*180/pi);
ang_2(monte_loop,:)=ang;
% new method
z2_new=A_new(:,2)*s(2,:)+(10^(-(SNR)/20))*v_noise(6*m,snap_num)/sqrt(2);
R_new=1/snap_num*z2_new*z2_new';
[aa bb cc]=svd(R_new);
Eo=aa(1:6,1);
Ex=aa(7:12,1);
Ey=aa(13:18,1);

wu=inv(Eo'*Eo)*Eo'*Ex;
wv=inv(Eo'*Eo)*Eo'*Ey;

[qq u]=eig(wu);
[ww v]=eig(wv);

u=angle(u)/(2*pi*d);
v=angle(v)/(2*pi*d);

if d>0.5
    bu=ceil(d*(-1-u)):floor(d*(1-u));
    bv=ceil(d*(-1-v)):floor(d*(1-v));
    temp_u=u+bu/d;
    temp_v=v+bv/d;
    [qu ii_u]=min(abs(temp_u-p_e1(1)));
    [qv ii_v]=min(abs(temp_v-p_e1(2)));
    u=temp_u(ii_u);
    v=temp_v(ii_v);
end
ang_2_new(monte_loop,:)=[asin(sqrt(u^2+v^2))*180/pi,atan(v/u)*180/pi];
end %end of monte_loop
ang_1_kt_ele((d+1)/2,(SNR+15)/5)=v_std(ang_1(:,1),monte_loop);
ang_1_kt_azi((d+1)/2,(SNR+15)/5)=v_std(ang_1(:,2),monte_loop);

ang_2_kt_ele((d+1)/2,(SNR+15)/5)=v_std(ang_2(:,1),monte_loop);
ang_2_kt_azi((d+1)/2,(SNR+15)/5)=v_std(ang_2(:,2),monte_loop);

%v

ang_1_v_ele((d+1)/2,(SNR+15)/5)=v_std(ang_1_new(:,1),monte_loop);
ang_1_v_azi((d+1)/2,(SNR+15)/5)=v_std(ang_1_new(:,2),monte_loop);

ang_2_v_ele((d+1)/2,(SNR+15)/5)=v_std(ang_2_new(:,1),monte_loop);
ang_2_v_azi((d+1)/2,(SNR+15)/5)=v_std(ang_2_new(:,2),monte_loop);
end %end of SNR
end % end of d

SNR=-10:5:30

figure
semilogy(SNR,ang_1_kt_ele(1,:),'-*');
hold on ;
semilogy(SNR,ang_1_kt_azi(1,:),'-square');
hold on ;

semilogy(SNR,ang_1_v_ele(1,:),'-diamond');
hold on ;
semilogy(SNR,ang_1_v_azi(1,:),'->');
xlabel('信噪比(dB)');
ylabel('标准差(deg)');
title('入射信号一空间到达角估计值的标准差与信噪比关系');
hold on ;
%
figure
semilogy(SNR,ang_2_kt_ele(1,:),'-*');
hold on ;
semilogy(SNR,ang_2_kt_azi(1,:),'-square');
hold on ;

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

%%%
figure
semilogy(SNR,ang_1_v_ele(2,:),'-*');
hold on;
semilogy(SNR,ang_1_v_ele(3,:),'-square');
hold on;
semilogy(SNR,ang_1_v_ele(4,:),'-diamond');
hold on;

semilogy(SNR,ang_2_v_ele(2,:),'->');
hold on;
semilogy(SNR,ang_2_v_ele(3,:),'-x');
hold on;
semilogy(SNR,ang_2_v_ele(4,:),'-^');
xlabel('信噪比(dB)');
ylabel('标准差(deg)');
title('不同阵元间距情形下俯仰角估计值的标准差与信噪比关系');
hold on;

figure
semilogy(SNR,ang_1_v_azi(2,:),'-*');
hold on;
semilogy(SNR,ang_1_v_azi(3,:),'-square');
hold on;
semilogy(SNR,ang_1_v_azi(4,:),'-diamond');
hold on;

semilogy(SNR,ang_2_v_azi(2,:),'->');
hold on;
semilogy(SNR,ang_2_v_azi(3,:),'-x');
hold on;
semilogy(SNR,ang_2_v_azi(4,:),'-^');
xlabel('信噪比(dB)');
ylabel('标准差(deg)');
title('不同阵元间距情形下方位角估计值的标准差与信噪比关系');
hold on;
💿 文件大小 12 K
👤 上传用户 foreverNewcomer
📂 所属分类人工智能/神经网络
🏷️ 相关标签

#DOA #极化 #参数估计 #阵列信号处理
⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -