music32.m

32阵元的MUSIC算法

%已知各信源的到达角度
%假设空间存在3个信号源，0为期望信号，1，2为干扰，各信号不相关；各信号功率相等
%采用32个阵元的线阵

sita0=50*pi/180;sita1=10*pi/180;sita2=30*pi/180;;   %各信号DOA
amp1=1;amp2=1;amp3=1;          %各信号幅度

t1=pi*sin(sita1);
a_sita1=[1,exp(i*t1),exp(i*2*t1),exp(i*3*t1),exp(i*4*t1),...
         exp(i*5*t1),exp(i*6*t1),exp(i*7*t1),exp(i*8*t1),exp(i*9*t1),exp(i*10*t1),...
         exp(i*11*t1),exp(i*12*t1),exp(i*13*t1),exp(i*14*t1),exp(i*15*t1),exp(i*16*t1),exp(i*17*t1),exp(i*18*t1),exp(i*19*t1),...
         exp(i*20*t1),exp(i*21*t1),exp(i*22*t1),exp(i*23*t1),exp(i*24*t1),exp(i*25*t1),...
         exp(i*26*t1),exp(i*27*t1),exp(i*28*t1),exp(i*29*t1),exp(i*30*t1),exp(i*31*t1)]';
t0=pi*sin(sita0);
a_sita0=[1,exp(i*t0),exp(i*2*t0),exp(i*3*t0),exp(i*4*t0),...
         exp(i*5*t0),exp(i*6*t0),exp(i*7*t0),exp(i*8*t0),exp(i*9*t0),exp(i*10*t0),...
         exp(i*11*t0),exp(i*12*t0),exp(i*13*t0),exp(i*14*t0),exp(i*15*t0),exp(i*16*t0),exp(i*17*t0),exp(i*18*t0),exp(i*19*t0),...
         exp(i*20*t0),exp(i*21*t0),exp(i*22*t0),exp(i*23*t0),exp(i*24*t0),exp(i*25*t0),...
         exp(i*26*t0),exp(i*27*t0),exp(i*28*t0),exp(i*29*t0),exp(i*30*t0),exp(i*31*t0)]';
t2=pi*sin(sita2);
a_sita2=[1,exp(i*t2),exp(i*2*t2),exp(i*3*t2),exp(i*4*t2),...
         exp(i*5*t2),exp(i*6*t2),exp(i*7*t2),exp(i*8*t2),exp(i*9*t2),exp(i*10*t2),...
         exp(i*11*t2),exp(i*12*t2),exp(i*13*t2),exp(i*14*t2),exp(i*15*t2),exp(i*16*t2),exp(i*17*t2),exp(i*18*t2),exp(i*19*t2),...
         exp(i*20*t2),exp(i*21*t2),exp(i*22*t2),exp(i*23*t2),exp(i*24*t2),exp(i*25*t2),...
         exp(i*26*t2),exp(i*27*t2),exp(i*28*t2),exp(i*29*t2),exp(i*30*t2),exp(i*31*t2)]';
% 生成方向矩阵

    N=5000;
    bit1=2*rem(unidrnd(2,1,N),2)-1;
    bit2=2*rem(unidrnd(2,1,N),2)-1;
    bit3=2*rem(unidrnd(2,1,N),2)-1;
 %生成N比特的用户数据 
 
    SNRdB=20;%暂定一个信噪比
    n0=1/(10^(SNRdB/10));
    sigma=sqrt(n0);
    
for nn=1:N;
   for mm=1:32;
      xf(mm,nn)=amp1*exp(i*pi*(mm-1)*sin(sita1))+...
         amp2*exp(i*pi*(mm-1)*sin(sita0))+...
         amp3*exp(i*pi*(mm-1)*sin(sita2))+ sigma*randn(1);
   end;
end;
% 产生观测数据

Cxx=0;
for ri=1:N;
   Cxx=Cxx+xf(:,ri)*xf(:,ri)';
end;
Rxd=Cxx/N;
% 计算协方差矩阵

iRxd=inv(Rxd);
sw1=iRxd*a_sita0;
sw2=a_sita0'*iRxd*a_sita0;
Wopt=sw1/sw2;
%计算波束形成的最佳权向量
% Wopt=R^(-1)a(sitad)/[a'(sitad)*R^(-1)*a(sitad)];  其中sitad为期望信号的DOA;
% '为哈密特转制

angle=0:0.1:90;   %在[0，90]范围内绘制波束图，步长为0.01度

for n=1:length(angle);
     sita=(n-1)*pi/1800;
   for l=1:32
      xd(l)=amp1*bit1(N)*exp(i*pi*(l-1)*sin(sita))+...
         amp2*bit1(N)*exp(i*pi*(l-1)*sin(sita))+...
         amp3*bit1(N)*exp(i*pi*(l-1)*sin(sita))+ sigma*randn(1);
            end;
   Xh=[xd(1),xd(2),xd(3),xd(4),xd(5),xd(6),xd(7),xd(8),xd(9),xd(10),xd(11),xd(12),xd(13),xd(14),xd(15),...
       xd(16),xd(17),xd(18),xd(19),xd(20),xd(21),xd(22),xd(23),xd(24),xd(25),xd(26),xd(27),xd(28),xd(29),xd(30),xd(31),xd(32) ]';
   y(n)=abs(Wopt'*Xh);        %阵列输出
end;
ymax=max(abs(y));
logy=20*log10(abs(y)/ymax);   %归一化，转换为分贝表示

plot(angle,logy);
hold on

% 在直角坐标下绘制输出方向图
%非自适应的算法，计算量小，但是不能准确的将期望最大化干扰置零
%同时由于主瓣宽度的限制，干扰如果离主瓣很近则完全无法抑制
%改变步长对输出方向图没有影响，只是改变其精度
%增加阵元数可以有效的减小主瓣和旁瓣宽度，提高精确度