% Figure 3.21: Error probability for QPSK and different beamforming techniques
% Elisabeth de Carvalho, April 2012


clear all

% number of channel realizations. This number needs to be increased when
% estimating a small probability
It = 50000;

% Number of transmit and receive antennas
M=4;

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


% Transmit power
P=1;

% initialize variables = error probablity for a given channel realization
Psimo = zeros(1,It);
Pmax  = zeros(1,It);
Peig1 = zeros(1,It);
Peig2 = zeros(1,It);
Peig3 = zeros(1,It);
Peig4 = zeros(1,It);


SNRidx = 0;
for SNRdB=SNRdBvalues
    SNRdB
    SNRidx = SNRidx + 1;
    SNR=10^(SNRdB/10);
    
    
    for kk=1:It
        
        
        % generate channel realization
        H     = ( randn(M)  + j*randn(M)   ) / sqrt(2);
        Hsimo = ( randn(M,1)+ j*randn(M,1) ) / sqrt(2);
        
        % SVD of channel matrix H
        % singular values if the diagonal of S are in descending order
        [U S V] = svd(H);
        
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% Capacity for SIMO/MISO
        SNRsimo = Hsimo'*Hsimo*SNR;
        Psimo(kk)= 0.5*erfc(sqrt(0.5*SNRsimo));

        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% Capacity for max Eigenmode transmission
        SNRmax = S(1,1)^2*SNR;
        Pmax(kk) =  0.5*erfc(sqrt(0.5*SNRmax));
        
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% Capacity for Eigenmode tranmission, waterfilling power distribution
        
        %% Waterfilling power allocation for 4 streams
        % Gamma: vector with values of gamma_i, i=1,2,3,4
        Gamma = 1./diag(S).^2/SNR;
        
        % Cut off value assuming all the streams are allocated with power
        gamma0 = (1 + sum(Gamma))/4;
        
        % vector serving to verify if 1/gamma0 is larger than gamma_i
        D = gamma0*ones(4,1) - Gamma ;
        
        % begin by verifying if 1/gamma0 is larger than gamma_i for the
        % waekest stream. If not, no power is allocated.
        if D(4) < 0             % test if power should be allocated to stream 4
            P4      = 0;                                    % no power to stream 4
            gamma0  = (1 + sum(Gamma(1:3)))/3;              % new cut-off value when 3 streams have power
            D       = gamma0*ones(3,1) - Gamma(1:3) ;       % new vector serving to verify if 1/gamma0 is larger than gamma_i
            
            if D(3) < 0         % test if power should be allocated to stream 3
                P3      = 0;                                % no power to stream 3
                gamma0  = (1 + sum(Gamma(1:2)))/2;          % new cut-off value when 2 streams have power
                D       = gamma0*ones(2,1) - Gamma(1:2) ;   % new vector serving to verify if 1/gamma0 is larger than gamma_i
                
                if D(2)<0      % test if power should be allocated to stream 2
                    
                    % power allocated to 1 stream
                    P2 = 0;                         % no power to stream 3
                    P1 = P;                         % all power allocated to stream 1
                    
                else
                    
                    % power allocated to 2 streams
                    P1 = (gamma0 - Gamma(1));       % power allocated to stream 1
                    P2 = (gamma0 - Gamma(2));       % power allocated to stream 2
                    
                end
                
            else
                
                % power allocated to 3 streams
                P1 = (gamma0 - Gamma(1));      % power allocated to stream 1
                P2 = (gamma0 - Gamma(2));      % power allocated to stream 2
                P3 = (gamma0 - Gamma(3));      % power allocated to stream 3
                
            end
            
        else
            
            % power allocated to 4 streams
            P1 = ( gamma0 - Gamma(1) );     % power allocated to stream 1
            P2 = ( gamma0 - Gamma(2) );     % power allocated to stream 2
            P3 = ( gamma0 - Gamma(3) );     % power allocated to stream 3
            P4 = ( gamma0 - Gamma(4) );     % power allocated to stream 4
            
        end
        
        SNR1 = P1*SNR*S(1,1)^2;      % SNR, stream 1
        SNR2 = P2*SNR*S(2,2)^2;      % SNR, stream 2
        SNR3 = P3*SNR*S(3,3)^2;      % SNR, stream 3
        SNR4 = P4*SNR*S(4,4)^2;      % SNR, stream 4
        
        
        
        Peig1(kk) = 0.5*erfc(sqrt(0.5*SNR1));
        Peig2(kk) = 0.5*erfc(sqrt(0.5*SNR2));
        Peig3(kk) = 0.5*erfc(sqrt(0.5*SNR3));
        Peig4(kk) = 0.5*erfc(sqrt(0.5*SNR4));
        
        
        
    end
    
    
    % BER: average over all channel realizations
    PPsimo(SNRidx) = mean(Psimo);
    PPmax (SNRidx) = mean(Pmax);
    PPeig1(SNRidx) = mean(Peig1);
    PPeig2(SNRidx) = mean(Peig2);
    PPeig3(SNRidx) = mean(Peig3);
    PPeig4(SNRidx) = mean(Peig4);
    
end


figure(2)
semilogy(SNRdBvalues, PPsimo,'k-.','linewidth',1)
hold on
semilogy(SNRdBvalues, PPmax,'k--','linewidth',1)
semilogy(SNRdBvalues, PPeig1,'k','linewidth',1)
semilogy(SNRdBvalues, PPeig2,'k','linewidth',1)
semilogy(SNRdBvalues, PPeig3,'k','linewidth',1)
semilogy(SNRdBvalues, PPeig4,'k','linewidth',1)
axis([SNRdBvalues(1) SNRdBvalues(end) 10e-8 1])
legend('SIMO/MISO','Max Eigenmode','Eigenmode')
text(15,0.3,'Eig 4')
text(13.5,0.015,'Eig 3')
text(12,0.00015,'Eig 2')
text(10.5,1.5*10e-7,'Eig 1')
title('Error probability for QPSK transmission')
xlabel('Average SNR \rho (dB)')
ylabel('Error probability')
grid
hold off

