cout_bijiao.asv

分别给出单天线和多天线下的MIMO－仿真

clc
clear all
Nt = 2; % the number of transmitted antennas
Nr = 2; % the number of transmitted antennas

T=4; % the number of symbol periods
Tp=Nt;
Td=T-Tp;
%Td(1,length(T))=T-Tp*ones(1,length(T));

SNR = 30; % the system signal-to-noise ratio with dB
snr = 10.^(0.1*SNR); % the signal-to-noise with normal scale
snrp=snr; 
cout=0:1:25;    %outage capacity
MM=100;    % 1000 channel
NN =100;% using 10000 Monte-Carlo runs
% Nr=Nt=2 rayleigh fading channel with perfect csi and without training symbol
pout_ray=zeros(1,length(cout));   % outage probability
nout_ray=zeros(1,length(cout)); 
%%%%% LMMSE estimator
%%%perfect csi
pouta_lmmse=zeros(1,length(cout));   % outage probability
nouta_lmmse=zeros(1,length(cout));                      % the number of capa(jj)<cout
capa1_lmmse = zeros(1,MM);  % perfect csi
A10_lmmse=zeros(Td*Nr,Td*Nr);
%%%%worst case noise lower bound
poutb_lmmse=zeros(1,length(cout));   % outage probability
noutb_lmmse=zeros(1,length(cout));                      % the number of capa(jj)<cout
capa2_lmmse= zeros(1,MM);    % worst case noise lower bound
%%%%inequality lower bound
poutc_lmmse=zeros(1,length(cout));   % outage probability
noutc_lmmse=zeros(1,length(cout));                      % the number of capa(jj)<cout
capa3_lmmse = zeros(1,MM);    % inequality lower bound
capa30_lmmse=zeros(1,MM);
capa300_lmmse=zeros(1,MM);

%%%%% ML estimator
%%%perfect csi
pouta_ml=zeros(1,length(cout));   % outage probability
nouta_ml=zeros(1,length(cout));                      % the number of capa(jj)<cout
capa1_ml = zeros(1,MM);  % perfect csi
%%%%worst case noise lower bound
poutb_ml=zeros(1,length(cout));   % outage probability
noutb_ml=zeros(1,length(cout));                      % the number of capa(jj)<cout
capa2_ml= zeros(1,MM);    % worst case noise lower bound
%%%%inequality lower bound
poutc_ml=zeros(1,length(cout));   % outage probability
noutc_ml=zeros(1,length(cout));                      % the number of capa(jj)<cout
capa3_ml = zeros(1,MM);    % inequality lower bound
capa30_ml=zeros(1,MM);
capa300_ml=zeros(1,MM);

tt = 0.5; % the correlation coefficients of transmitterwith exponential correlated , equal to zero when the channel is transmitted uncorrelated
% generating the correlated matrix at Tx
K =0.2; % Rician fading factor, and it will be qqual to zero when the channel is Rayleigh fading
for i =1:Nt
for j=1:Nt
Rt(i,j) = tt^(abs(i-j));
end
end

rr = 0.5; % the correlation coefficients of receiver with exponential correlated , equal to zero when the channel is received uncorrelated
% generating the correlated matrix at Rx
for i =1:Nr
for j=1:Nr
Rr(i,j) = rr^(abs(i-j));
end
end

% hm=zeros(Nt*Nr,1);   % channel mean 
hm=sqrt(K/(1+K))*ones(Nr*Nt,1);
hh=sqrt(1/(1+K))*kron(Rt.',Rr);       % channel covariance
vv=eye(Tp*Nr);

capa_ray = zeros(1,length(MM));  % Nr=Nt=2 rayleigh fading channel with perfect csi and without training symbol 

nn=eye(Td*Nr);       % noise covariance

nn_ray=eye(T*Nr);     % % Nr=Nt=2 rayleigh fading channel with perfect csi and without training symbol
A_ray=zeros(T*Nr,T*Nr); 



A20_lmmse=zeros(Td*Nr,Td*Nr);
AA20_lmmse=zeros(Td*Nr,Td*Nr);
A20_ml=zeros(Td*Nr,Td*Nr);
AA20_ml=zeros(Td*Nr,Td*Nr);

A30=zeros(Td*Nr,Td*Nr);
AA30=zeros(Td*Nr,Td*Nr);
aa300=0;

p=eye(Tp*Nr);       % training symbol
pp=inv(p'*p)*p';    % wei ni 
nn=eye(Td*Nr);       % noise covariance

%%%%% LMMSE estimator
F_lmmse=sqrt(snrp)*hh*p'*inv(snrp*p*hh*p'+vv);
f0_lmmse=(eye(size(F_lmmse*p))-sqrt(snrp)*F_lmmse*p)*hm;
% %%%ML estimator
F_ml=1/sqrt(snrp)*pp;
f0_ml=zeros(Nt*Nr,1);

for jj=1:MM
waitbar(jj/MM);
      %%% % Nr=Nt=2 rayleigh fading channel with perfect csi and without training symbol
      Hw_ray=sqrt(1/2)*(randn(Nr*Nt,1) + j*randn(Nr*Nt,1));          % channel 
      % %%% Rician Fading channel correlated 
      % Hw =sqrt(1/2)*(randn(Nr*Nt,1) + j*randn(Nr*Nt,1)); 
      Hw = hh^(1/2)*Hw_ray;
      Hw= hm + sqrt(1/(1+K))*Hw;
      v=sqrt(1/2)*(randn(Nr*Tp,1) + j*randn(Nr*Tp,1));          % training noise
      z=sqrt(snrp)*p*Hw+v;                                        % received symbol
      % %%%%%%LMMSE estimation 
     g_lmmse=F_lmmse*z+f0_lmmse;         
     B_lmmse=sqrt(snrp)*hh*p'*F_lmmse'*inv(F_lmmse*(snrp*p*hh*p'+vv)*F_lmmse');
     C_lmmse=eye(size(B_lmmse*F_lmmse*p))-sqrt(snrp)*B_lmmse*F_lmmse*p;
     d_lmmse=B_lmmse*g_lmmse+C_lmmse*hm-B_lmmse*f0_lmmse;
      % %%%%%ML estimation
     g_ml=F_ml*z+f0_ml;     
     B_ml=sqrt(snrp)*hh*p'*F_ml'*inv(F_ml*(snrp*p*hh*p'+vv)*F_ml');
     C_ml=eye(size(B_ml*F_ml*p))-sqrt(snrp)*B_ml*F_ml*p;
     d_ml=B_ml*g_ml+C_ml*hm-B_ml*f0_ml;
     
          for kk =1:NN
             X_ray =sqrt(1/(2*Nt))* (randn(Nt,T) + j* randn(Nt,T));   % rayleigh fading with perfect csi and without training symbol
             XX_ray=kron(X_ray.',eye(Nr)); 
             A_ray=A_ray+XX_ray*Hw_ray*Hw_ray'*XX_ray';   % rayleigh fading with perfect csi and without training symbol
             
             X =sqrt(1/(2*Nt))* (randn(Nt,Td) + j* randn(Nt,Td)); 
             XX=kron(X.',eye(Nr)); 
    
             % %%%%%%LMMSE estimation 
             A10_lmmse=A10_lmmse+XX*Hw*Hw'*XX';   % perfect csi
             A20_lmmse=A20_lmmse+XX*d_lmmse*d_lmmse'*XX';    % worst case noise lower bound
             AA20_lmmse=AA20_lmmse+snrp*XX*C_lmmse*hh*XX'+nn;
             A30_lmmse=A30_lmmse+XX*Hw*Hw'*XX';  % inequality lower bound
             aa300_lmmse=aa300_lmmse+1/T*log2(det(eye(Nt*Nr)+snr*XX'*inv(nn)*XX*C_lmmse*hh));
             % %%%%% ML estimation
             A10_ml=A10_ml+XX*Hw*Hw'*XX';   % perfect csi
             A20_ml=A20_ml+XX*d_ml*d_ml'*XX';    % worst case noise lower bound
             AA20_ml=AA20_ml+snrp*XX*C_ml*hh*XX'+nn;
             A30_ml=A30_ml+XX*Hw*Hw'*XX';  % inequality lower bound
             aa300_ml=aa300_ml+1/T*log2(det(eye(Nt*Nr)+snr*XX'*inv(nn)*XX*C_ml*hh));
         end
    
    A_ray=1/NN*A_ray;   %% Nr=Nt=2 rayleigh fading channel with perfect csi and without training symbol
    capa_ray(jj)=1/T*log2(det(eye(T*Nr)+snr*A_ray*inv(nn_ray))); 
       for  ii=1:length(cout)
           if capa_ray(jj)<=cout(ii)
           nout_ray(ii)=nout_ray(ii)+1
           end  
       end
       
    % %%%%%%LMMSE estimation
    %%%% % perfect csi
    A10_lmmse=1/NN*A10_lmmse;    
    capa1_lmmse(jj)=1/T*log2(det(eye(Td*Nr)+snr*A10_lmmse*inv(nn)));
     for  ii=1:length(cout)
         if capa1_lmmse(jj)<=cout(ii)
           nouta_lmmse(ii)=nouta_lmmse(ii)+1;
         end
     end
   % %%%%%worst case noise lower bound
   A20_lmmse=1/NN*A20_lmmse;      
   AA20_lmmse=1/NN*AA20_lmmse;
   capa2_lmmse(jj)=1/T*log2(det(eye(Td*Nr)+snr*A20_lmmse*inv(AA20_lmmse)));
     for  ii=1:length(cout)    
         if capa2_lmmse(jj)<=cout(ii)
           noutb_lmmse(ii)=noutb(ii)+1;
         end
     end
   % %%%%inequality lower bound  
   A30_lmmse=1/NN*A30_lmmse;       
   capa30_lmmse(jj)=1/T*log2(det(eye(Td*Nr)+snr*A30_lmmse*inv(nn)));
   aa300_lmmse=1/NN*aa300_lmmse;
   capa300_lmmse(jj)=aa300_lmmse;
   capa3_lmmse(jj)=capa30_lmmse(jj)-capa300_lmmse(jj);
    for  ii=1:length(cout)  
        if capa3_lmmse(jj)<=cout(ii)
           noutc_lmmse(ii)=noutc(ii)+1;
        end
    end
    
    % %%%%%%ML estimation
        %%%% % perfect csi
    A10_ml=1/NN*A10_ml;    
    capa1_ml(jj)=1/T*log2(det(eye(Td*Nr)+snr*A10_ml*inv(nn)));
     for  ii=1:length(cout)
         if capa1_ml(jj)<=cout(ii)
           nouta_ml(ii)=nouta_ml(ii)+1;
         end
     end
   % %%%%%worst case noise lower bound
   A20_ml=1/NN*A20_ml;      
   AA20_ml=1/NN*AA20_ml;
   capa2_ml(jj)=1/T*log2(det(eye(Td*Nr)+snr*A20_ml*inv(AA20_ml)));
     for  ii=1:length(cout)    
         if capa2_ml(jj)<=cout(ii)
           noutb_ml(ii)=noutb(ii)+1;
         end
     end
   % %%%%inequality lower bound  
   A30_ml=1/NN*A30_ml;       
   capa30_ml(jj)=1/T*log2(det(eye(Td*Nr)+snr*A30_ml*inv(nn)));
   aa300_ml=1/NN*aa300_ml;
   capa300_ml(jj)=aa300_ml;
   capa3_ml(jj)=capa30_ml(jj)-capa300_ml(jj);
    for  ii=1:length(cout)  
        if capa3_ml(jj)<=cout(ii)
           noutc_ml(ii)=noutc(ii)+1;
        end
    end
end
pout_ray=1/MM*nout_ray;
% %%%%%%LMMSE estimation
pouta_lmmse=1/MM*nouta_lmmse;
poutb_lmmse=1/MM*noutb_lmmse;
poutc_lmmse=1/MM*noutc_lmmse;
% %%%%%%ML estimation
pouta_ml=1/MM*nouta_ml;
poutb_ml=1/MM*noutb_ml;
poutc_ml=1/MM*noutc_ml;
ww=waitbar(jj/MM);
close(ww)

% figure
plot(cout,pout_ray,'ko-')
hold on
% %%%%%%LMMSE estimation
plot(cout,pouta_lmmse,'ks-')
hold on
plot(cout,poutb_lmmse,'kv-')
hold on
plot(cout,poutc_lmmse,'k*-')
hold on
% %%%%%%ML estimation
plot(cout,pouta_ml,'bs-')
hold on
plot(cout,poutb_ml,'bv-')
hold on
plot(cout,poutc_ml,'b*-')
hold on
grid on
% title('outage probability of MIMO system's outage capacity')
xlabel('Cout [bit/channel use]')
ylabel('outage probability ') 

% legend('Rayleigh 2x2','perfect csi','worst case','inequality')