uniformestimatesdoaparkarray.doc

均匀园阵doa估计music算法源程序

clc;
clear all;
Pd=2000;
Fd=1;
Fs=4*Fd;
R=0.5;
Delay=5;
M=4;
x1=randint(Pd,1,M); 
x2=randint(Pd,1,M);
x3=randint(Pd,1,M);
y1=modmap(x1,Fd,Fs,'qask',M);
y2=modmap(x2,Fd,Fs,'qask',M);
y3=modmap(x3,Fd,Fs,'qask',M);
[rcv_a1,ti]=rcosflt(y1,Fd,Fs,'fir/sqrt/Fs',R,Delay);
[rcv_a2,ti]=rcosflt(y2,Fd,Fs,'fir/sqrt/Fs',R,Delay);
[rcv_a3,ti]=rcosflt(y3,Fd,Fs,'fir/sqrt/Fs',R,Delay);
s1=amodce(rcv_a1,10,'qam'); 
s2=amodce(rcv_a2,10,'qam');
s3=amodce(rcv_a3,10,'qam');
save sig3 s1 s2 s3;
m=16;
p=3;
%俯仰角为入射方向与圆阵轴线之间的夹角,方位角为x轴到入射线在圆阵上投影的夹角
st1=10; % 俯仰角
dt1=10; % 方位角
st2=40;
dt2=40;
st3=70;
dt3=70;
st=[st1;st2;st3];
dt=[dt1;dt2;dt3];
nn=1024; 
SN1=20; 
SN2=20;
SN3=20;
sn=[SN1;SN2;SN3];
degrad=pi/180;
load sig3
tt=1:nn;
S=[s1(tt),s2(tt),s3(tt)].';
% noise=randn(m,nn)+j*randn(m,nn); 
% noise=randn(m+1,nn)+j*randn(m+1,nn); 
Ps=S*S'/nn;
ps=diag(Ps);
refp=2*10.^(sn/10);
tmp=sqrt(refp./ps);
S2=diag(tmp)*S;
% 圆半径为半波长
A=zeros(m,p);
% A=zeros(m+1,p);
% A[0:0,:]=[1 1 1];
% k=[0:m]';
k=[0:m-1]';
for t=1:p
A(:,t)=exp(j*pi*(sin(st(t)*degrad)*cos(dt(t)*degrad)*cos(2*pi*k/m)+sin(st(t)*degrad)*sin(dt(t)*degrad)*sin(2*pi*k/m)));
end
AA=ones(m+1,p);
AA(2:m+1,=A;
% X=AA*S2+noise;
X=AA*S2;
Rxx=X*X'/nn;
RXX=zeros(9,9);
for rr=1:9
? ? RXX(1,rr)=Rxx(1,rr);
end
for tt=2:9
? ? RXX(tt,1)=Rxx(1,tt+8);
end
for rr=1:9
? ? for tt=1:9
? ?? ???if (tt-rr==0)
? ?? ?? ?? ?RXXX(rr,tt)=RXX(1,1);
? ?? ???elseif(tt-rr==1)
? ?? ?? ?? ?RXXX(rr,tt)=RXX(1,2);
? ?? ???elseif(tt-rr==2)
? ?? ?? ?? ?RXXX(rr,tt)=RXX(1,3);? ?
? ?? ???elseif(tt-rr==3)
? ?? ?? ?? ?RXXX(rr,tt)=RXX(1,4);??
? ?? ???elseif(tt-rr==4)
? ?? ?? ?? ?RXXX(rr,tt)=RXX(1,5); 
? ?? ???elseif (tt-rr==5)
? ?? ?? ?? ?RXXX(rr,tt)=RXX(1,6);
? ?? ???elseif(tt-rr==6)
? ?? ?? ?? ?RXXX(rr,tt)=RXX(1,7);
? ?? ???elseif(tt-rr==7)
? ?? ?? ?? ?RXXX(rr,tt)=RXX(1,8);? ?
? ?? ???elseif(tt-rr==8)
? ?? ?? ?? ?RXXX(rr,tt)=RXX(1,9);??
? ?? ? elseif(tt-rr==-1)
? ?? ?? ?? ?RXXX(rr,tt)=RXX(2,1);
? ?? ???elseif(tt-rr==-2)
? ?? ?? ?? ?RXXX(rr,tt)=RXX(3,1);? ?
? ?? ? elseif(tt-rr==-3)
? ?? ?? ?? ?RXXX(rr,tt)=RXX(4,1);??
? ?? ???elseif(tt-rr==-4) 
? ?? ?? ?? ?RXXX(rr,tt)=RXX(5,1);? ???
? ?? ???elseif(tt-rr==-5)
? ?? ?? ?? ?RXXX(rr,tt)=RXX(6,1);
? ?? ???elseif(tt-rr==-6)
? ?? ?? ?? ?RXXX(rr,tt)=RXX(7,1);? ?
? ?? ? elseif(tt-rr==-7)
? ?? ?? ?? ?RXXX(rr,tt)=RXX(8,1);??
? ?? ???else 
? ?? ?? ?? ?RXXX(rr,tt)=RXX(9,1);??
? ?? ???end
? ? end
end
[U,s,v]=svd(RXXX);
Vs=U(:,1:p);
Vu=U(:,p+1:9);
for st1=0:90
for dt1=0:90
k=[0:7]';
AAA=exp(j*pi*(sin(st1*degrad)*cos(dt1*degrad)*cos(2*pi*k/16)+sin(st1*degrad)*sin(dt1*degrad)*sin(2*pi*k/16)));
AAAA=[1 AAA.'];
AAAA=AAAA.';
WW=AAAA'*Vu*Vu'*AAAA;
Pmusic(st1+1,dt1+1)=abs(25./WW);%角谱
end 
end 
sst=0:90;
ddt=0:90;
mesh(sst,ddt,Pmusic);
title('MUSIC空间谱');
xlabel('方位角'); 
ylabel('俯仰角');
zlabel('角谱');
grid on;