% Figure 2.20: Capacity with outage 10% for MIMO, SIMO, MISO and SISO channels 
% Elisabeth de Carvalho, April 2012
% Modified December 2015

clear all

% number of channel realization. This number needs to be increased when
% estimating a small outage probability
It = 500000;

fprintf('Total number of channel realizations: %d\n', It);

% SNR range in dB
SNRdBvalues = 0:2:50;

% Selected transmission rate: 10bits/s
R=10;

% initialize variables = number of channel realizations for which the
% selected rate exceeds the maximal achievable rate
N_mimo4x4 = zeros(1,length(SNRdBvalues));
N_mimo2x2 = zeros(1,length(SNRdBvalues));
N_simo1x4 = zeros(1,length(SNRdBvalues));
N_miso4x1 = zeros(1,length(SNRdBvalues));
N_siso1x1 = zeros(1,length(SNRdBvalues));


for kk=1:It
    
    if mod(kk,10000) == 0
        fprintf('Number of channel realizations: %d\n', kk);
    end
    
    % generate channel realization
    Hmimo4x4 = ( randn(4) + j*randn(4) )/sqrt(2);
    Hmimo2x2 = ( randn(2) + j*randn(2) )/sqrt(2);
    Hsimo1x4 = ( randn(4,1) + j*randn(4,1) )/sqrt(2);
    Hmiso4x1 = ( randn(1,4) + j*randn(1,4) )/sqrt(2);
    Hsiso1x1 = ( randn + j*randn )/sqrt(2);

    
    SNRidx = 0;
    for SNRdB = SNRdBvalues
        
        SNR = 10.^(SNRdB./10);  % linear scale
        SNRidx = SNRidx + 1;
                
        % number of channel realizations for which the selected rate exceeds the maximal achievable rate
        % MIMO 4x4
        Cmimo4x4 =  log2(real(det( eye(4) + SNR/4*Hmimo4x4*Hmimo4x4' )));
        N_mimo4x4(SNRidx) = N_mimo4x4(SNRidx) + (1-sign(Cmimo4x4 - R))/2;
        
        % MIMO 2x2
        Cmimo2x2 =  log2(real(det( eye(2) + SNR/2*Hmimo2x2*Hmimo2x2' )));
        N_mimo2x2(SNRidx) = N_mimo2x2(SNRidx) + (1-sign(Cmimo2x2 - R))/2;
        
        % SIMO 1x4
        Csimo1x4 =  log2(real(det( 1 + SNR*norm(Hsimo1x4)^2 )));
        N_simo1x4(SNRidx) = N_simo1x4(SNRidx) + (1-sign(Csimo1x4 - R))/2;
        
        % MISO 4x1
        Cmiso4x1 =  log2(real(det( 1 + SNR/4*norm(Hmiso4x1)^2 )));
        N_miso4x1(SNRidx) = N_miso4x1(SNRidx) + (1-sign(Cmiso4x1 - R))/2;
        
        % SISO 
        Csiso1x1 =  log2(real(det( 1 + SNR*norm(Hsiso1x1)^2 )));
        N_siso1x1(SNRidx) = N_siso1x1(SNRidx) + (1-sign(Csiso1x1 - R))/2;
        
    end
    
end

Pout_mimo4x4 = N_mimo4x4/It
Pout_mimo2x2 = N_mimo2x2/It
Pout_simo1x4 = N_simo1x4/It
Pout_miso4x1 = N_miso4x1/It
Pout_siso1x1 = N_siso1x1/It


%plot
figure(1)
semilogy(SNRdBvalues,Pout_mimo4x4,'r');
hold on
semilogy(SNRdBvalues,Pout_mimo2x2,'g');
semilogy(SNRdBvalues,Pout_simo1x4,'b');
semilogy(SNRdBvalues,Pout_miso4x1,'k');
semilogy(SNRdBvalues,Pout_siso1x1,'m');
xlabel('Average SNR \rho (dB)')
ylabel('Outage Probability')
title('Outage pobability vs SNR for an i.i.d slow Rayleigh fading channel')
grid
hold off

