clear; close all; clc
rng default
M = 4;                 % Modulation order
k = log2(M);            % Bits per symbol
EbNoVec = (0:6)';       % Eb/No values (dB)
numSymPerFrame = 300000;   % Number of QAM symbols per frame
modul = comm.PSKModulator(M, 'BitInput', true);
berEstSoft = zeros(size(EbNoVec)); 
trellis = poly2trellis(7,[171 133]);
tbl = 96;
rate = 3/4;
spect = distspec(trellis);
encoders = comm.ConvolutionalEncoder(trellis,...
    'PuncturePatternSource', 'Property', 'PuncturePattern', [1; 1; 0; 1; 0; 1]);
decoders = comm.ViterbiDecoder(trellis,'TracebackDepth',tbl,...
    'TerminationMethod','Continuous','InputFormat','Unquantized',...
    'PuncturePatternSource', 'Property', 'PuncturePattern', [1; 1; 0; 1; 0; 1]);
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, numErrsHard, numBits] = deal(0);
    
    while numErrsSoft < 100 && numBits < 1e7
        % Generate binary data and convert to symbols
        dataIn = randi([0 1], numSymPerFrame*k, 1);
        
        % Convolutionally encode the data
        dataEnc = step(encoders, dataIn);
        % QAM modulate
        txSig = step(modul, dataEnc);
        % Pass through AWGN channel
        rxSig = awgn(txSig, snrdB, 'measured');
        
        % Demodulate the noisy signal using hard decision (bit) and
        % soft decision (approximate LLR) approaches.
        
        demods = comm.PSKDemodulator(M, 'BitOutput', true, ...
        'DecisionMethod', 'Approximate log-likelihood ratio', 'VarianceSource', 'Property', 'Variance', noiseVar);
        rxDataSoft = step(demods, rxSig);
    
        % Viterbi decode the demodulated data
        dataSoft = step(decoders, rxDataSoft);
               
        % Calculate the number of bit errors in the frame. Adjust for the
        % decoding delay, which is equal to the traceback depth.
        numErrsInFrameSoft = biterr(dataIn(1:end-tbl), dataSoft(tbl+1:end));
        
        % Increment the error and bit counters
        numErrsSoft = numErrsSoft + numErrsInFrameSoft;
        numBits = numBits + numSymPerFrame*k;
    end
    
    % Estimate the BER for both methods
    berEstSoft(n) = numErrsSoft/numBits;
end
%% Theoretical curves
spect = distspec(trellis, 7);
soft_bertool = bercoding(EbNoVec,'conv','soft',1/2,spect); % BER bound
figure(1)
semilogy(EbNoVec, soft_bertool.','-o',EbNoVec,berEstSoft.','-o', 'LineWidth', 1.5)
grid on
hold on
legend('1/2 (theory)','3/4 (simulation)','location','best')
grid on
xlabel('Eb/No (dB)')
ylabel('Bit Error Rate')