viterbi.m

数字通信第四版原书的例程

    if expen_flag
        for j = 1 : length(pre_state)
            jj = pre_state(j) - 1;           % index number - 1 is the state.
            cur_sta = cur_sta_pre(pre_state(j),:)';
            indx_j = (pre_state(j) - 1) * n_std_sta;
            for num_N = 1 : N
                expen_work(num_N) = max(find(tran_prob(1,:) <= code(i, num_N)));
            end;
            for num_k = 1 : K2
                inp = inp_pre(num_k, :)';
                if isempty(C)
                    tran_indx = pre_state(j) + n_std_sta * (num_k - 1);
                    nex_sta = A(tran_indx, :)';
                    out = B(tran_indx, :)';
                else
                    out = rem(C * cur_sta + D * inp,2);
                    nex_sta = rem(A * cur_sta + B * inp, 2);
                end;
                nex_sta_de = bi2de(nex_sta') + 1;
                % find the expense by the transfer probability
                expen_0 = find(out' <= 0.5);
                expen_1 = find(out' > 0.5);
                loca_exp = sum([tran_prob(2,expen_work(expen_0)) 0])...
                                        +sum([tran_prob(3,expen_work(expen_1)) 0]);
                tmp = (nex_sta_de-1)*n_std_sta + pre_state(j);
                if isnan(expense(trace_num, tmp))
                    expense(trace_num, tmp) = loca_exp;
                    solu(trace_num, nex_sta_de + indx_j) = num_k;
                elseif expense(trace_num, tmp) < loca_exp
                    expense(trace_num, tmp) = loca_exp;
                    solu(trace_num, nex_sta_de + indx_j) = num_k;
                end;
            end;
        end;
    else
        for j = 1 : length(pre_state)
            jj = pre_state(j) - 1;           % index number - 1 is the state.
            cur_sta = cur_sta_pre(pre_state(j),:)';
            indx_j = (pre_state(j) - 1) * n_std_sta;
            for num_k = 1 : K2
                inp = inp_pre(num_k, :)';
                if isempty(C)
                    tran_indx = (num_k - 1) * n_std_sta + pre_state(j);
                    nex_sta = A(tran_indx, :)';
                    out = B(tran_indx, :)';
                else
                    out = rem(C * cur_sta + D * inp,2);
                    nex_sta = rem(A * cur_sta + B * inp, 2);
                end;
                nex_sta_de = bi2de(nex_sta') + 1;
                loca_exp = sum(rem(code(i, :) + out', 2));
                tmp = (nex_sta_de-1)*n_std_sta + pre_state(j);
                if isnan(expense(trace_num, tmp))
                    expense(trace_num, tmp) = loca_exp;
                    solu(trace_num, nex_sta_de + indx_j) = num_k;
                elseif expense(trace_num, tmp) > loca_exp
                    expense(trace_num, tmp) = loca_exp;
                    solu(trace_num, nex_sta_de + indx_j) = num_k;
                end;
            end;
        end;
    end;

    aft_state = [];
    for j2 = 1 : n_std_sta
        if max(~isnan(expense(trace_num, [1-n_std_sta : 0] + j2*n_std_sta)))
            aft_state = [aft_state, j2];
        end
    end;

    % go back one step to re-arrange the lines.
    if (plot_flag(1))
        plot_curv_x = [];
        plot_curv_y = [];
        for j = 1 : n_std_sta
            tmp = expense(trace_num, (j-1)*n_std_sta+1: j*n_std_sta);
            if expen_flag
                expen_rec(j) = max(tmp(find(~isnan(tmp))));
            else
                expen_rec(j) = min(tmp(find(~isnan(tmp))));
            end;
            for j2 = 1 : length(tmp)
                if ~isnan(tmp(j2))
                    plot_curv_x = [plot_curv_x, i-1+.1, i, NaN];
                    plot_curv_y = [plot_curv_y, ...
                        flipped(j2), ...
                        flipped(j), NaN];
                end;
            end;
        end;
    end;

    %find out the plot lines:

    if plot_flag(1)
        % plot only; no computation function here.
        plot_curv_x = [get(plot_y, 'XData'), NaN, plot_curv_x];
        plot_curv_y = [get(plot_y, 'YData'), NaN, plot_curv_y];
        for j = 1 : length(aft_state)
            xx=[xx; text(i, flipped(aft_state(j)), num2str(expen_rec(aft_state(j))))];
            set(xx,'Color',[0 1 1]);
        end;
%        plot(plot_curv_x, plot_curv_y, 'y-');
        set(plot_y, 'XData', plot_curv_x, 'YData', plot_curv_y);
        if i < leng
            drawnow
        end;
    end;

    sol = expense(trace_num,:);
    % decision making.
    if i >= leng
        trace_eli = rem(trace_num, leng) + 1;
        % strike out the unnecessary.
        for j_k = 1 : leng - 2
            j_pre = rem(trace_num - j_k - 1 + leng, leng) + 1;
            sol = vitshort(expense(j_pre, :), sol, n_std_sta, expen_flag);
        end;

        tmp = (ones(n_std_sta,1) * expense(trace_eli, [starter+1:n_std_sta:PowPowM]))';
        sol = sol + tmp(:)';

        if expen_flag
            loc_exp =  max(sol(find(~isnan(sol))));   
        else
            loc_exp =  min(sol(find(~isnan(sol))));
        end
        dec = find(sol == loc_exp);
        dec = dec(1);
        dec = rem((dec - 1), n_std_sta);

        inp = de2bi(solu(trace_eli, starter*n_std_sta+dec+1)-1, K);
        if isempty(C)
            tran_indx = pre_state(j) + n_std_sta * (num_k - 1);
            out = B(tran_indx, :)';
        else
            cur_sta = cur_sta_pre(starter+1, :)';
            out = rem(C * cur_sta + D * inp(:),2);
        end;

        msg = [msg; inp];
        codd = [codd; out'];
        [n_msg, m_msg] = size(msg);
        if nargout > 1
            if expen_flag
                % find the expense by the transfer probability
                expen_0 = find(out' <= 0.5);
                expen_1 = find(out' > 0.5);
                loc_exp = sum([tran_prob(2,expen_work(expen_0)) 0])...
                          +sum([tran_prob(3,expen_work(expen_1)) 0]);
            else
                loc_exp = sum(rem(code(n_msg, :) + out', 2));
            end
            expen = [expen; loc_exp];
        end;

        if plot_flag(1)
            plot_curv_x = [get(plot_r, 'XData'), n_msg-1+.1, n_msg, NaN];
            plot_curv_y = [get(plot_r, 'YData'), flipped(starter+1), flipped(dec+1), NaN];
            set(plot_r, 'XData', plot_curv_x, 'YData', plot_curv_y,'LineWidth',2);
            drawnow
        end;

        starter = dec;
    end; %(if i >= leng)
    pre_state = aft_state;
    aft_state = [];
  end;

  for i =  1 : leng-1
    sol = solution;
    trace_eli = rem(trace_num + i, leng) + 1;
    if i < leng-1
        sol(1:n_std_sta) = expense(trace_num, 1:n_std_sta);
        for j_k = 1 : leng - 2 - i
            j_pre = rem(trace_num - j_k - 1 + leng, leng) + 1;
            sol = vitshort(expense(j_pre, :), sol, n_std_sta, expen_flag);
        end;
        tmp = (ones(n_std_sta,1) * expense(trace_eli, [starter+1:n_std_sta:PowPowM]))';
        sol = sol + tmp(:)';
        if expen_flag
            loc_exp =  max(sol(find(~isnan(sol))));   
        else
            loc_exp =  min(sol(find(~isnan(sol))));
        end
        dec = find(sol == loc_exp);
        dec = dec(1);
        dec = rem((dec - 1), n_std_sta);
    else
        dec = 0;
    end;

    inp = de2bi(solu(trace_eli, starter*n_std_sta+dec+1)-1, K);
    cur_sta = de2bi(starter, num_state);
    out = rem(C*cur_sta' + D * inp', 2);
    msg = [msg; inp];
    codd = [codd; out'];
    [n_msg, m_msg] = size(msg);

    if nargout > 1
        if expen_flag
            % find the expense by the transfer probability
            expen_0 = find(out' <= 0.5);
            expen_1 = find(out' > 0.5);
            loc_exp = sum([tran_prob(2,expen_work(expen_0)) 0])...
                      +sum([tran_prob(3,expen_work(expen_1)) 0]);
        else
            loc_exp = sum(rem(code(n_msg, :) + out', 2));
        end
        expen = [expen; loc_exp];
    end;

    if plot_flag(1)
        plot_curv_x = [get(plot_r, 'XData'), n_msg-1+.1, n_msg, NaN];
        plot_curv_y = [get(plot_r, 'YData'), flipped(starter+1), flipped(dec+1), NaN];
        set(plot_r, 'XData', plot_curv_x, 'YData', plot_curv_y,'LineWidth',2);
        drawnow
    end;
    starter = dec;
  end;


  % cut the extra message length
  [n_msg, m_msg] = size(msg);
  msg(n_msg-M+1:n_msg, :) = [];
  if plot_flag(1)
    set(xx,'Color',[0 1 1]);
   hold off
  end;
end;
% end of VITERBI.M