% Figure 3.28: Outage probability of OSTBC and QOSTBC for a selected rate
% R=2 bits per transmission and a 4x1 MISO i.i.d. Rayleigh slow fading channel.
% Elisabeth de Carvalho, April 2012


clear all

% number of channel realizations. This number needs to be increased when
% estimating small probabilities
It = 100000;

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


% transmission rate: 4 bits per tranmission
RR = 2;


SNRidx=0;
for SNRdB=SNRdBvalues
    SNRdB
    SNRidx = SNRidx + 1;    % SNR index
    SNR = 10^(SNRdB/10);    % SNR in linear scale
    
    % initialize variables: number of channel realizations where RR exceeds
    % the maximal achievable rate 
    
    % OSTBC: 1 variable per stream
    Nostbc1=0;
    Nostbc2=0;
    Nostbc3=0;
    Nostbc4=0;
    
    % OSTBC: 1 variable per stream
    Nqostbc1=0;
    Nqostbc2=0;
    Nqostbc3=0;
    Nqostbc4=0;
    
    Nala=0;     % Alamouti
    Nmiso=0;    % MISO
    Nsiso=0;    % SISO

    
    % iterations on channel realizations
    for kk=1:It
        
        % generate channels
        H  = (randn(1,4) + j * randn(1,4)) / sqrt(2);        % 1x4 MISO channel
        h1 = H(1,1); h2 = H(1,2); h3 = H(1,3); h4= H(1,4);   % subchannels
        
        %%%% Orthogonal STBC
        % equivalent channel
        HH_o = [h1 h2 h3 h4;
            h2 -h1 h4 -h3;
            h3 -h4 -h1 h2;
            h4 h3 -h2 -h1;
            h1' h2' h3' h4';
            h2' -h1' h4' -h3';
            h3' -h4' -h1' h2';
            h4' h3' -h2' -h1'];
        % compute outage events for each stream
        R = inv(HH_o'*HH_o);
        SNR1 = SNR/4/R(1,1);                        % post-processing SNR stream 1
        C1 = real(log2(1 + SNR1 ))/2;               % Capacity for stream 1
        Nostbc1 = Nostbc1 + .5*(1-sign(C1-RR));     % number of outage events
        SNR2 = SNR/4/R(2,2);                        % post-processing SNR stream 2
        C2 = real(log2(1 + SNR2 ))/2;               % Capacity for stream 2
        Nostbc2 = Nostbc2 + .5*(1-sign(C2-RR));     % number of outage events
        SNR3 = SNR/4/R(3,3);                        % post-processing SNR stream 3
        C3 = real(log2(1 + SNR3 ))/2;               % Capacity for stream 3
        Nostbc3 = Nostbc3 + .5*(1-sign(C3-RR));     % number of outage events
        SNR4 = SNR/4/R(4,4);                        % post-processing SNR stream 4
        C4 = real(log2(1 + SNR4 ))/2;               % Capacity for stream 4
        Nostbc4 = Nostbc4 + .5*(1-sign(C4-RR));     % number of outage events
        
        %%%%  Quasi-Orthogonal STBC
        % equivalent channel
        HH_qo = [h1  -h2  -h3   h4;
            h2'  h1' -h4' -h3';
            h3' -h4'  h1' -h2';
            h4   h3   h2   h1];
        % compute outage events for each stream
        R = inv(HH_qo'*HH_qo + eye(4)/(SNR/4) );
        SNR1 = SNR/4/R(1,1) - 1;                    % post-processing SNR stream 1
        C1 = real(log2(1 + SNR1 ));                 % Capacity for stream 1
        Nqostbc1 = Nqostbc1 + .5*(1-sign(C1-RR));   % number of outage events
        SNR2 = SNR/4/R(2,2) -1;                     % post-processing SNR stream 2
        C2 = real(log2(1 + SNR2 ));                 % Capacity for stream 2
        Nqostbc2 = Nqostbc2 + .5*(1-sign(C2-RR));   % number of outage events
        SNR3 = SNR/4/R(3,3) - 1;                    % post-processing SNR stream 3
        C3 = real(log2(1 + SNR3 ));                 % Capacity for stream 3
        Nqostbc3 = Nqostbc3 + .5*(1-sign(C3-RR));   % number of outage events
        SNR4 = SNR/4/R(4,4) - 1;                    % post-processing SNR stream 4
        C4 = real(log2(1 + SNR4 ));                 % Capacity for stream 4
        Nqostbc4 = Nqostbc4 + .5*(1-sign(C4-RR));   % number of outage events
        
        
        %%%%  Alamouti STBC, transmission from antenna 1 and antenna 2
        SNRala = SNR/2*(abs(h1)^2 + abs(h2)^2);
        Cala = (log2(1 + SNRala ))   ;
        Nala = Nala + .5*(1-sign(Cala-RR));
        
        
        %%%%  MISO
        Cmiso = real(log2( 1 + SNR/4*H*H' ));
        Nmiso = Nmiso + .5*(1-sign(Cmiso-RR));
        
        
        %%%%  SISO
        SNRsiso = SNR*abs(h1)^2;
        Csiso = real(log2(1 + SNRsiso ))   ;
        Nsiso = Nsiso + .5*(1-sign(Csiso-RR));
        
    end
    
    
    % Outage probabilities
    Postbc(SNRidx)  = (Nostbc1  + Nostbc2  + Nostbc3  + Nostbc4)/It/4
    Pqostbc(SNRidx) = (Nqostbc1 + Nqostbc2 + Nqostbc3 + Nqostbc4)/It/4
    Pala(SNRidx)    = Nala/It
    Pmiso(SNRidx)   = Nmiso/It
    Psiso(SNRidx)   = Nsiso/It
    
    
end


figure(1)
semilogy(SNRdBvalues, Postbc,'-*','linewidth',1)
hold on
semilogy(SNRdBvalues, Pqostbc,'-o','linewidth',1)
semilogy(SNRdBvalues, Pala,'-+','linewidth',1)
semilogy(SNRdBvalues, Psiso,'--','linewidth',1)
semilogy(SNRdBvalues, Pmiso,'-','linewidth',1)
axis([0 35 10e-7 1])
grid
title('Outage probability for 1 stream - 4x1 MISO i.i.d. Rayleigh slow fading channel')
xlabel('Average SNR \rho (dB)')
ylabel('Outage probability')
legend('O STBC','QO STBC, MMSE','Alamouti STBC','SISO','MISO',3)
hold off
