coutlmmse.m

分别给出单天线和多天线下的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
% 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)); 
%%%perfect csi
pouta=zeros(1,length(cout));   % outage probability
nouta=zeros(1,length(cout));                      % the number of capa(jj)<cout
%%%%worst case noise lower bound
poutb=zeros(1,length(cout));   % outage probability
noutb=zeros(1,length(cout));                      % the number of capa(jj)<cout
%%%%inequality lower bound
poutc=zeros(1,length(cout));   % outage probability
noutc=zeros(1,length(cout));                      % the number of capa(jj)<cout

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

MM=100;    % 1000 channel
NN =100;% using 10000 Monte-Carlo runs

% 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 
capa1 = zeros(1,MM);  % perfect csi
capa2 = zeros(1,MM);    % worst case noise lower bound
capa3 = zeros(1,MM);    % inequality lower bound
capa30=zeros(1,MM);
capa300=zeros(1,MM);

for ii = 1:length(cout)
waitbar(ii/length(cout));
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);  

% capa1 = zeros(1,length(T));  % perfect csi
A10=zeros(Td*Nr,Td*Nr);

%capa2 = zeros(1,length(T));    % worst case noise lower bound
A20=zeros(Td*Nr,Td*Nr);
AA20=zeros(Td*Nr,Td*Nr);


% capa3 = zeros(1,length(T));    % inequality lower bound
% capa30=zeros(1,length(T));
% capa300=zeros(1,length(T));
A30=zeros(Td*Nr,Td*Nr);
AA30=zeros(Td*Nr,Td*Nr);
aa300=0;

%%%%% LS estimator
% p=eye(Nt*Nr);       % training symbol
% pp=inv(p'*p)*p';    % wei ni 
% F=1/sqrt(snrp)*pp;
% f0=zeros(Nt*Nr,1);     
%%%%% LMMSE estimator
p=eye(Nt*Nr);       % training symbol
pp=inv(p'*p)*p';    % wei ni 
F=sqrt(snrp)*hh*p'*inv(snrp*p*hh*p'+vv');
f0=(eye(size(F*pp))-sqrt(snrp)*F*pp)*hm;
nn=eye(Td*Nr);       % noise covariance

      for jj=1: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;                                        % receive
     g=F*z+f0;                                                   % estimation

     B=eye(Nt*Nr);
     C=inv(eye(Nt*Nr)+snrp*hh*p'*inv(vv)*p);
     d=B*g+C*hm-B*f0;

             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)); 
    
             A10=A10+XX*Hw*Hw'*XX';   % perfect csi
    
             A20=A20+XX*d*d'*XX';    % worst case noise lower bound
             AA20=AA20+snrp*XX*C*hh*XX'+nn;
    
             A30=A30+XX*Hw*Hw'*XX';  % inequality lower bound
             aa300=aa300+1/T*log2(det(eye(Nt*Nr)+snr*XX'*inv(nn)*XX*C*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))); 
       if capa_ray(jj)<=cout(ii)
           nout_ray(ii)=nout_ray(ii)+1
       end
     
     A10=1/NN*A10;     % perfect csi
     capa1(jj)=1/T*log2(det(eye(Td*Nr)+snr*A10*inv(nn)));
       if capa1(jj)<=cout(ii)
           nouta(ii)=nouta(ii)+1;
       end

     A20=1/NN*A20;      % worst case noise lower bound
     AA20=1/NN*AA20;
     capa2(jj)=1/T*log2(det(eye(Td*Nr)+snr*A20*inv(AA20)));
       if capa2(jj)<=cout(ii)
           noutb(ii)=noutb(ii)+1;
       end

     A30=1/NN*A30;        % inequality lower bound
     capa30(jj)=1/T*log2(det(eye(Td*Nr)+snr*A30*inv(nn)));
     aa300=1/NN*aa300;
     capa300(jj)=aa300;
     capa3(jj)=capa30(jj)-capa300(jj);
       if capa3(jj)<=cout(ii)
           noutc(ii)=noutc(ii)+1;
       end
    end
% capa1(ii)=1/MM*capa1(ii);
% capa2(ii)=1/MM*capa2(ii); 
% capa3(ii)=1/MM*(capa30(ii)-capa300(ii));  
pout_ray(ii)=1/MM*nout_ray(ii);
pouta(ii)=1/MM*nouta(ii);
poutb(ii)=1/MM*noutb(ii);
poutc(ii)=1/MM*noutc(ii);
 
ww= waitbar(ii/length(cout));
end
close(ww)
% figure
plot(cout,pout_ray,'ko-')
hold on
plot(cout,pouta,'ks-')
hold on
plot(cout,poutb,'kv-')
hold on
plot(cout,poutc,'k*-')
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')