function ratioInfo = vitsimexample(EbNo,M)
% Example of a convolutional encoder and Viterbi decoder simulation
%
%   RATIOINFO = VITSIMEXAMPLE(EBNO,M) is a complete simulation of a binary
%   phase shift keying (BPSK) or Quaternary PSK (QPSK) over an additive white
%   gaussian channel (AWGN) using convolutional encoding and the Viterbi
%   decoding algorithm with hard decision decoding. VITSIMEXAMPLE is an
%   extension of the VITDEMO demonstration and further illustrates how to use
%   the new convolutional trellis generator, POLY2TRELLIS, encoder, CONVENC, and
%   decoder, VITDEC.  It also demonstrates the use of other functions: DMODCE,
%   DDEMODCE, SYMERR, BITERR, RANDSRC and AWGN.
%
%   VITSIMEXAMPLE calculates the union upper bound of the bit error rate for
%   each given Eb/No values, EBNO.  VITSIMEXAMPLE iterates on the vector
%   containing the Eb/No values, EBNO, and simulates the communcations link for
%   each of the EBNO values until 200 errors are generated.  The resulting bit
%   error rate is then plotted along with the performance bounds.
%
%   References
%   "Digital Communications", J. G. Proakis, 1983.
%   "Optimal Decoding of Convolutional Codes", PhD Dissertation,
%    J. P. Oldenwalder, 1970

%   Copyright 1996-2000 The MathWorks, Inc.
%   $Revision: 1.2 $ $Date: 2000/08/31 15:20:17 $

% --- Argument handling

% Default arguments and input error handling
if nargin < 2
    M = 2;
end
if nargin < 1
    error('Eb/No values were not passed');
end

% M can be 2 or 4 for BPSK or QPSK modulation
if (M ~= 2) & (M ~= 4)
    error('M must be 2 or 4')
end

% --- Setup for the simulation

% Define modulation, linear SNR and number of bits per symbol (k).
modulation = 'psk'; k = log2(M);
linSNRpB = 10.^(EbNo(:).*0.1);

% Calculate probability of error based on errors distribution vs. distance for
% rate 1/2 code, K = 7, (Oldenwalder)
codeRate = 1/2; constlen = 7;
dist = [10:2:16];  numErr = [36 211 1404 11633];

% Estimate the expected performance of viterbi algorithm with ideal viterbi
% decoding using union bound for hard decision decoding. (Proakis)
% Bound calculations are only good for BPSK or QPSK with Gray coding.  Use
% a different expresion for expBER for M > 4 or for other modulations.
expBER = 0.5.*erfc(sqrt(linSNRpB*codeRate));
expVitBER = vitproberrdemo(expBER,dist,numErr);
expVitBER(expVitBER>0.9) = NaN;

% Plot the theoretical results for BER.
semilogy(EbNo,expVitBER,'m-');
titleStr = sprintf('Performance of %s (M=%d) and Conv. Code (R=%g, K=%d) with Hard Decision', ...
    upper(modulation), M, codeRate, constlen);
title(titleStr);
xlabel('Eb/No (dB)');   ylabel('BER');
legend('Theoretical BER',0);
grid on;  drawnow;   hold on;

% Set number of symbols per iteration, max and min number of iterations,
% and expected number of symbol errors
symbPerIter = 1024; min_iters = 10; max_iters = 300000000; expSymErrs = 200;

% Initalize sampling times, oversampling rate (=1).
Fd = 1; Fs = 1; N = Fs/Fd;

% Set random number seeds for uniform and Gaussian noise
rand('state', 987654321);  randn('state', 123456789);

% Generate trellis for convolutional encoder & Viterbi decoder for G = 171 133
% Set traceback length = 64
codegen = [171 133]; tblen = 64;
trel = poly2trellis(constlen, codegen);

% Adjust SNR for use by AWGN.  Adjust for coded bits and baseband signals
adjSNR = EbNo - 10*log10(1/codeRate) - 10*log10(N) + 10*log10(k);

% Create Gray encoding and decoding arrays
grayencod	= bitxor([0:M-1],floor([0:M-1]/2));
[dummy graydecod] = sort(grayencod);  graydecod = graydecod - 1;

% --- Outer loop: Iterate on the Eb/No values

% Drive the simulation for each of the SNR values calculated above
for(idx2 = [1:length(EbNo)])

    % Initialize the leftover vector with all zeros to create a decoder
    % 'prehistory' of 0's.
    msg_orig_lo = zeros(tblen,1);   initIdx = 1;

    % Reset the encoder/decoder states and persistent data
    stateEnc = [];  metric = []; stateDec = []; in = [];

    % Initialize the bit error rate and error counters.
    errBitInfo = [];  ratBitInfo = [];
	totErr = 0;       totBits = 0;

    % --- Inner Loop: Iterate until simulation exit criteria is met

    % Exit loop only when minimum number of iterations have completed and  either
    % the number of errors exceeds 'expSymErrs' or the maximum iterations have
    % completed
    idx = 1;
    while(((totErr <= expSymErrs) & idx <= max_iters) | idx <= min_iters)

        % Generate binary random numbers
        msg_orig = randsrc(symbPerIter,1,[0 1]);

        % Convolutionally encode the message, saving encoder state between iterations
        [msg_enc_bi, stateEnc] = convenc(msg_orig(:), trel, stateEnc);

        % Change binary symbols to quaternary symbols and gray encode to ensure
        % one bit transition per symbol
        msg_enc  = bi2de(reshape(msg_enc_bi, size(msg_enc_bi,2)*k,size(msg_enc_bi,1) / k)');
        msg_gr_enc = grayencod(msg_enc+1);

        % Digitally modulate and upsample the signal, shifting signal phase by PI/4
        msg_tx  = dmodce(msg_gr_enc, Fd, Fs, modulation, M);

        % Add Gaussian noise to the signal.
        msg_rx  = awgn(msg_tx, adjSNR(idx2), 'measured', [], 'dB');

        % Demodulate, detect, and downsample the signal, shifting signal phase by PI/4
        msg_gr_demod = ddemodce(msg_rx, Fd, Fs, modulation, M);

        % Gray decode and change demodulated quaternary/binary signals to binary
        msg_demod = graydecod(msg_gr_demod+1);
        msg_demod_bi = de2bi(msg_demod,k)';

        % Use the Viterbi algorithm to decode the recieved signal.  Save the
        % trellis and metric states from iteration to iteration in
        % 'metric', 'stateDec', and 'in'.
        [msg_dec, metric, stateDec, in] = vitdec(msg_demod_bi(:), trel, ...
			tblen, 'cont', 'hard', metric, stateDec, in);

        % Add leftover symbols from last iteration to those from this iteration
        msg_orig_w_lo = [msg_orig_lo; msg_orig];

        % Compare the input and outputs to determine BER and error count
        size_msg_dec = size(msg_dec,1) - initIdx + 1;
        [errBitInfo ratBitInfo] = biterr(msg_dec(initIdx:end), ...
            msg_orig_w_lo(initIdx:size_msg_dec+initIdx-1));

		% Acculumulate bit count and bit error statistics after each iteration
		totErr = totErr + errBitInfo;
		totBits = totBits + size_msg_dec;

        % Save leftover bits from this iteration for the next iteration
        msg_orig_lo = msg_orig_w_lo(end-tblen+1:end);

        % Increment for next iteration, set the initial index for comparison in
        % subsequent iterations and pause to allow <ctrl-C> break
        idx = idx + 1;   initIdx = 1+tblen;    pause(.1);
    end

    % Compute the BER for all the iterations together
    errorsInfo(idx2) = totErr;
    ratioInfo(idx2)  = totErr/totBits;

    % Plot the simulated BER results for the current SNR per bit
    semilogy(EbNo(idx2),ratioInfo(idx2),'bo');
    legend('Theoretical BER','Simulated BER',0);
	drawnow;
end
hold off;