clear; close all; clc
rng default
M = [4, 8, 16, 64]; % Modulation order
EbNoVec = (0:5)'; % Eb/No values (dB)
numSymPerFrame = 100000; % Number of QAM symbols per frame
berEstSoft = zeros(size(EbNoVec)); 
trellis = poly2trellis(7,[171 133]);
tbl = 32;
rate = 1/2;
decoders = comm.ViterbiDecoder(trellis,'TracebackDepth',tbl,...
'TerminationMethod','Continuous','InputFormat','Unquantized');
for m = 1:length(M)
    k = log2(M(m)); % Bits per symbol
    if M(m) <= 8
        modul = comm.PSKModulator(M(m), 'BitInput', true);
    end
    for n = 1:length(EbNoVec)
        % Convert Eb/No to SNR
        snrdB = EbNoVec(n) + 10*log10(k*rate);
        % Noise variance calculation for unity average signal power.
        noiseVar = 10.^(-snrdB/10);
        % Reset the error and bit counters
        [numErrsSoft_exact, numErrsHard, numBits] = deal(0);
        [numErrsSoft_approx, numErrsHard, numBits] = deal(0);
        
        while (numErrsSoft_exact < 100 OR numErrsSoft_approx < 100)... 
            && numBits < 1e8
            % Generate binary data and convert to symbols
            dataIn = randi([0 1], numSymPerFrame*k, 1);
            
            % Convolutionally encode the data
            dataEnc = convenc(dataIn, trellis);
            
            % QAM modulate
            if M(m) <= 8
                txSig = step(modul, dataEnc);
            else
                txSig = qammod(dataEnc, M(m), 'InputType','bit',...
                               'UnitAveragePower',true);
            end
            % Pass through AWGN channel
            rxSig = awgn(txSig, snrdB, 'measured');
            
            % Demodulate the noisy signal using hard decision (bit) and
            % soft decision (approximate LLR) approaches.       
            if M(m) <= 8
                demods_approx = comm.PSKDemodulator(M(m), ...
                    'BitOutput', true, ...
                    'DecisionMethod', ...
                    'Approximate log-likelihood ratio',...
                    'VarianceSource', 'Property', 'Variance', noiseVar);
                demods_exact = comm.PSKDemodulator(M(m), ...
                    'BitOutput', true, ...
                    'DecisionMethod', 'Log-likelihood ratio',...
                    'VarianceSource', 'Property', 'Variance', noiseVar);
                rxDataSoft_exact = step(demods_exact, rxSig);
                rxDataSoft_approx = step(demods_approx, rxSig);
            else 
                
                rxDataSoft_exact = qamdemod(rxSig, M(m), ...
                    'OutputType','llr', ...
                    'UnitAveragePower',true,'NoiseVariance',noiseVar);
                rxDataSoft_approx = qamdemod(rxSig, M(m), ...
                    'OutputType','approxllr', ...
                    'UnitAveragePower',true,'NoiseVariance',noiseVar);
            end
            % Viterbi decode the demodulated data
            dataSoft_exact  = step(decoders, rxDataSoft_exact );
            dataSoft_approx = step(decoders, rxDataSoft_approx);
            
            % Calculate the number of bit errors in the frame. 
            % Adjust for the decoding delay, 
            % which is equal to the traceback depth.
            numErrsInFrameSoft_exact = biterr(dataIn(1:end-tbl), ...
                dataSoft_exact(tbl+1:end));
            numErrsInFrameSoft_approx = biterr(dataIn(1:end-tbl), ...
                dataSoft_approx(tbl+1:end));
            
            % Increment the error and bit counters
            numErrsSoft_exact = numErrsSoft_exact + ...
                                numErrsInFrameSoft_exact;
            numErrsSoft_approx = numErrsSoft_approx + ...
                                 numErrsInFrameSoft_approx;
            
            numBits = numBits + numSymPerFrame*k;
        end
        
        % Estimate the BER for both methods
        berEstSoft_exact(n, m) = numErrsSoft_exact/numBits;
        berEstSoft_approx(n, m) = numErrsSoft_approx/numBits;
    end
end
semilogy(EbNoVec, berEstSoft_exact(:, 1),'r-o', ...
         EbNoVec, berEstSoft_exact(:, 2),'k-o',...
         EbNoVec, berEstSoft_exact(:, 3),'b-o', ...
         EbNoVec, berEstSoft_exact(:, 4),'c-o',...
         EbNoVec, berEstSoft_approx(:, 1),'r->', ...
         EbNoVec, berEstSoft_approx(:, 2),'k->',...
         EbNoVec, berEstSoft_approx(:, 3),'b->', ...
         EbNoVec, berEstSoft_approx(:, 4),'c->','LineWidth', 1.5)
hold on
legend('QPSK, Exact LLR', ...
       '8PSK, Exact LLR', ...
       '16-QAM, Exact LLR', ...
       '64-QAM, Exact LLR',...
       'QPSK, Approx. LLR', ...
       '8PSK, Approx. LLR', ...
       '16-QAM, Approx. LLR', ...
       '64-QAM, Approx. LLR', ...
       'location','best')
grid
title('Convolutional codes 1/2, AWGN')
xlabel('Eb/No (dB)')
ylabel('Bit Error Rate')