function [err, rat] = sttc4Psk(stMap, nextStMap, bin, nrState, G, N, ch, snrInDb)
% Space-time trellis encoder/decoder for QPSK modulation
% 2 bits per second per Hertz
%
% 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 and AWGN
s = awgn((s1 + s2), snrInDb);

% Received signal length
inLen = length(s);

% Initializations 
prevSt = zeros(nrState, inLen);
metric = zeros(1, nrState);

% Viterbi decoding; find the surviving paths over time interval inLen
for k = 2:inLen - 1,      
   for l = 1:stPerMAry,            
      for m = 1:mAry,            
         minDist = 1e6;
         for n = 0:inBranch - 1,
                                   
            % Find the distance of the received vector over all possible branches
            dist = (ch(1) * constel(stMap(stPerMAry*n + l, 1) + 1) + ...
               ch(2) * constel(nextStMap(stPerMAry*n + l, m) + 1) - s(k));
            % Magnitude    
            dist = dist * conj(dist);
               
            % Determine the minimum distance and the suviving paths
            if dist + metric(stPerMAry*n + l) < minDist
               minDist = dist + metric(stPerMAry*n + l);
               nextMetric(m + mAry*(l - 1)) = minDist;
               prevSt(m + mAry*(l - 1), k) = stPerMAry*n + l - 1; 
            end  % if
         end  % for n    
      end  % for m
   end  % for l
   metric = nextMetric;
end  % for k

% Trace back to find the estimated states
xTmp = zeros(1, inLen + 1);

for k = inLen:-1:1,
    xTmp(k) = prevSt(xTmp(k + 1) + 1, k);
    decodedSt(k) = 4*xTmp(k) + mod(xTmp(k + 1), 4) + 1;
end

% Find the estimated bits sequence
for i = 1:inLen - maxDelay
   zDec(2*i - 1) = bin(decodedSt(i), col);
   zDec(2*i) = bin(decodedSt(i), col - 1);
end

% Count the frame error 
err = 0;
[num, rat] = biterr(zDec, dSource);
err = ceil(rat);