v_paper_5.m

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

clear all;
close all;
clc;

j_2p=j*2*pi;
degrad=pi/180;
snapshot=1200;
m=3;


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

A(:,1)=v_a(10,30,15,90);
A(:,2)=v_a(20,40,25,90);
A(:,3)=v_a(30,60,35,90);

p=[v_p(10,30),v_p(20,40),v_p(30,60)];
mis_pos=[0.1,0.2,0.3];

for i=1:3
    mis_q(i)=exp(j*2*pi*(p(1,i)*mis_pos(1)+p(2,i)*mis_pos(2)+p(3,i)*mis_pos(3)));
end
A(4:6,:)=A(4:6,:)*diag(mis_q);


for SNR=0:5:40
%SNR=20;
    for monte_loop=1:200
 %  monte_loop=1
% receive signals

z1=A(:,1)*s(1,:)+10^(-SNR/20)*v_noise(6,snapshot)/sqrt(2);
z2=A(:,2)*s(2,:)+10^(-SNR/20)*v_noise(6,snapshot)/sqrt(2);
z3=A(:,3)*s(3,:)+10^(-SNR/20)*v_noise(6,snapshot)/sqrt(2);

R1=1/snapshot*z1*z1';
R2=1/snapshot*z2*z2';
R3=1/snapshot*z3*z3';
 
[u1 v1]=eig(R1);
[u2 v2]=eig(R2);
[u3 v3]=eig(R3);

Es(:,1)=u1(:,1);
Es(:,2)=u2(:,1);
Es(:,3)=u3(:,1);

for i=1:3
    re_ang(:,i)=abs(angle(cross(Es(1:3,i),conj(Es(4:6,i)))))/(2*pi);
end
for i=1:3
    re_q(i)=mean(re_ang(:,i));
end


%%%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^
% single vector sensor

A2=A*[exp(j_2p*5/25*.20),0,0; 0,exp(j_2p*7/25*.20),0; 0,0,exp(j_2p*9/25*.20)];
%% the first signal
z1=A(:,1)*s(1,:);
z2=A2(:,1)*s(1,:);
z=[z1;z2]+10^(-SNR/20)*v_noise(12,snapshot)/sqrt(2);
Rr=1/snapshot*z*z';
[u v ]=eig(Rr);
a1=u(1:6,1:1);
a2=u(7:12,1:1);
w=inv(a1'*a1)*a1'*a2;
[T WW ]=eig(w);
re_A=a1*inv(T);
mis_p(:,1)=cross(re_A(1:3,1),conj(re_A(4:6,1)))/(norm(re_A(1:3,1))*norm(re_A(4:6,1)));
%%%% the second signal
z1=A(:,2)*s(2,:);
z2=A2(:,2)*s(2,:);
z=[z1;z2]+10^(-SNR/20)*v_noise(12,snapshot)/sqrt(2);
Rr=1/snapshot*z*z';
[u v ]=eig(Rr);
a1=u(1:6,1:1);
a2=u(7:12,1:1);
w=inv(a1'*a1)*a1'*a2;
[T WW ]=eig(w);
re_A=a1*inv(T);
mis_p(:,2)=cross(re_A(1:3,1),conj(re_A(4:6,1)))/(norm(re_A(1:3,1))*norm(re_A(4:6,1)));
% the third signal 
z1=A(:,3)*s(3,:);
z2=A2(:,3)*s(3,:);
z=[z1;z2]+10^(-SNR/20)*v_noise(12,snapshot)/sqrt(2);
Rr=1/snapshot*z*z';
[u v ]=eig(Rr);
a1=u(1:6,1:1);
a2=u(7:12,1:1);
w=inv(a1'*a1)*a1'*a2;
[T WW ]=eig(w);
re_A=a1*inv(T);
mis_p(:,3)=cross(re_A(1:3,1),conj(re_A(4:6,1)))/(norm(re_A(1:3,1))*norm(re_A(4:6,1)));
%%%%%%%
p_2=abs(mis_p);
%%%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

% compute dz
%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%8888888888888
% the first method
c1=p(1,1)* (p(2,2)*re_q(1)-p(2,1)*re_q(2)) / (p(2,2)*p(1,1)-p(2,1)*p(1,2));
c2=p(2,1)* (p(1,3)*re_q(2)-p(1,2)*re_q(3))/  (p(1,3)*p(2,2)-p(1,2)*p(2,3));
c3=p(1,1)* (p(2,2)*p(3,1)-p(2,1)*p(3,2))/    (p(2,2)*p(1,1)-p(2,1)*p(1,2));
c4=p(2,1)* (p(1,3)*p(3,2)-p(1,2)*p(3,3))/    (p(1,3)*p(2,2)-p(1,2)*p(2,3));
dz_1(monte_loop)=(c1+c2-re_q(1))/(c3+c4-p(3,1));

cy1=(p(1,3)*p(3,2)-p(1,2)*p(3,3))*(re_q(1)*(p(3,2)*p(1,1)-p(3,1)*p(1,2))-p(1,1)*(p(3,2)*re_q(1)-p(3,1)*re_q(2)));
cy2=p(3,1)*(p(3,2)*p(1,1)-p(3,1)*p(1,2))*(p(1,3)*re_q(2)-p(1,2)*re_q(3));
cy3=(p(1,3)*p(3,2)-p(1,2)*p(3,3))*(p(2,1)*(p(3,2)*p(1,1)-p(3,1)*p(1,2))-p(1,1)*(p(3,2)*p(2,1)-p(3,1)*p(2,2)));
cy4=p(3,1)*(p(3,2)*p(1,1)-p(3,1)*p(1,2))*(p(1,3)*p(2,2)-p(1,2)*p(2,3));
dy_1(monte_loop)=(cy1-cy2)/(cy3-cy4)

cz1=(p(2,3)*p(3,2)-p(2,2)*p(3,3))*(   re_q(1)*(p(3,2)*p(2,1)-p(3,1)*p(2,2))  -p(2,1)*(p(3,2)*re_q(1)-p(3,1)*re_q(2))   );
cz2=p(3,1)*(p(3,2)*p(2,1)-p(3,1)*p(2,2))  *(p(2,3)*re_q(2)-p(2,2)*re_q(3));
cz3=(p(3,2)*p(2,1)-p(3,1)*p(2,2))*   (  p(1,1)*(p(2,3)*p(3,2)-p(2,2)*p(3,3))  -p(3,1)*(p(2,3)*p(1,2)-p(2,2)*p(1,3)));
cz4=p(2,1)*(p(2,3)*p(3,2)-p(2,2)*p(3,3))*(p(3,2)*p(1,1)-p(3,1)*p(1,2));
dx_1(monte_loop)=(cz1-cz2)/(cz3-cz4)

%the second method

% compute dz
c1=p_2(1,1)* (p_2(2,2)*re_q(1)- p_2(2,1)*re_q(2)) /    (p_2(2,2)*p_2(1,1)-p_2(2,1)*p_2(1,2));
c2=p_2(2,1)* (p_2(1,3)*re_q(2)- p_2(1,2)*re_q(3))/     (p_2(1,3)*p_2(2,2)-p_2(1,2)*p_2(2,3));
c3=p_2(1,1)* (p_2(2,2)*p_2(3,1)-p_2(2,1)*p_2(3,2))/    (p_2(2,2)*p_2(1,1)-p_2(2,1)*p_2(1,2));
c4=p_2(2,1)* (p_2(1,3)*p_2(3,2)-p_2(1,2)*p_2(3,3))/    (p_2(1,3)*p_2(2,2)-p_2(1,2)*p_2(2,3));
dz_2(monte_loop)=(c1+c2-re_q(1))/(c3+c4-p_2(3,1))

cy1=(p_2(1,3)*p_2(3,2)-p_2(1,2)*p_2(3,3))*(re_q(1)*(p_2(3,2)*p_2(1,1)-p_2(3,1)*p_2(1,2))-p_2(1,1)*(p_2(3,2)*re_q(1)-p_2(3,1)*re_q(2)));
cy2=p_2(3,1)*(p_2(3,2)*p_2(1,1)-p_2(3,1)*p_2(1,2))*(p_2(1,3)*re_q(2)-p_2(1,2)*re_q(3));
cy3=(p_2(1,3)*p_2(3,2)-p_2(1,2)*p_2(3,3))*(p_2(2,1)*(p_2(3,2)*p_2(1,1)-p_2(3,1)*p_2(1,2))-p_2(1,1)*(p_2(3,2)*p_2(2,1)-p_2(3,1)*p_2(2,2)));
cy4=p_2(3,1)*(p_2(3,2)*p_2(1,1)-p_2(3,1)*p_2(1,2))*(p_2(1,3)*p_2(2,2)-p_2(1,2)*p_2(2,3));
dy_2(monte_loop)=(cy1-cy2)/(cy3-cy4)


cz1=(p_2(2,3)*p_2(3,2)-p_2(2,2)*p_2(3,3))*(   re_q(1)*(p_2(3,2)*p_2(2,1)-p_2(3,1)*p_2(2,2))  -p_2(2,1)*(p_2(3,2)*re_q(1)-p_2(3,1)*re_q(2))   );
cz2=p_2(3,1)*(p_2(3,2)*p_2(2,1)-p_2(3,1)*p_2(2,2))  *(p_2(2,3)*re_q(2)-p_2(2,2)*re_q(3));
cz3=(p_2(3,2)*p_2(2,1)-p_2(3,1)*p_2(2,2))*   (  p_2(1,1)*(p_2(2,3)*p_2(3,2)-p_2(2,2)*p_2(3,3))  -p_2(3,1)*(p_2(2,3)*p_2(1,2)-p_2(2,2)*p_2(1,3)));
cz4=p_2(2,1)*(p_2(2,3)*p_2(3,2)-p_2(2,2)*p_2(3,3))*(p_2(3,2)*p_2(1,1)-p_2(3,1)*p_2(1,2));
dx_2(monte_loop)=(cz1-cz2)/(cz3-cz4)
end
bias_x_1((SNR+5)/5)=max(abs((dx_1)-0.1));
bias_y_1((SNR+5)/5)=max(abs((dy_1)-0.2));
bias_z_1((SNR+5)/5)=max(abs((dz_1)-0.3));

x_1((SNR+5)/5)=v_std(dx_1,monte_loop);
y_1((SNR+5)/5)=v_std(dy_1,monte_loop);
z_1((SNR+5)/5)=v_std(dz_1,monte_loop);

bias_x_2((SNR+5)/5)=max(abs((dx_2)-0.1));
bias_y_2((SNR+5)/5)=max(abs((dy_2)-0.2));
bias_z_2((SNR+5)/5)=max(abs((dz_2)-0.3));

x_2((SNR+5)/5)=v_std(dx_2,monte_loop);
y_2((SNR+5)/5)=v_std(dy_2,monte_loop);
z_2((SNR+5)/5)=v_std(dz_2,monte_loop);
end


SNR=0:5:40
figure
plot(SNR,x_1,'-*');
hold on;
plot(SNR,x_2,'-s');
xlabel('信噪比(dB)');
ylabel('标准偏差');
%title('X轴向位置偏差与信噪比关系');


figure
plot(SNR,y_1,'-*');
hold on;
plot(SNR,y_2,'-s');
xlabel('信噪比(dB)');
ylabel('标准偏差');
%title('Y轴向位置偏差与信噪比关系');

figure
plot(SNR,z_1,'-*');
hold on;
plot(SNR,z_2,'-s');
xlabel('信噪比(dB)');
ylabel('标准偏差');
%title('Z轴向位置偏差与信噪比关系');