demo6_3_dimult.m

《最佳阵列处理》一书中第六章的例程

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% demo6_3_DImult.m
% K. Bell 12/2/99
% updated by K. Bell 11/20/00
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Multiple Directional Interference
% Example 6.3.1, Fig 6.23, 6.24

%************************
% Array
%************************
clear all
N = 10;                                 % Elements in array
d = 0.5;                                % sensor spacing half wavelength wrt wc
D = [-(N-1)/2:1:(N-1)/2].';
BWNN = 2/(N*d);
u=[-1:0.001:1];
uI = [0.99:-0.01:0.01];
nI = length(uI);
nu=length(u);
vv = exp(j*2*pi*d*D*u);
vI1 = exp(j*2*pi*d*D*uI);
vI2 = exp(-j*2*pi*d*D*uI);

sigma_i = 10^(10/10);

%************************
% Source
%************************
us = 0.0;
AS = exp(j*2*pi*d*D*us);

sigma_n = 1;
sigma_s = 1;
SINRin = sigma_s/(2*sigma_i+sigma_n);

wc = AS/N;
Rs = sigma_s*AS*AS';

T = zeros(1,nI);
Tc = zeros(1,nI);
Tn = zeros(1,nI);
A = zeros(1,nI);
An = zeros(1,nI);
Ac = zeros(1,nI);
Bc = 10*log10(abs(wc'*vv).^2);

for n=1:nI
   
   Snc = sigma_i*vI1(:,n)*vI1(:,n)'+sigma_i*vI2(:,n)*vI2(:,n)';
   Sn = Snc+sigma_n*eye(N);
   V = [vI1(:,n) vI2(:,n)];
   PIo = eye(N)-V*inv(V'*V)*V';
   wn = PIo*AS/real(AS'*PIo*AS);
   Sninv = inv(Sn);
   w = Sninv*AS/real(AS'*Sninv*AS);
   T(n) = 10*log10(real(w'*w));
   Tc(n) = 10*log10(real(wc'*wc));
   Tn(n) = 10*log10(real(wn'*wn));
   I(n) = 10*log10(max(1e-15,real(w'*Snc*w)));
   Ic(n) = 10*log10(max(1e-15,real(wc'*Snc*wc)));
   In(n) = 10*log10(max(1e-15,real(wn'*Snc*wn)));
   A(n) = 10*log10(real((w'*Rs*w)/(w'*Sn*w))/SINRin);
   Ac(n) = 10*log10(real((wc'*Rs*wc)/(wc'*Sn*wc))/SINRin);
   An(n) = 10*log10(real((wn'*Rs*wn)/(wn'*Sn*wn))/SINRin);
   B = 10*log10(abs(w'*vv).^2);
   Bn = 10*log10(abs(wn'*vv).^2);
   
   figure(1)
   h=plot(u,B,'y');
   hold on
   hc=plot(u,Bc,'g');
   hn=plot(u,Bn,'m');
   plot(us*[1 1],[-60 30],'c')
   plot(uI(n)*[1 1],[-60 30],'r')
   plot(-uI(n)*[1 1],[-60 30],'r')
   xlabel('u')
   ylabel('Beampattern (dB)')
   title(['Interference Power ' num2str(10*log10(sigma_i)) ' dB'])
   grid on
   hold off
   legend([hc h hn],'Conventional','MVDR','NUll Steering')
   axis([-1 1 -50 20])
   
   figure(2)
   subplot(3,1,1)
   h=plot(uI(1:n),A([1:n]),'y');
   hold on
   hc=plot(uI(1:n),Ac([1:n]),'g');
   hn=plot(uI(1:n),An([1:n]),'m');
   
   hold off
   axis([0 1 -10 40])
   grid on
   xlabel('u')
   ylabel('Array gain (dB)')
   legend([hc h hn],'Conventional','MVDR','NUll Steering',4)
   
   subplot(3,1,3)
   plot(uI(1:n),T([1:n]),'y')
   hold on
   plot(uI(1:n),Tc([1:n]),'g')
   plot(uI(1:n),Tn([1:n]),'m')
   hold off
   axis([0 1 -20 20])
   grid on
   xlabel('u')
   ylabel('Sensitivity (noise power) (dB)')
   
   subplot(3,1,2)
   plot(uI(1:n),I([1:n]),'y')
   hold on
   plot(uI(1:n),Ic([1:n]),'g')
   plot(uI(1:n),In([1:n]),'m')
   hold off
   axis([0 1 -80 20])
   grid on
   xlabel('u')
   ylabel('Interference power (dB)')
   drawnow
   %pause
end