% Figure 3.10: Outage Probability of ZF and MMSE receivers for a slow
% fading channel for rate per stream equal to 2, 4 and 6 bits per
% transmission. ZF and MMSE receivers are not equivalent at high SNR.
% Elisabeth de Carvalho, April 2012

clear all

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

% Number of transmit and receive antennas
M=4;

% SNR range in dB
SNRdBvalues = [0:2:40];

% Selected transmission rate: 2, 4 and 6 bits per transmission
Rvalues = [2 4 6];

% initialize variables = number of channel realizations for which the
% selected rate exceeds the maximal achievable rate
Nzf     = zeros(3,length(SNRdBvalues));
Nmmse   = zeros(3,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;
        
        
        %% Zero-Forcing
        % number of channel realizations for which the selected rate exceeds the maximal achievable rate
        Dzf = inv(H'*H)/SNR*M;
        SNRzf = (1./diag(Dzf));         % output SNR (vector)
        Czf = real( log2(1 + SNRzf(1)) );  % maximal achievable rate
        for RR = 2:2:6
            Nzf(RR/2, SNRidx) = Nzf(RR/2,SNRidx) + .5*(1-sign(Czf-Rvalues(RR/2)));
        end
        
        %% MMSE
        % 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)));
        for RR = 2:2:6
            Nmmse(RR/2, SNRidx) = Nmmse(RR/2,SNRidx) + .5*(1-sign(Cmmse-Rvalues(RR/2)));
        end

    end
    
    
    
end

% Outage Probability
Pzf = Nzf/It
Pmmse = Nmmse/It

% plot
figure(1)
for RR=1:3
    semilogy(SNRdBvalues, Pzf(RR,:),'k-s','linewidth',1)
    hold on
    semilogy(SNRdBvalues, Pmmse(RR,:),'b-+','linewidth',1)
end
title('Outage probability for 1 stream - ZF and MMSE receivers - i.i.d. slow fading channel')
xlabel('Average SNR \rho (dB)')
ylabel('Outage probability')
text(14.5,.11,'Rate = 2')
text(23,.13,'Rate = 4')
text(33,.15,'Rate = 6')
legend('Zero-Forcing Receivers','MMSE Receivers',3)
grid
hold off



