% Figure 3.11: Outage Probability of nonordered SIC for each layer of a 4x4
% MIMO system. Transmission rate is 4 bits per transmission. 
% Elisabeth de Carvalho, April 2012

clear all

% number of channel realizations. This number needs to be increased when
% estimating a small probability
It = 50000;
fprintf('Total number of channel realizations: %d\n', It);

% Number of transmit and receive antennas
M=4;

% SNR range in dB
SNRdBvalues = [10:2:30];

% Selected transmission rate: 2, 4 and 6 bits per transmission
RR = 4;

% initialize variables = number of channel realizations for which the
% selected rate exceeds the maximal achievable rate
Nsic1   = zeros(1,length(SNRdBvalues));
Nsic2   = zeros(1,length(SNRdBvalues));
Nsic3   = zeros(1,length(SNRdBvalues));
Nsic4   = 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
    H = (randn(M)+j*randn(M))/sqrt(2);
    
    SNRidx = 0;
    for SNRdB = SNRdBvalues
        
        SNR = 10.^(SNRdB./10);  % linear scale
        SNRidx = SNRidx + 1;
        
        %% Layer 1
        % number of channel realizations for which the selected rate exceeds the maximal achievable rate
        Dmmse = inv(H'*H + eye(M)/SNR*M);
        SNRmmse = (1./diag(Dmmse))*SNR/M -ones(M,1);
        Cmmse = (log2( 1 + SNRmmse(1)));
        Nsic1(SNRidx) = Nsic1(SNRidx) + .5*(1-sign(Cmmse-RR));
        
        %% Layer 2
        % number of channel realizations for which the selected rate exceeds the maximal achievable rate
        H2 = H(:,2:M);
        Dmmse = inv(H2'*H2 + eye(M-1)/SNR*M);
        SNRmmse = (1./diag(Dmmse))*SNR/M -ones(M-1,1);
        Cmmse = (log2( 1 + SNRmmse(1)));
        Nsic2(SNRidx) = Nsic2(SNRidx) + .5*(1-sign(Cmmse-RR));
        
        %% Layer 3
        % number of channel realizations for which the selected rate exceeds the maximal achievable rate
        H3 = H(:,3:M);
        Dmmse = inv(H3'*H3 + eye(M-2)/SNR*M);
        SNRmmse = (1./diag(Dmmse))*SNR/M -ones(M-2,1);
        Cmmse = (log2( 1 + SNRmmse(1)));
        Nsic3(SNRidx) = Nsic3(SNRidx) + .5*(1-sign(Cmmse-RR));
        
        %% Layer 4
        % number of channel realizations for which the selected rate exceeds the maximal achievable rate
        H4 = H(:,4:M);
        Dmmse = inv(H4'*H4 + eye(M-3)/SNR*M);
        SNRmmse = (1./diag(Dmmse))*SNR/M -ones(M-3,1);
        Cmmse = (log2( 1 + SNRmmse(1)));
        Nsic4(SNRidx) = Nsic4(SNRidx) + .5*(1-sign(Cmmse-RR));
        
    end
    
end

% Outage Probability
Psic1 = Nsic1/It
Psic2 = Nsic2/It
Psic3 = Nsic3/It
Psic4 = Nsic4/It

% plot
figure(1)
semilogy(SNRdBvalues, Psic1,'k')
hold on
semilogy(SNRdBvalues, Psic2,'b')
semilogy(SNRdBvalues, Psic3,'m')
semilogy(SNRdBvalues, Psic4,'r')
title('Outage probability of nonordered SIC for rate 4 bits/transmission')
xlabel('Average SNR \rho (dB)')
ylabel('Outage probability')
legend('Layer 1','Layer 2','Layer 3','Layer 4',3)
grid
hold off



