SDD, HDD, AutoEncoder Comparison

clear all
close all

n = 7; % Codeword length (No. bits per Codeword)
k = 4; % Message length  (No. bits per Message)
M = 2; % BPSK
dmin = 3; % minimum distance
EbN0dB = 0:0.5:8; % in dB

% No. Transmitted Bits
noBits = k*log2(M)*1e4; % (must be a multiple of k and log2(M))

% Uncoded Bitstream (Serial)
uncoded_bits = randi([0 1],[noBits,1]);

% Generator Matrix (G) and Parity Check Matrix (H) of (7,4) Hamming Code
[H,G] = hammgen(log2(n+1)); % 3 = log2(n+1);

% Channel Encoding
channel_coding_bits = encode(uncoded_bits, n, k, 'hamming/binary');

% List of codewords for (n, k) Hamming code
U_ref_matrix = encode(dec2bin(0:(2^k-1),k)-'0', n, k, 'hamming/binary');
% UBPSK = real(pskmod(U_ref_matrix,M,InputType="bit"));

% MPSK Modulation with Gray Mapping, i.e., InputType="bit"
% modulate only source coding (no channel coding) output
modulated_uncoded_MPSK= pskmod(uncoded_bits, M, InputType="bit");

% modulate channel coding output
modulated_coded_MPSK  = pskmod(channel_coding_bits, M, InputType="bit");

% figure;
% scatterplot(modulated_coded_MPSK);

% Unit Energy AWGN
AWGN = @(x) sqrt(1/2)*(randn(size(x))+1i*randn(size(x)));

for iSNR = 1:length(EbN0dB)
    EbN0 = db2pow(EbN0dB(iSNR));
    EcN0 = EbN0*(k/n);

    %% Uncoded MPSK
    % Received Signal given a specific EbN0
    received_uncoded_MPSK = modulated_uncoded_MPSK + 1/sqrt(EbN0*log2(M))*AWGN(modulated_uncoded_MPSK);

    %     figure;
    %     scatterplot(received_uncoded_MPSK);
    
    % MPSK Demodulation (Maximum Likelihood Detection)
    demodulated_uncoded_MPSK = pskdemod(received_uncoded_MPSK, M, OutputType='bit');

    % Bit Error Rate
    BER_uncoded_sim(iSNR)= mean(xor(uncoded_bits, demodulated_uncoded_MPSK));

    %% Coded MPSK
    % Received Signal given a specific EbN0
    received_coded_MPSK = modulated_coded_MPSK + 1/sqrt(EcN0*log2(M))*AWGN(modulated_coded_MPSK);

    % MPSK Demodulation (Maximum Likelihood Detection)
    demodulated_coded_MPSK  = pskdemod(received_coded_MPSK, M, OutputType='bit');

    % Channel Decoding (Hard Decision Decoding - HDD)
    channel_decoding_bits_hard = decode(demodulated_coded_MPSK, n, k, 'hamming/binary');

    % Bit Error Rate (HDD)
    BER_hard_decision_sim(iSNR) = mean(xor(uncoded_bits, channel_decoding_bits_hard));

    % Channel Symbol Error Probability (HDD)
    p = mean(channel_coding_bits(demodulated_coded_MPSK == 0) == 1);

    % Bit Error Rate (Theory - HDD)
    BER_hard_decision_ana(iSNR) = p - p*(1-p)^(n-1); % For (7, 4) Hamming Only

    %     if (M >= 4)
    %         p = 2*qfunc( sqrt(2*EcN0*log2(M))*sin(pi/M) )/log2(M);
    %     else
    %         p = qfunc( sqrt(2*EcN0) );
    %     end

    %% Analytical results
    % Uncoded MPSK
    if (M >= 4) % Larger M, higher deviation from Simulation
        BER_uncoded_ana(iSNR) = 2*qfunc( sqrt(2*EbN0*log2(M))*sin(pi/M) )/log2(M);
    else
        BER_uncoded_ana(iSNR) = qfunc( sqrt(2*EbN0) );
    end

    %% Soft Decision Decoding
    %     cipSoftM  = reshape(real(received_coded_MPSK),n,[]).';
    %
    %     % Measure Correlation, Extract Codeword with Highest Correlation
    %     [~,idx]   = max(cipSoftM*UBPSK.',[],2);
    %     %
    %
    %     MuRev_soft= reshape(dec2bin(idx-1,k).'-'0',[],1);
    %     % Bit Error Rate (SDD)
    %
    %     BER_soft_decision_sim(iSNR) = mean(xor(uncoded_bits,MuRev_soft));

end

% Minimum distance of (7,4) Hamming Code: dmin = 3
BER_Hamming_hard_MATLAB = bercoding(EbN0dB, 'Hamming', 'hard', n);
BER_block_hard_3_MATLAB = bercoding(EbN0dB, 'block', 'hard', n, k, dmin);

BER_uncoded_MATLAB = berawgn(EbN0dB, 'psk', M, 'nondiff');

BER_soft_decision_MATLAB = bercoding(EbN0dB,'block', 'soft', n, k, dmin);

%%

figure;

semilogy(EbN0dB, BER_uncoded_ana, 'b-', ...
    'DisplayName', 'Uncoded (ana.)',...
    'LineWidth',3); hold on;
% semilogy(EbN0dB, BER_uncoded_sim, 'bo:', ...
%     'DisplayName', 'Uncoded (sim.)'); hold on;
semilogy(EbN0dB, BER_uncoded_MATLAB,'b^:', ...
    'DisplayName', 'Uncoded (MATLAB)'); hold on;

semilogy(EbN0dB,BER_hard_decision_ana,'r-', ...
    'DisplayName', '(7,4) Hamming (hard, ana.)',...
    'LineWidth',3); hold on;
% semilogy(EbN0dB,BER_hard_decision_sim,'b^:', ...
%     'DisplayName', '(7,4) Hamming (hard, sim.)'); hold on;
semilogy(EbN0dB,BER_Hamming_hard_MATLAB,'ro:', ...
    'DisplayName', '(7,4) Hamming (hard, MATLAB)'); hold on;
% semilogy(EbN0dB,BER_block_hard_3_MATLAB,'bv:'); hold on;

% semilogy(EbN0dB,BER_soft_decision_sim,'r*:'); hold on;
semilogy(EbN0dB,BER_soft_decision_MATLAB,'g-', ...
    'DisplayName', '(7,4) Hamming (sof, MATLAB)',...
    'LineWidth',3); hold on;

% maximum likelihood (ML) decoding
% https://www.mathworks.com/help/comm/ug/autoencoders-for-wireless-communications.html
hammingBLER = load('codedBLERResults');
simParams.EbNoVec = 0:0.5:8;
semilogy(hammingBLER.simParams.EbNoVec, hammingBLER.hammingML74BLER,'--vg', ...
    'DisplayName', 'AutoEncoder (Machine learning)')

% % maximum likelihood (ML) decoding
% load codedBLERResults.mat
% semilogy(simParams.EbNoVec,hammingML74BLER,'--d')

%
% semilogy(EbN0dB,BER_ML,'b^-'); hold on;

legend('Location','Best');
xlabel('Eb/N0 (dB)');
ylabel('BER');
axis([0 max(EbN0dB) 1e-8 1]);

%% MPSK Constellation
% figure;
% plot(sigRx(1:1e3),'.'); hold on;