function [I,dmin] = va(v,f,A,SpS)
% va - Equalization using the Viterbi algorithm 
%
% Given an observed sequence v, the tap gains of the FIR filter f, and
% the alphabet of symbols A, i.e., a vector containing al possible input
% symbols, the most likely sequence of input symbols is computed using
% the Viterbi algorithm. The minimum distance dis also returned. The
% parameter SpS denotes the number of samples per symbol.
%
% I = va( v, f, A, SpS)

% for backwards compatibility, if SpS is not supplied it defaults to 1
if (nargin == 3)
    SpS = 1;
end
% pad filter with zeros so that its length is a multiple of SpS
Lf = ceil(length(f)/SpS)*SpS;        % Filter order
f = [f zeros(1, Lf-length(f))];
% pad input v with zeros so that its length is a multiple of SpS minus
% one, i.e., the padded length of v and f are consistent with the original
% number of symbol: the length of v is related to the length of the input
% and the length of f via: length(v) = length(symbols)*SpS + length(f) - 1.
N = ceil((length(v)+1)/SpS)*SpS - 1; 
v = [v zeros(1, N-length(v))];

% Initialization
M = length(A);            % Alphabet size
L = Lf/SpS-1;             % Number of "additional" symbol periods spanned by ISI
N = (length(v)-Lf+1)/SpS;    % Number of input symbols

short = zeros(M^L,N);     % Create storage for large arrays
survivor = short;

% Set up tables defining for each possible transition:
%
%	- Initial state number, init (element of {0,1,...,M^L-1})
%	- Final state number, final (element of {0,1,...,M^L-1})
%     	- new symbol, ns (element of A)
%	- reconstructed signal value, v_hat
%	- register contents after transition
%
%
% Enumerate all possible transition and assign corresponding input symbols
% through functions dec2baseM and baseM2dec
%

% initialize vectors
ns    = zeros( 1,M^(L+1) );
init  = zeros( 1,M^(L+1) );
final = zeros( 1,M^(L+1) );
v_hat = zeros( M^(L+1), SpS );
contents = zeros( M^(L+1), L );

for i = 1:M^(L+1)
	b = dec2baseM(i-1,M,L+1); 	% Find the base M representation of i-1
    s = A(b + 1); 	            % corresponding symbols
    
	ns(i) = s(1);			% new symbol
	init(i) = baseM2dec(b(2:L+1),M,L)+1; 	% ID number of initial state
	final(i) = baseM2dec(b(1:L),M,L)+1;  	% ID number of final state
    for j=1:SpS
	    v_hat(i,j) = s*f(j:SpS:end).';		% reconstructed signal, SpS samples per symbol
    end
	contents(final(i),1:L) = s(1:L); 	% register contents 
					% oldest bit is in L-th position
end
%
% The registers are initially empty. The first L symbols load the registers
% already producing outputs. Thus, we can compute the distance for each of 
% the M^L possible sequences and assign each of these to the proper state.
%
last = zeros(M^L, L*SpS);
dist = zeros(1, M^L);

for i = 1:M^L
	s   = fliplr(contents(i,1:L));
	out = tdl(s,f,0,SpS);
	first = out(1:L*SpS);
	last(i,1:L*SpS) = out(SpS*L+1:2*L*SpS);
	dist(i) = norm(first-v(1:L*SpS))^2;
	survivor(i,1:L) = s;
end

%
% The main recursion of the Viterbi Algorithm
%
for i = L+1:N
	new_dist = Inf*ones(1,M^L);   	% Initialize with large values

	for t = 1:M^(L+1)		% Loop over all transitions
		update = dist(init(t)) + norm(v((i-1)*SpS+1:i*SpS)-v_hat(t,:))^2; 
								% Compute
					% distance resulting from transition
		if new_dist(final(t)) > update
			new_dist(final(t))    = update;  
			short(final(t),1:i-1) = survivor(init(t),1:i-1); 
			short(final(t),i)     = ns(t);
		end				      
	end %for t ...
	dist     = new_dist;
	survivor = short;
end % for i ...
%
% Now, there are no new input symbols. However, the symbols in the register
% still produce observations. The reconstructed values for these observations
% depend on the register contents at time N and have been computed earlier
% (last(i,:)).
%
for i = 1:M^L
	dist(i) = dist(i) + norm(last(i,1:L*SpS)-v(N*SpS+1:(N+L)*SpS))^2;
end

% Find the shortest path
[dmin,ind] = min(dist);
dmin = sqrt(dmin);

I = survivor(ind,1:N);
