simplifiedrbmusic.m

本程序对于圆阵列的DOA估计空间平滑算法

clc;
clear;
sig_num = 1000;%%信号的快拍数
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'];
s=x;
q=2;
p=8;%智能天线阵元数
h=3;%最大模式数
P=2*h+1;
fc=2000*10^6;%*(10^6); %hz
lemda=3*10^8/fc;
banjing=lemda/(4*sin(pi/p));
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;
snr=10;
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
xnc=x11+c1;
xn=xnc;
for i=0:p./2-1
   xn(i+p./2+1,:)=xnc(p-i,:);
end

for i=0:p-1
    for m=-h:h
        v(i+1,m+h+1)=exp(j*2*pi*i*m./P);
    end
end
W=(1/sqrt(8))*v;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%music算法实现过程
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fr=W';
rx=xn*xn'./sig_num;
ry=Fr*rx*Fr';
r=real(ry);
[eigenvector,D]=eig(r,'nobalance');         
R_eigenvalue=diag(abs(D))';    
[eigenvalue,index]=sort(R_eigenvalue);     
R_eigenvalue_decreasing=fliplr(eigenvalue);
eigenvector=eigenvector(:,index);
noisespace1=fliplr(eigenvector);
noisespace=noisespace1(:,q+1:p-1);
bbb=1;
for alpha1=0:100
    kk=1;
    for theta1=0:90
    for n=0:p./2-1
        aa(n+1)=exp(j*kesei*sin(theta1*jh)*cos(alpha1*jh-gama(n+1)));
    end
    for i=0:p./2-1
        aa(p-i)=conj(aa(i+1));
    end
    bb=Fr*aa.';
    ww=(bb.')*noisespace*(noisespace.')*bb;
    smusic(kk,bbb)=abs(1./ww);
    kk=kk+1;
    end
    bbb=bbb+1;
end
alpha=0:100;
theta=0:90;
figure(1);
contour(alpha,theta,10*log10(smusic),2,'kd');
title('MUSIC功率谱');
xlabel('方位角 单位：度'); 
ylabel('俯仰角 单位：度');
grid on;