fig6_70.m

《最优阵列处理》一书第六章的MATLAB例程

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 6.70
% Lillian Xu 12/18/2000
% updated by K. Bell 1/18/08
% function called: sinc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
close all

N = 10;
BWNN = 4/N;
n = (-(N-1)/2:(N-1)/2)';
ud = 0.1;
signalRange = [-ud:1/500:ud];
INRrange = [10 30];
SNRrange = [-10:2:40];
LNRrange = [-20:2:50];

ui1 = -0.30;
ui2 = 0.30;

Vi1 = exp(j*n*pi*ui1);
Vi2 = exp(j*n*pi*ui2);

% MPDR
C1 = exp(j*n*pi*[0]);
f1 = [1];

% directional
u1 = 0.0866;
uc = [0, -u1, u1];
C2 = exp(j*n*pi*uc);
f2 = [1;sin(-(N/2)*pi*u1)./(N*sin(-.5*pi*u1));sin((N/2)*pi*u1)./(N*sin(.5*pi*u1))];

% derivative
coef = 0.8;
C3 = [ones(N,1),j*n,-n.^2];
f3 = [1;0;coef*(1-N^2)/12];

% Taylor
SLL=20;
nBar=2;
Ne = 2;

scale = N*0.5;
A = 1/pi*acosh(10^(SLL/20));
uroot=[1:1:floor((N-1)/2)]/scale;     % (N-1)/2 for odd and N/2 - 1 for N even
for m = 1:nBar-1
    Vn = nBar*sqrt( (A^2+(m-0.5)^2)/(A^2+(nBar-0.5)^2) );
    uroot(m) = Vn/scale;
end
if rem(N,2) == 0 % even
    roots = [0 uroot -uroot 1];
else
    roots = [0 uroot -uroot];
end

Vr = exp(j*n*pi*roots);
w0 = inv(Vr')*[1;zeros(N-1,1)];

wdBar = w0/norm(w0)^2;
PcWdBar = eye(N) - wdBar*inv(wdBar'*wdBar)*wdBar';
for h1 = 1:N
    for h2 = 1:N
        Rs(h1,h2) = 2*ud*sinc( (h1-h2)*ud );
    end
end
RsBar = PcWdBar*Rs*PcWdBar;
[U,lambda,V] = svd(RsBar);
Cs = V(:,1:Ne);					%select 3 eigenvectors
C4 = [wdBar,Cs];
f4 = [1;zeros(Ne,1)];

Gain1 = zeros(length(INRrange),length(SNRrange),length(LNRrange));
Gain2 = zeros(length(INRrange),length(SNRrange),length(LNRrange));
Gain3 = zeros(length(INRrange),length(SNRrange),length(LNRrange));
Gain4 = zeros(length(INRrange),length(SNRrange),length(LNRrange));

k1 = 1;
for INR = 10.^(INRrange/10)
    disp(['loop ' int2str(k1) ' of 2 ...'])
    k2 = 1;
    for SNR = 10.^(SNRrange/10)
        k3 = 1;
        SINRi = SNR/(1+2*INR);
        for LNR = 10.^(LNRrange/10)
            for ua = signalRange
                Vs = exp(j*n*pi*ua);
                Ss = SNR*Vs*Vs';
                Sn = eye(N) + INR*Vi1*Vi1'+ INR*Vi2*Vi2';
                Sx = Ss + Sn + LNR*eye(N);     % Diagonal Loading
                % MPDR
                W = inv(Sx)*C1*inv(C1'*inv(Sx)*C1)*f1;           %consider White Noise only
                SINR0 = real((W'*Ss*W)/(W'*Sn*W));
                Gain1(k1,k2,k3) = Gain1(k1,k2,k3)+SINR0/SINRi;
                
                % Directional
                W = inv(Sx)*C2*inv(C2'*inv(Sx)*C2)*f2;           %consider White Noise only
                SINR0 = real((W'*Ss*W)/(W'*Sn*W));
                Gain2(k1,k2,k3) = Gain2(k1,k2,k3)+SINR0/SINRi;
                
                % Derivative
                W = inv(Sx)*C3*inv(C3'*inv(Sx)*C3)*f3;           %consider White Noise only
                SINR0 = real((W'*Ss*W)/(W'*Sn*W));
                Gain3(k1,k2,k3) = Gain3(k1,k2,k3)+SINR0/SINRi;
                
                % Taylor
                W = inv(Sx)*C4*inv(C4'*inv(Sx)*C4)*f4;           %consider White Noise only
                SINR0 = real((W'*Ss*W)/(W'*Sn*W));
                Gain4(k1,k2,k3) = Gain4(k1,k2,k3)+SINR0/SINRi;
            end  % end of ua
            k3 = k3 + 1;
        end	% end of LNR
        k2 = k2 + 1;
    end		% end of SNR
    k1 = k1 + 1;
end			% end of INR

Gain1 = 10*log10( abs(Gain1)/length(signalRange) ) ;
Gain2 = 10*log10( abs(Gain2)/length(signalRange) ) ;
Gain3 = 10*log10( abs(Gain3)/length(signalRange) ) ;
Gain4 = 10*log10( abs(Gain4)/length(signalRange) ) ;

for k1 = 1:size(Gain1,1)
    for k2 = 1:size(Gain1,2)
        [A1(k1,k2),I(k1,k2)] = max(Gain1(k1,k2,:));
        [A2(k1,k2),I(k1,k2)] = max(Gain2(k1,k2,:));
        [A3(k1,k2),I(k1,k2)] = max(Gain3(k1,k2,:));
        [A4(k1,k2),I(k1,k2)] = max(Gain4(k1,k2,:));
    end
end

figure
plot(SNRrange,A1(1,:),'-',SNRrange,A2(1,:),'--',...
   SNRrange,A3(1,:),'-.',SNRrange,A4(1,:),':');

hold on 
xlabel('{\it SNR} (dB)')
ylabel('Optimal gain (dB)')
legend('MPDR','DIR,u_c=[0,0.0866,-0.0866],g=[1;B_C;B_C]',...
   'DER,C=[1,d(0),d(0)],\alpha=0.8','QP,Taylor(-20dB,{\it n}=3),Ne=2','location','southwest')
%title('performance comparison with optimal LNR,ui=-0.3 & +0.3, INR=10dB(each), N=10')
axis([-10 40 5 25])
grid on

figure
plot(SNRrange,A1(2,:),'-',SNRrange,A2(2,:),'--',...
   SNRrange,A3(2,:),'-.',SNRrange,A4(2,:),':');

hold on 
xlabel('{\it SNR} (dB)')
ylabel('Optimal gain (dB)')
legend('MPDR','DIR,u_c=[0,0.0866,-0.0866],g=[1;B_C;B_C]',...
   'DER,C=[1,d(0),d(0)],\alpha=0.8','QP,Taylor(-20dB,{\it n}=3),Ne=2','location','southwest')
%title('performance comparison with optimal LNR,ui=-0.3 & +0.3, INR=30dB(each), N=10')
axis([-10 40 20 45])
grid on