clear,clc
m=8;                                      % sensors
n=2;                                      % sources
theta=[-20 0];                            % in angle
d=1/2;                                    % 1/2 lambada
N=100;                                    % samples
L=100;                                    % resolution in [-90' 90']
meann=0;                                  % mean of noise
varn=1;                                   % variance of noise
SNR=10;                                  
INR=20;                                  
rvar1=sqrt(varn) * 10^(SNR/20);           % variance of signal
rvar2=sqrt(varn) * 10^(INR/20);           % variance of interference

% generate the source signals
s=[rvar1*exp(j*2*pi*(50*0.001*[0:N-1]+rand))
   rvar2*exp(j*2*pi*(100*0.001*[0:N-1]+rand))];
% generate the A matrix
A=exp(-j*2*pi*d*[0:m-1].'*sin(theta*pi/180));
% generate the noise component
e=sqrt(varn/2)*(randn(m,N)+j*randn(m,N));
% generate the ULA data
Y=A*s+e;
% compute R
R=Y*Y'/N;
% inverse of R
IR=inv(R);

% beamform using the MVDR method
a1=A(:,1);
w=IR*a1/real(a1'*IR*a1);
beam=zeros(1,L);
for i = 1 : L
   a=exp(-j*2*pi*d*[0:m-1].'*sin(-pi/2 + pi*(i-1)/L));
   beam(i)=20*log10(abs(w'*a));
end
angle=-90:180/L:(90-180/L);
plot(angle,beam,'b');
hold on;

%*****************************************************************
m=8;  
n=2;  
theta=[-20 0]; 
d=1/2;         
N=100; 
L=100; 
meann=0; 
varn=1;  
SNR=10;
INR=30; 
rvar1=sqrt(varn) * 10^(SNR/20); 
rvar2=sqrt(varn) * 10^(INR/20); 

s=[rvar1*exp(j*2*pi*(50*0.001*[0:N-1]+rand))
   rvar2*exp(j*2*pi*(100*0.001*[0:N-1]+rand))];
A=exp(-j*2*pi*d*[0:m-1].'*sin(theta*pi/180));
e=sqrt(varn/2)*(randn(m,N)+j*randn(m,N));
Y=A*s+e;
R=Y*Y'/N;
IR=inv(R);
% beamform using the MVDR method
a1=A(:,1);
w=IR*a1/real(a1'*IR*a1);
beam=zeros(1,L);
for i = 1 : L
   a=exp(-j*2*pi*d*[0:m-1].'*sin(-pi/2 + pi*(i-1)/L));
   beam(i)=20*log10(abs(w'*a));
end
angle=-90:180/L:(90-180/L);
plot(angle,beam,'r');

%****************************************************************
m=8;  
n=2;  
theta=[-20 0]; 
d=1/2;         
N=100; 
L=100; 
meann=0; 
varn=1;  
SNR=10; 
INR=40;
rvar1=sqrt(varn) * 10^(SNR/20);
rvar2=sqrt(varn) * 10^(INR/20); 

s=[rvar1*exp(j*2*pi*(50*0.001*[0:N-1]+rand))
   rvar2*exp(j*2*pi*(100*0.001*[0:N-1]+rand))];
A=exp(-j*2*pi*d*[0:m-1].'*sin(theta*pi/180));
e=sqrt(varn/2)*(randn(m,N)+j*randn(m,N));
Y=A*s+e;
R=Y*Y'/N;
IR=inv(R);
% beamform using the MVDR method
a1=A(:,1);
w=IR*a1/real(a1'*IR*a1);
beam=zeros(1,L);
for i = 1 : L
   a=exp(-j*2*pi*d*[0:m-1].'*sin(-pi/2 + pi*(i-1)/L));
   beam(i)=20*log10(abs(w'*a));
end
angle=-90:180/L:(90-180/L);
plot(angle,beam,'g');

xlabel('angle');
ylabel('幅度響應/dB');
gtext({'\bf\leftarrowINR=20';'\bf\leftarrowINR=30';'\bf\leftarrowINR=40'})