function [err, rat] = sttc4PskUseMex(stMap, nextStMap, bin, nrState, G, N, ch, snrInDb)
% Space-time trellis encoder/decoder for QPSK modulation 
% 2 bits per second per Hertz
% This function uses Matlab's C-MEX for decoding part
%
% stMap     : Current states trellis 
% nextStMap : Next states trellis
% bin       : Decimal to binary map matrix
% nrState   : Number of state
% G         : Generator matrix
% N         : Number of data (bits)
% ch        : Fading channel coefficients
% snrInDb   : SNR in dB
%
% err       : Number of frame error
% rat       : Ratio of bits error
%
%
% Copyright Bagawan S. Nugroho, 2001-2006


% Initialize the parameters
mAry = 4;                     % m-ary of PSK
inBranch = mAry;              % Number of incoming branch in the trellis
stPerMAry = nrState/mAry;     % Number of state divided by m-ary
powOf2 = log2(nrState);       % 2^power equals to number of state

% Other paramaters based on the generator matrix
[row col] = size(G);
maxDelay = ceil(col/2);

% QPSK constellation map
constel = [1 1i -1 -1i];

% Generate random binary data within the frame
dSource = round(rand(1, N));

% Interlace bits
tempData = reshape(dSource, 2, N/2);

% Add leading and trailing zero padding accordingly
% (The formula of padding is obtained by experiments) 
leadPad = zeros(1, maxDelay - 1);
tailPad = zeros(1, maxDelay);

a = [leadPad tempData(1,:) tailPad];
b = [leadPad tempData(2,:) tailPad];

% Length of the new bits order
bitsLen = length(a);

% Bits arrangement to built the trellis
D = cat(1, a(1, maxDelay:bitsLen), b(1, maxDelay:bitsLen)); 

for i = (maxDelay - 1):-1:1
   D = cat(1, D, a(1, i:bitsLen + i - maxDelay), b(1, i:bitsLen + i - maxDelay));   
end

% x(1,:) the is current states; x(2,:) is the next states
x = zeros(2, bitsLen + 1 - maxDelay);

for i = 1:col
   x(1,:) = mod((x(1,:) + D(i,:).*G(1, col - i + 1)), 4);
   x(2,:) = mod((x(2,:) + D(i,:).*G(2, col - i + 1)), 4);
end

% QPSK constellation mapping
s1 = constel(x(1,:) + 1).* ch(1);
s2 = constel(x(2,:) + 1).* ch(2);

% Receiver input with AWGN
s = awgn((s1 + s2), snrInDb);

% Received signal length
inLen = length(s);

% Viterbi decoding - Using Matlab's C-MEX 
zDec = stViterbi(s, stMap, nextStMap, nrState, bin, ch);

% Count the frame error; bitErrFast is C-MEX 
err = 0;
[num, rat] = bitErrFast(zDec, dSource);
err = ceil(rat);