ucaesprit.m

本程序主要是对于圆阵列下的角度进行估计

clc;
clear;
sig_num = 1000;%%信号的快拍数
rr=1;
snr=10;     
rou=0;
x1=randint(sig_num,1,4);%产生sig_sum*1的矩阵，数值在0-3之间
x2=randint(sig_num,1,4);
x3=randint(sig_num,1,4);
x222=rou*x1+sqrt(1-(abs(rou))^2)*x3;%产生相关的两路信号
x2=round(x222);
for i=1:sig_num
    if x2(i)>=4 
        x2(i)=3;
    elseif x2(i)<=0
        x2(i)=0;
    end
end
y1=dmod(x1,1,0.1,4,'psk',4);
y1=y1(1:sig_num);
y2=dmod(x2,1,0.1,4,'psk',4);
y2=y2(1:sig_num);
x=[y1';y2'];
coeff=corr2(y1,y2);
s=x;
q=2;
p=8;%智能天线阵元数
h=3;%最大模式数
P=2*h+1;
fc=200;%*(10^6); %hz
lemda=3*10^8/fc;
banjing=lemda/(4*sin(pi/p));
jh=pi./180;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%均匀圆阵转化为均匀线阵
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
jh=pi/180;
alpha=[10,30];
theta=[20,40];
kesei=2*pi*banjing./lemda;%%
j=sqrt(-1);
a=[];
n=[0:(p-1)]';
gama=2*pi*n/p;
for i=1:q
    for n=0:(p-1)
        a(n+1,i)=exp(j*kesei*sin(theta(i)*jh)*cos(alpha(i)*jh-gama(n+1)));
    end
end
x11=a*s;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c1=zeros(p,sig_num);
vn=var(x11.')./(10^(snr/10));
n11=randn(p,sig_num)+j*randn(p,sig_num);
a1=(vn./var(n11.'));
b1=((n11.'-ones(size(n11,2),1)*mean(n11.'))).';
for i=1:p
    c1(i,:)=a1(i)*b1(i,:);
end
xn=x11+c1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=-h:h
    cv1(i+h+1)=j^(-abs(i));
end
cv=diag(cv1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=0:p-1
    for m=-h:h
        w(i+1,m+h+1)=exp(-j*2*pi*m*i./p);
    end
end
V=1./sqrt(p)*w;
Fe=cv*V';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=-h:h
    alph(i+h+1)=2*pi*i/P;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=-h:h
    for m=-h:h
        v(m+h+1,i+h+1)=exp(j*2*pi*i*m./P);
    end
end
W=(1/sqrt(P))*v;
Fr=W'*Fe;
y=Fr*xn;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%music算法实现过程
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r1=xn*xn'./sig_num;
ry=Fr*r1*Fr';
r=real(ry);
[u,sss]=eig(r);
signalspace=u(:,1:q);
noisespace=u(:,((q+1):P));
ii=1;
for i=h:(-1):-h
    if i<1
        coo(ii)=1;
    else
        coo(ii)=(-1)^i;
    end
    ii=ii+1;
end
co=diag(coo);
so=co*W*signalspace;
me=P-2;
s1=so(1:me,:);
s2=so(2:(me+1),:);
s3=so(3:(me+2),:);
ii=1;
for i=-(h-2):h
    dd(ii)=(-1)^i;    
    ii=ii+1;
end
D=diag(dd);
ini1=eye(me);
for i=1:me
    for m=1:me
        ini(i,me-m+1)=ini1(i,m);
    end
end
E=[s1 D*ini*conj(s1)];
gama=zeros(2*(h-1)+1,1);
ii=1;
for i=-(h-1):h-1
    gama(ii)=i;
    ii=ii+1;
end
gam=(lemda/(pi*banjing))*(diag(gama));
ks=inv(E'*E)*E'*gam*s2;
leng=length(ks(:,1));
kiss=ks((1:(leng/2)),:);
[U,S]=eig(kiss);
miu=diag(S);
for i=1:size(miu)
    thetafind(i)=(asin(abs(miu(i))))/jh;
    alphafind(i)=angle(miu(i))/jh;
end
alphafind
thetafind