% simulateDiversity - script to simulate receiver diversity in Rayleigh fading channels
%
% For a given SNR and diversity order L, simulate the transmission of
% random BPSK symbols with the objective to estimate the bit error
% probability.
%
% There are L narrowband channels between the transmitter and the receiver.
% These channels are modeled as having a Rayleigh distributed gain. The
% channels also introduce white Gaussian noise.
%
% One symbol transmission is simulated at a time; each symbol is sent over
% a new set of channels to guarantee that averaging over the channel
% statistics occurs.
%
% Besides the bit error rate, the outage probability is estimated as well.

%% Simulation parameters
% Set SNR and diversity order
SNRdB = 10;
LL = 4;

% if SNR drops below this level, then an outage occurs
OutageLeveldB = 0;

% stop criterion for simulation
MAX_ERRORS = 400;      % stop after this many errors
MAX_SYMBOLS = 100000;  % or this many symbols

VERBOSE = true;        % print results

% derived parameters
SNR = 10^(SNRdB/10);
OutageLevel = 10^(OutageLeveldB/10);

%% Simulation loop
% initialization
isDone = false;
numErrors = 0;
numSymbols = 0;
numOut = 0;

while ~isDone
   % transmitter
   symbol = 2*(rand > 0.5) - 1; % equally likely +/-1 symbols
   
   % channel
   % complex Gaussian channel gains - variance 1
   hh = 1/sqrt(2)*(randn(1,LL) + 1i*randn(1,LL));
   
   % complex Gaussian noise - variance 1/SNR since a^2=1
   ww = 1/sqrt(2*SNR)*(randn(1,LL) + 1i*randn(1,LL));
   
   % matched filter output
   rxSig = symbol*hh + ww;
   
   % receiver - consists of combiner and decision
   % a really bad combiner, uses just one channel :-(
   rxCombined = rxSig(1)*conj(hh(1)); 
   
   % instantaneous SNR for this receiver
   instSNR = SNR*abs(hh(1)).^2;
   
   % decision
   estSymbol = sign(real(rxCombined));
   
   % update counts
   if estSymbol ~= symbol
       numErrors = numErrors + 1;
   end
   
   if instSNR < OutageLevel
       numOut = numOut + 1;
   end
   
   numSymbols = numSymbols+1;
   
   % update termination criterion
   isDone = numErrors >= MAX_ERRORS || numSymbols >= MAX_SYMBOLS;
   
   if VERBOSE
      if mod(numSymbols, 10000) == 0
          fprintf('Transmitted %d symbols, found %d errors (%6.3g percent)\n', ...
              numSymbols, numErrors, 100*numErrors/numSymbols);
      end
   end
end

%% Report results
if VERBOSE
   fprintf('\n\nTransmitted symbols: %d\n', numSymbols);
   fprintf('Number of symbol Errors: %d\n', numErrors);
   fprintf('Symbol Error Rate: %6.3g\n', numErrors/numSymbols);
   fprintf('95-percent confidence interval: %6.3g\n', 1.96/numSymbols*sqrt(numErrors*(1-numErrors/numSymbols)));
   fprintf('\nNumber of Outages: %d\n', numOut);
   fprintf('Outage Rate: %6.3g\n', numOut/numSymbols);
   fprintf('95-percent confidence interval: %6.3g\n', 1.96/numSymbols*sqrt(numOut*(1-numOut/numSymbols)));
end

