%% Compare BLER Performance of Autoencoders with Uncoded QPSK
% Set simulation parameters.

simParams.EbNoVec = 0:0.5:10;
simParams.MinNumErrors = 200;
simParams.MaxNumFrames = 300;
simParams.NumSymbolsPerFrame = 10000;
simParams.SignalPower = 1;

%
% Simulate a (2,2) autoencoder with energy normalization. The block size is 2 bits. The channel symbol consists of 2 samples, where odd and even numbered samples are modulated on the in-phase and quadrature components of the complex symbol. Each block is carried with one APSK symbol.

load trainedNet_n2_k2_energy txNet rxNet
ae22eBLER = helperAEWAutoencoderBLER(txNet,rxNet,simParams);

%
% Simulate a (4,4) autoencoder with energy normalization. The block size is 4 bits. The channel symbol consists of 4 samples, where odd and even numbered samples are modulated on the in-phase and quadrature components of the complex symbol. Each block is carried with two APSK symbols.

load trainedNet_n4_k4_energy txNet rxNet
ae44eBLER = helperAEWAutoencoderBLER(txNet,rxNet,simParams);

%
% Simulate a (8,8) autoencoder with energy normalization. The block size is 8 bits. The channel symbol consists of 8 samples, where odd and even numbered samples are modulated on the in-phase and quadrature components of the complex symbol. Each block is carried with four APSK symbols.

load trainedNet_n8_k8_energy txNet rxNet
ae88eBLER = helperAEWAutoencoderBLER(txNet,rxNet,simParams);

%
% Simulate a (2,2) uncoded QPSK system. The block size is 2 bits.

n = 2;
k = 2;
qpsk22BLER = qpskBLER(n,k,simParams);

%
% Simulate a (8,8) uncoded QPSK system. The block size is 8 bits.

n = 8;
k = 8;
qpsk88BLER = qpskBLER(n,k,simParams);

%
%% Plot Results

qpsk22BLERTh = 1-(1-berawgn(simParams.EbNoVec,'psk',4,'nondiff')).^2;
semilogy(simParams.EbNoVec,qpsk22BLERTh,':*')
hold on
semilogy(simParams.EbNoVec,qpsk88BLER,'--*')
qpsk88BLERTh = 1-(1-berawgn(simParams.EbNoVec,'psk',4,'nondiff')).^8;
semilogy(simParams.EbNoVec,qpsk88BLERTh,':o')
semilogy(simParams.EbNoVec,ae22eBLER,'-o')
semilogy(simParams.EbNoVec,ae44eBLER,'-d')
semilogy(simParams.EbNoVec,ae88eBLER,'-s')
hold off
ylim([1e-5 1])
grid on
xlabel('E_b/N_o (dB)')
ylabel('BLER')
legend('Uncoded QPSK (2,2)','Uncoded QPSK (8,8)','Theoretical Uncoded QPSK (8,8)',...
  'Autoencoder (2,2)','Autoencoder (4,4)','Autoencoder (8,8)','Location','southwest')
title('BLER performance of R=1 Autoencoders')

% Uncomment next line to save results.

%save uncodedBLERResults ae22eBLER ae44eBLER ae88eBLER qpsk22BLER qpsk88BLER simParams

%
%% Compare BLER Performance of Autoencoders with Hamming Coded QPSK

simParams.EbNoVec = -2:0.5:8;

%
% Simulate a (8,8) uncoded QPSK system. The block size is 8 bits.

n = 4;
k = 4;
qpsk44BLER = qpskBLER(n,k,simParams);

%
% Simulate a (7,4) hard decision decoded Hamming code with QPSK modltion.

n = 7;
k = 4;
hammingHard74BLER = hammingBLER(n,k,'hard',simParams);

%
% Simulate a (7,4) soft decision decoded Hamming code with QPSK modltion.

n = 7;
k = 4;
hammingML74BLER = hammingBLER(n,k,'soft',simParams);

%
% Simulate a (7,4) autoencoder with energy normalization.

load trainedNet_n7_k4_energy txNet rxNet
ae74eBLER = helperAEWAutoencoderBLER(txNet,rxNet,simParams);

%
%% Plot Results

figure
qpsk44BLERTh = 1-(1-berawgn(simParams.EbNoVec,'psk',4,'nondiff')).^4;
semilogy(simParams.EbNoVec,qpsk44BLERTh,':*')
hold on
semilogy(simParams.EbNoVec,qpsk44BLER,':o')
semilogy(simParams.EbNoVec,hammingHard74BLER,'--s')
semilogy(simParams.EbNoVec,ae74eBLER,'-')
semilogy(simParams.EbNoVec,hammingML74BLER,'--d')
hold off
ylim([1e-5 1])
grid on
xlabel('E_b/N_o (dB)')
ylabel('BLER')
legend('Theoretical Uncoded QPSK (4,4)','Uncoded QPSK (4,4)','Hamming (7,4) Hard Decision',...
  'Autoencoder (7,4)','Hamming (7,4) ML','Location','southwest')
title('BLER comparison of (7,4) Autoencoder')

% Uncomment next line to save results.

%save codedBLERResults ae74eBLER hammingHard74BLER hammingML74BLER qpsk44BLER simParams

%
%% Helper Functions

function BLER = qpskBLER(n,k,simParams)
%qpskBLER Block error rate (BLER) simulation for QPSK modulation
%   BLER = qpskBLER(n,k,simParams) estimates the BLER performance of a QPSK
%   system over an AWGN channel using Monte Carlo simulations. k is the
%   number of bits in a block and n is twice the number of QPSK output
%   symbols. simParams is the simulation parameters structure.

EbNoVec = simParams.EbNoVec;

numSymbolsPerFrame = simParams.NumSymbolsPerFrame;
BLER = zeros(size(EbNoVec));
for EbNoIdx = 1:length(EbNoVec)
  EbNo = EbNoVec(EbNoIdx);
  chan = comm.AWGNChannel("BitsPerSymbol",2, ...
    "EbNo", EbNo, "SamplesPerSymbol", 1, "SignalPower", 1);

  numBlockErrors = 0;
  frameCnt = 0;
  while (numBlockErrors < simParams.MinNumErrors) ...
      && (frameCnt < simParams.MaxNumFrames)
    d = randi([0 3],numSymbolsPerFrame,1);
    dBlocks = reshape(d, n/2, []);

    x = pskmod(d, 4);

    y = chan(x);

    dHat = pskdemod(y, 4);
    dHatBlocks = reshape(dHat, n/2, []);

    numBlockErrors = numBlockErrors + sum(any(dBlocks ~= dHatBlocks,1));
    frameCnt = frameCnt + 1;
  end
  numInfoBitsPerFrame = numSymbolsPerFrame*2;
  numBlocks = frameCnt*numInfoBitsPerFrame/k;
  BLER(EbNoIdx) = numBlockErrors / numBlocks;
end
end

function BLER = hammingBLER(n,k,method,simParams)
%hammingBLER Block error rate (BLER) simulation for Hamming codes
%   BLER = hammingBLER(n,k,method,simParams) estimates the BLER performance
%   of a (n,k) Hamming coded QPSK system over an AWGN channel using Monte
%   Carlo simulations. method is the decoding method and must be one of
%   'soft' and 'hard'. simParams is the simulation parameters structure.

EbNoVec = simParams.EbNoVec;
R = k/n;

numSymbolsPerFrame = simParams.NumSymbolsPerFrame;
numBitsPerFrame = numSymbolsPerFrame * k;
BLER = zeros(size(EbNoVec));
parfor EbNoIdx = 1:length(EbNoVec)
  EbNo = EbNoVec(EbNoIdx) + 10*log10(R);
  chan = comm.AWGNChannel("BitsPerSymbol",2, ...
    "EbNo", EbNo, "SamplesPerSymbol", 1, "SignalPower", 1);
  mod = comm.PSKModulator(4,'BitInput',true);
  demod = comm.PSKDemodulator(4,'BitOutput',true);
  if strcmp(method, 'soft')
    demod.DecisionMethod = 'Log-likelihood ratio';
    demod.Variance = 1;
    softDecision = true;
  else
    softDecision = false;
  end

  numBlockErrors = 0;
  frameCnt = 0;
  while (numBlockErrors < simParams.MinNumErrors) ...
      && (frameCnt < simParams.MaxNumFrames) %#ok<PFBNS>
    d = randi([0 1],numBitsPerFrame,1);
    dBlocks = reshape(d, k, []);

    e = encode(d,n,k,'hamming');
    
    x = mod(e);

    y = chan(x);

    eHat = demod(y);

    if softDecision
      dHat = decodeHammingML(eHat,n,k);
    else
      dHat = decode(eHat,n,k,'hamming');
    end
    
    dHatBlocks = reshape(dHat, k, []);
    
    numBlockErrors = numBlockErrors + sum(any(dBlocks ~= dHatBlocks,1));
    frameCnt = frameCnt + 1;
  end
  numBlocks = frameCnt*numBitsPerFrame/k;
  BLER(EbNoIdx) = numBlockErrors / numBlocks;
end
end

function y = decodeHammingML(x,n,k)
%ydecodeHammingML Maximum likelihood decoder for Hamming codes
%   Y = decodeHammingML(X,n,k) decodes received symbols, X, of an (n,k)
%   Hamming code using maximum likelihood decoding and produces data bit
%   estimates, Y.

m = n-k;
h = hammgen(m);
gen = gen2par(h);

messages = de2bi(0:2^k-1, 4);
codeWords = 1-2*rem(messages * gen,2);

numCodeWords = length(x) / n;
xIdx = 1:n;
yIdx = 1:k;
y = zeros(numCodeWords*k,1);
for p=1:numCodeWords
  rcvdCodeWord = x(xIdx,1);
  decisionVar = codeWords * rcvdCodeWord;
  [~,cwIdx] = max(decisionVar);
  y(yIdx,1) = messages(cwIdx,:)';
  xIdx = xIdx + n;
  yIdx = yIdx + k;
end
end