sviterbi.m

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

          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 %expen_flag = 0
%    t0 = clock;
    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 = pre_state(j) + (num_k -1) * K2;
            tran_indx = pre_state(j) + (num_k -1) * n_std_sta;
            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(trace_num,:) + 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;
%    disp(etime(clock, t0))
  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;

  %%%%% begin plot related %%%%%
  if (x(1) > 0) & plot_flag_test
    % plot only; no computation function here.
    flipped = zeros(1, n_std_sta);
    for i = 1 : n_std_sta 
      flipped(i) =  bi2de(fliplr(de2bi(i-1, M)));
    end;
    for j2 = 1 : length(pre_state)
      jj = pre_state(j2) - 1;
      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, fig_position-1+.1, fig_position, NaN];
            plot_curv_y = [plot_curv_y, ...
                           flipped(j2), ...
                           flipped(j), NaN];
          end;
        end;
      end;
    end;
    for j2 = 1 : length(aft_state)
      jj = aft_state(j2) - 1;
      set(handles(jj+5), 'Position', [fig_position, bi2de(fliplr(de2bi(jj,num_state))), 0],...
           'String',  num2str(expen_rec(jj+1)));
      tmp = get(handles(jj+5), 'Position');
      dot_x = [dot_x, tmp(1)];
      dot_y = [dot_y, tmp(2)];
    end;
    set(handles(2), 'XData', plot_curv_x, 'YData', plot_curv_y);
    set(handles(4), 'XData', dot_x, 'YData', dot_y);
    if ~trace_flag
      drawnow;
    end;
  end; 
    %%%%% end plot related %%%%%

  % decision making.
%  t0 = clock;
  if trace_flag
    sol = expense(trace_num,:);
    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);
    num_k = solu(trace_eli, starter*n_std_sta+dec+1);
    inp = inp_pre(num_k, :)';
    if isempty(C)
%        tran_indx = pre_state(j) + (num_k -1) * K2;
        tran_indx = starter + 1 + (num_k -1) * n_std_sta;
        out = B(tran_indx, :)';
    else
        cur_sta = cur_sta_pre(starter+1, :)';
        out = rem(C * cur_sta + D * inp,2);
    end;

    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(trace_eli, :) + out', 2));
    end

    %%%%%begin plot related %%%%%
    if plot_flag & plot_flag_test
      plot_curv_x = [get(handles(3), 'XData'), NaN];
      plot_curv_y = [get(handles(3), 'YData'), NaN];
      if length(plot_curv_x) > 3 * plot_flag
        leng_x = length(plot_curv_x);
        limi    = leng_x - 3 * plot_flag + 1;
        plot_curv_x = plot_curv_x(limi:leng_x);
        plot_curv_y = plot_curv_y(limi:leng_x);
      end;

      set(handles(3), 'XData', [plot_curv_x, fig_position-leng+.1, fig_position-leng+1],...
                      'YData', [plot_curv_y, bi2de(fliplr(de2bi(starter, num_state))),...
                       bi2de(fliplr(de2bi(dec, num_state)))]);
      drawnow
    end;
    %%%%% end plot related %%%%%
    starter = dec;
    output = [inp(:); loca_exp];
  else %(if trace_flag)
    output = zeros(K+1, 1);
    starter = 0;
  end; %(if trace_flag)
%  disp(['decision etime ', num2str(etime(clock, t0))])      

  trace_num = rem(trace_num, leng);
  sys = [x(1); fig_position; trace_num; trace_flag; starter; expense(:); solu(:); code(:); output(:)];
%  disp(['total etime ', num2str(etime(clock, t2))])
elseif flag == 3 % output
  if u(length(u)) > 0.2
    K = tran_func(2, size(tran_func, 2));
    len_x = length(x);
    sys = x(len_x - K: len_x);
  end;
elseif flag == 0
  [A, B, C, D, N, K, M] = gen2abcd(tran_func);
  num_state = M;
  n_std_sta = 2^M;
  PowPowM = n_std_sta^2;

  [n_tran_prob, m_tran_prob] = size(tran_prob);
  if n_tran_prob == 3
    expen_flag = 1;     % expense flag == 0; BSC code. == 1, with TRAN_PROB.
  else
    expen_flag = 0;
  end;

  % veraible to record the weight trace back to the first state.
  % the column number is the state number - 1.
  % the first line is the previous. The second line is the current.
  expense = ones(leng, PowPowM) * NaN;
  expense(leng, [1:n_std_sta]) = zeros(1, n_std_sta);
  solu = zeros(leng, PowPowM);

  % the column number is the state + 1.
  % the contents is the state it transfered from.
  % the starter is always single number.
  starter = 0;

  % the solution in each of the transfer as above.
  % ith row is the transition input (msg) from (i-1)th row of trace to ith row
  % of trace. When i = 1, the transfer is from starter to trace.

  if plot_flag > 0
    x0 = -Inf;
  else
    x0 = 0;
  end;
  code = zeros(leng, N);
  y = zeros(K+1, 1);

  % add output storage.
  x0 = [x0; 0; 0; 0; starter; expense(:); solu(:); code(:); y];
%        |  |  |  |      |       |        |        \- output
%        \  \  \  \      \       \        \ decode
%         \  \  \  \      \start  \-weight, expense
%          \  \  \  \-trace_flag
%           \  \  \-trace_num
%            \  \-fig_position
%             \-figure handle
  sys = [0; length(x0); K+1; N+1; 0; 0; 1];
  ts = [-1, 0];
elseif flag == 9
  if x(1) > 0
    plot_flag_test = get(0,'child');
    if isempty(plot_flag_test)
      plot_flag_test = 0;
    elseif isempty(find(plot_flag_test==x(1)))
      plot_flag_test = 0;
    else
      plot_flag_test = 1;
      handles = get(x(1), 'UserData');
      if length(handles) < (5+2*n_std_sta)
        plot_flag_test = 0;
      end;
    end;
    if plot_flag_test
      plot_curv_x = get(handles(4), 'XData');
      plot_curv_y = get(handles(4), 'YData');
      set(handles(4), 'Erasemode','Normal',...
                     'XData', plot_curv_x,...
                    'YData', plot_curv_y);
      drawnow
    end;
  end;
  sys = [];
else
  sys = [];
end;