%已知各信源的到達角度
%假設空間存在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

% 在直角坐標下繪制輸出方向圖
%非自適應的算法，計算量小，但是不能準確的將期望最大化干擾置零
%同時由于主瓣寬度的限制，干擾如果離主瓣很近則完全無法抑制
%改變步長對輸出方向圖沒有影響，只是改變其精度
%增加陣元數可以有效的減小主瓣和旁瓣寬度，提高精確度