function [alpha, beta] = f_b(a, b)

T = size(b, 2);
Q = size(a,1);
a(Q-1,Q-1)=1;

alpha = zeros(Q,T);
beta = zeros(Q,T);

%FORWARD
t = 1;
alpha(:,t) = a(1,:)*b(:,1);
alpha(1,1) = 1;
% alpha(:,t) = zeros(Q,1);
% alpha(2,1) = 1;
for t = 2 : 1 : T
    alpha(:,t) = b(:,t) .* (a'*alpha(:,t-1));
    c=sum(alpha(:,t));
    d = c + (c==0);
    alpha(:,t) = alpha(:,t) / d;
end

%BACKWARD
beta(:,T) = zeros(Q,1);
beta(2:(Q-1),T)=a(2:(Q-1),Q);
%beta(:,T) = a(:,Q);
for t = T-1 : -1 : 1
    beta(:,t) = zeros(Q,1);
    beta(2:(Q-1),t) = a(2:(Q-1),2:(Q-1)) * (b(2:(Q-1),t+1).*beta(2:(Q-1),t+1));
    c=sum(beta(:,t));
    d = c + (c==0);
    beta(:,t) = beta(:,t) / d;
end
