
function [s,v] = viterbi(x,hmm);

ns= length(hmm.p0);
n = length(x);

v = viterbi_values(x,hmm);

s = zeros(1,n);
s(n) = argmax(v(:,n));
for t=n-1:-1:1,
  s(t) = argmax(v(:,t)+log(hmm.Pt(:,s(t+1))));
end;

% ------------------------------
% viterbi values

function [v] = viterbi_values(x,hmm)

n = length(x);
ns= length(hmm.p0);
v = zeros(ns,n);

v(:,1) = log(hmm.p0' .* hmm.Po(:,x(1)));
for t=2:n,
  %v(:,t) = max( v(:,t-1)*ones(1,ns)+log(hmm.Pt), [], 2)+log(hmm.Po(:,x(t)));
  v(:,t) = max( v(:,t-1)*ones(1,ns)+log(hmm.Pt), [], 1)'+log(hmm.Po(:,x(t)));
end;

% ------------------------------

function [s] = argmax(v)

[vt,I] = sort(-v);
s = I(1);
