%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BPSK Simulator                            %
% Version 3.4                               %
% By Oscar RODRIGUEZ, Keio University       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Description:                              %
% This simulator creates some random data   %
% bits, simulates their transmission on a   %
% noisy rayleigh fading channel using       %
% Binary Phase Shift Keying and detected    %
% with a Coherent Detector and also with a  %
% maximal ratio combining receiver with 2   %
% receivers. It repeats the simulation for  %
% various signal to noise ratios, and plots %
% the total bit error for the transmission  %
% for each noise levels.                    %
% The simulated results are compared with   %
% the theoretical analysis for BPSK.        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% License terms:                            %
% This code is hereby released to the       %
% public domain. You may use it as you wish %
% for any purpose without any restrictions. %
%                                           %
% You may be able to find the latest        %
% version for this simulator at:            %
% http://www.rapapaing.com/~o-rodrig        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear; 

% CONFIGURABLE VALUES FOR THE SIMULATION:
symbolCount = 10000;                % The total amount of bits to send, the larger this is, the longer the simulation will take.
simStart = -2;                      % Simulation starting point. Eb/No in Db
simEnd = 30;                        % Simulation end point. Eb/No in Db
simGranularity = 0.5;               % Simulation granularity.

% DO NOT EDIT ANYTHING BELOW THIS LINE.
% How many cycles are used to send each symbol.
symbolCycles = 1;
% These are our sample points, equally distributed in decibels...
simDataDb = (simStart:simGranularity:simEnd)';

% These are our sample points, for EbN0
simEbN0 = 10.^(simDataDb ./ 10);

% These are our sample points, for noise energies
simData = 1 ./ simEbN0;

% And this is where we'll store our results
simErrorAWGN = 0 * simData;
simErrorFading = 0 * simData;
simErrorDiversity = 0 * simData;

% Calculate the theoretical data...
theoryAWGN = 0.5 * erfc(sqrt(simEbN0));
theoryFading = 0.5 * (1 - sqrt((simEbN0) ./ (1 + simEbN0)));
theoryDiversity = ((1 ./ (4 * (simEbN0 ./ 2))) .^ 2) .* 3;

% Set up some common values and arrays...
startTime = 1;
totalTime = symbolCount * symbolCycles;
t = (startTime:1:totalTime)';

% Create some random noise...
noise = (sqrt(2) / 2) * (randn(size(t,1),1) + i * randn(size(t,1),1));
noise2 = (sqrt(2) / 2) * (randn(size(t,1),1) + i * randn(size(t,1),1));

% Create our Rayleigh channels...
fading = (sqrt(2) / 2) * (randn(size(t,1),1) + i * randn(size(t,1),1));
fading2 = (sqrt(2) / 2) * (randn(size(t,1),1) + i * randn(size(t,1),1));

% Start our simulations!
for r=(1:size(simData))
    % Setup our values...
    noiseEnergy = simData(r); 
    noiseSigma = sqrt(noiseEnergy);
    mistakenSymbolsAWGN = 0;
    mistakenSymbolsFading = 0;
    mistakenSymbolsDiversity = 0;

    % Create some random bits...
    symbols = sign(rand(symbolCount, 1) - 0.5);
    
    % Setup our data array...
    data = 0 * t;
    for i=(0:symbolCount-1)
        for j=(1:symbolCycles)
            data(i*symbolCycles+j) = symbols(i+1);
        end
    end
    
    % Encode using PSK...
    encodedSignal = data;
    
    % On the receiving end, not only the signal, but also some noise was
    % detected...
    decodedSignalAWGN = encodedSignal + (noiseSigma * noise);
    decodedSignalFading = (fading .* encodedSignal) + (noiseSigma * noise);
    decodedSignalFading2 = (fading2 .* encodedSignal) + (noiseSigma * noise2);
    decodedSignalDiversity = (conj(fading) .* decodedSignalFading + conj(fading2) .* decodedSignalFading2);
    
    % Extract the detected bits...
    decodedSymbolsAWGN = sign(real(decodedSignalAWGN));
    decodedSymbolsFading = sign(real(decodedSignalFading ./ fading));
    decodedSymbolsDiversity = sign(real(decodedSignalDiversity));
    for i=(1:size(decodedSymbolsAWGN))
        % Count the amount of mistaken symbols...
        if (decodedSymbolsAWGN(i) ~= symbols(i))
            mistakenSymbolsAWGN = mistakenSymbolsAWGN + 1;
        end
        if (decodedSymbolsFading(i) ~= symbols(i))
            mistakenSymbolsFading = mistakenSymbolsFading + 1;
        end
        if (decodedSymbolsDiversity(i) ~= symbols(i))
            mistakenSymbolsDiversity = mistakenSymbolsDiversity + 1;
        end
    end

    % Now store the results
    simErrorAWGN(r) = mistakenSymbolsAWGN / symbolCount;
    simErrorFading(r) = mistakenSymbolsFading / symbolCount;
    simErrorDiversity(r) = mistakenSymbolsDiversity / symbolCount;
end

figure('Name', sprintf('Simulation results for %g bits', symbolCount));

% Plot the results...
semilogy(simDataDb,simErrorAWGN,'rx',simDataDb,theoryAWGN,'r',simDataDb,simErrorFading,'bx',simDataDb,theoryFading,'b',simDataDb,simErrorDiversity,'gx',simDataDb,theoryDiversity,'g');
set(gca,'XTick',simStart:2:simEnd);
set(gca,'YMinorTick','on');
set(legend('BPSK simulation results for AWGN','BPSK theoretical error probability for AWGN','BPSK simulation results for fading channel','BPSK theoretical error probability for fading channel','BPSK simulation results for fading channel with MRC','BPSK theoretical error probability approximation for fading channel with MRC','Location','SouthOutside'),'Interpreter','none')
xlim([simStart,simEnd]);
ylim([1/symbolCount,1]);
title('Bit error probability for BPSK');
xlabel('Eb/No (dB)')
ylabel('Bit error probability')