zuoye23.m

array processing code for communication engineering

clear all;
k=0.5;
N=10;
K=10*N;
m_carlo=50;
dir_signal=0;

count=0;
for i=0:0.05:2
    count=count+1;
    dir_gan(count)=i*0.4;
end

leg=count;

SNR=20;

bartlett=zeros(1,leg);
capon=zeros(1,leg);
music=zeros(1,leg);
rmusic=zeros(1,leg);

for mon=1:m_carlo
    mon
    
    for xx=1:leg
        v = zeros(N,1);
        u = zeros(N,1);
        vrr=zeros(N,1);
        vii=zeros(N,1);
        n = zeros(N,1);
        vr=exp(j*pi*dir_signal);
        vi=exp(j*pi*dir_gan(xx));
        
        for i = 1:N
             v(i) = exp(sqrt(-1)*(-(N+1)/2+i)*2*pi*k*dir_signal);
             u(i) = exp(sqrt(-1)*(-(N+1)/2+i)*2*pi*k*dir_gan(xx));
             vrr(i) = vr^(i-1);
             vii(i) = vi^(i-1);
         end 
        inr_fact=sqrt(10^(SNR/10));
        snr_fact = sqrt(10^(SNR/10));
        
        signal_power = snr_fact.*exp(j*2*pi*rand);
        gan_power = inr_fact.*exp(j*2*pi*rand);
        Rx=zeros(N);
        Rxr=zeros(N);
               
        for i_temp = 1:K 
        signal_power_w = signal_power.*((randn+sqrt(-1)*randn)/(sqrt(2)));
        gan_power_w = gan_power.*((randn+sqrt(-1)*randn)/(sqrt(2)));
        n=(randn(N,1)+sqrt(-1)*randn(N,1))/(sqrt(2)); 
        x_temp=signal_power_w*v+gan_power_w*u+n;
        xr=signal_power_w*vrr+gan_power_w*vii+n;
        
        Rx=Rx+x_temp*x_temp';
        Rxr=Rxr+xr*xr';
        end
        Rx=Rx/K;
        Rxr=Rxr/K;
        [V,D]=eigs(Rx,2);
        [Vr,Dr]=eigs(Rxr,2);
        
        count=0;
           for mm=-0.25:0.001:0.25
               count=count+1;
               uu(count)=mm;
               vy=zeros(N,1);
               z=exp(j*pi*mm);
               Vz=zeros(N,1);
                       
               for i = 1:N
                   vy(i)=exp(sqrt(-1)*(-(N+1)/2+i)*2*pi*k*mm);
                   Vz(i)=z^(i-1);
               end
               P_bartlett(count)=abs((vy'*Rx*vy)/(vy'*vy));
               Q_capon(count)=abs(vy'*inv(Rx)*vy);
               Q_music(count)=abs(vy'*(eye(10)-V*V')*vy);
               Q_rmusic(count)=abs(Vz'*(eye(10)-Vr*Vr')*Vz);
           end
           
           [tmp1,tmp2]=max(P_bartlett);
           bartlett_tmp=abs(uu(tmp2))^2;
           bartlett(xx)=bartlett(xx)+bartlett_tmp;
           
           
           [tmp1,tmp2]=min(Q_capon);
           capon_tmp=abs(uu(tmp2))^2;
           capon(xx)=capon(xx)+capon_tmp;
           
           [tmp1,tmp2]=min(Q_music);
           music_tmp=abs(uu(tmp2))^2;
           music(xx)=music(xx)+music_tmp;
           
           [tmp1,tmp2]=min(Q_rmusic);
           rmusic_tmp=(abs(uu(tmp2)))^2;
           rmusic(xx)=rmusic(xx)+rmusic_tmp;
       end
   end
   
bartlett=log10(bartlett/m_carlo);
capon=log10(capon/m_carlo);
music=log10(music/m_carlo);
rmusic=log10(rmusic/m_carlo);
   
figure(2);

plot(0:0.05:2,10*bartlett,'b');
hold on;
plot(0:0.05:2,10*capon,'--r');
hold on;
   
plot(0:0.05:2,10*music,'g-.');
hold on;
   
plot(0:0.05:2,10*rmusic,'k:');
hold on;
grid on;
hold off;
legend('bartlett','capon','music','rootmusic');
xlabel('Δu/HPBW');
ylabel('10log(RMSE)');   
title('SNR=20db,K=100,');