rice1.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;
SNR = 0:5:70; % the system signal-to-noise ratio with dB

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
snr = 10.^(0.1*SNR); % the signal-to-noise with normal scale
snrp=snr; 

% 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);
nn=eye(Td*Nr);       % noise covariance
nn_ray=eye(T*1);     %  Nt=Nr=1 rayleigh fading with perfect csi and without training symbol

capa_awgn= zeros(1,length(SNR));   % AWGN channel capacity

capa_ray = zeros(1,length(SNR));  % rayleigh fading with perfect csi and without training symbol
A_ray=zeros(T,T);   % rayleigh fading with perfect csi and without training symbol

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

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

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

for ii = 1:length(SNR)
waitbar(ii/length(SNR));

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

      for jj=1:MM

      %%% % Nr=Nt=1 rayleigh fading channel with perfect csi and without training symbol
      Hw_ray=sqrt(1/2)*(randn(1,1) + j*randn(1,1));         % channel 

      % %%% Rician Fading channel
      Hw =sqrt(1/2)*(randn(Nr*Nt,1) + j*randn(Nr*Nt,1)); 
      Hw = hh^(1/2)*Hw;
      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(ii))*p*Hw+v;                                        % receive
      g=F*z+f0;                                                   % estimation

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

           for kk =1:NN
           X_ray =sqrt(1/2)* (randn(1,T )+ j* randn(1,T));   % rayleigh fading with perfect csi and without training symbol
           XX_ray=kron(X_ray.',eye(1)); 
           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(ii)*XX*C*hh*XX'+nn;
    
           A30=A30+XX*Hw*Hw'*XX';  % inequality lower bound
           aa300=aa300+1/T*log2(det(eye(Nt*Nr)+snr(ii)*XX'*inv(nn)*XX*C*hh));
           end
       
      A_ray=1/NN*A_ray;   % rayleigh fading with perfect csi and without training symbol
      capa_ray(ii)=capa_ray(ii)+1/T*log2(det(eye(size(A_ray*inv(nn_ray)))+snr(ii)*A_ray*inv(nn_ray)));
     
      A10=1/NN*A10;     % perfect csi
      capa1(ii)=capa1(ii)+1/T*log2(det(eye(Td*Nr)+snr(ii)*A10*inv(nn)));

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

      A30=1/NN*A30;        % inequality lower bound
      capa30(ii)=capa30(ii)+1/T*log2(det(eye(Td*Nr)+snr(ii)*A30*inv(nn)));
      aa300=1/NN*aa300;
      capa300(ii)=aa300+capa300(ii);
      end
capa_awgn(ii)=log2(1+snr(ii));     % AWGN channel capacity
capa_ray(ii)=1/MM*capa_ray(ii);    % rayleigh fading with perfect csi and without training symbol 
capa1(ii)=1/MM*capa1(ii);
capa2(ii)=1/MM*capa2(ii); 
capa3(ii)=1/MM*(capa30(ii)-capa300(ii));   
 
ww= waitbar(ii/length(SNR));
end
close(ww)
% figure
plot(SNR,capa_awgn,'k+-')
hold on
plot(SNR,capa_ray,'ko-')
hold on
plot(SNR,capa1,'ks-')
hold on
plot(SNR,capa2,'kv-')
hold on
plot(SNR,capa3,'k*-')
hold on
grid on
axis([0 70 0 20])
% title('The Ergodic Capacity of MIMO Systems')
xlabel('SNR [dB]')
ylabel('Ergodic Capacity [bit/channel use]') 

legend('AWGN 1x1','Rayleigh 1x1','perfect csi','worst case','inequality')