sviterbi.m

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

function [sys, x0, str, ts] = sviterbi(t, x, u, flag, tran_func, leng, tran_prob, plot_flag);
%SVITERBI SIMULINK file for convolution decoding using viterbi algorithm.
%       This file is designed to be used in a SIMULINK S-Function block.
%       The function requires the system inputs
%       code_param = [N, K, M, num_state, num_std_sta]; % prepared by simviter
%       tran_func  = [A, B; C, D];                      % prepared by simviter
%       leng                                            % memory length
%       tran_prob                                       % transfer probability
%          % when it is not a three row matrix, take the code as hard decision one.
%       plot_flag                                       % should plot or not.
%          % when it is not a positive integer, don't plot.
%          % when it is a positive integer, keep the plot have the time length
%          % as required.
%       This function keeps a number of discrete-time states:
%       figure number. -Inf: to be opened; 0 not need to plot; positive: handle
%       figure posit.  Current figure position.
%       trace_num.     current trace number.
%       trace_flag.    flag to indicate that computation is not under leng.
%       stater.        variable in keeping the very beinging of state.
%       trace.         a LENG-by-NUM_STD_STA matrix storage the traces.
%       solut.         a LENG-by-NUM_STD_STA matrix storage the possible msg.
%       expense.       a 2-by-NUM_STD_STA matrix storage the expense.
%

%       Wes Wang
%       Copyright (c) 1995-96 by The MathWorks, Inc.

if (flag == 2) % refresh discrete-time states
%  t2 = clock;
  % the major processing routine.
  if u(length(u)) < .2
    % in the case of no signal, no processing.
    return;
  end;

  % otherwise, there is a signal, have to process.
  u = u(:)';

  % initial parameters.
  [A, B, C, D, N, K, M] = gen2abcd(tran_func);
  num_state = M;
  n_std_sta = 2^M;
  PowPowM = n_std_sta^2;
  K2 = 2^K;
  inp_pre = de2bi([0:K2-1]', K);
  cur_sta_pre = de2bi([0:n_std_sta-1], M);

  %%%%%begin plot related %%%%%
  if x(1) < 0
    % initialize the figure.
    [sl_name, block] = get_param;

    [n_b, m_b]= size(block);
    if n_b < 1
      error('Cannot delete block while simulating')
    elseif n_b > 1
      error('Something wrong in get_param, You don''t have the current SIMULINK')
    end;

    % findout if the graphics window exist
    ind = find(sl_name == setstr(10));
    dash = '_';
    sl_name(ind)=dash(ones(size(ind)));

    % in the case of the figure exists
    Figures = get(0, 'Child');
    new_figure = 1;
    i = 1;
    while ((new_figure) & (i <= length(Figures)))
      if strcmp(get(Figures(i), 'Type'), 'figure')
        if strcmp(get(Figures(i), 'Name'), sl_name)
          h_fig = Figures(i);
          handles = get(h_fig,'UserData');
          new_figure = 0;
          h_axis = handles(1);
          set(0,'CurrentF',h_fig);
          set(h_axis,'Xlim', [0, plot_flag-1], 'Ylim', [0, n_std_sta-1]);
          delete(get(h_axis,'child'));
          handles = handles(1);
        end;
      end;
      i = i + 1;
    end;
    if new_figure
      h_fig = figure;
      set(h_fig, 'Name', sl_name);
      handles = axes('Xlim',[0, plot_flag-1], 'Ylim',[0, n_std_sta-1]);
    end;
    handles = get(h_fig, 'Child');
    set(handles, 'Visible', 'Off');
    set(handles(1),'NextPlot','add');
    set(h_fig,'NextPlot','add', 'Clipping', 'off');
    handles(2) = plot(0,0,'y-','Erasemode','normal');
    handles(3) = plot(0,0,'r-','Erasemode','normal','LineWidth',2);
    handles(4) = plot(-10,-10, 'go', 'Erasemode','normal');
    for i = 1 : n_std_sta
      handles(i+4) = text(0,0,'0');
      set(handles(i+4), 'Color',[0,1 1]);
    end;
    for i = 0 : n_std_sta - 1
      handles(5+n_std_sta+i) = text(-plot_flag/8, i, ['State ', num2str(bi2de(fliplr(de2bi(i, num_state))))]);
    end;
    handles(5+2*n_std_sta) = text(plot_flag/2.1, - n_std_sta/15, 'Time');
    set(h_fig, 'UserData',handles);
    set(h_axis,'NextPlot','new');
    set(h_fig,'NextPlot','new');
    sys(1) = h_fig;
    x(1) = h_fig;
  end;% 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;
  %%%%% end plot related %%%%%

  [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.
    expen_work = zeros(1, N); % work space.
    tran_prob(2:3, :) = log10(tran_prob(2:3, :));
  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.
  starter = x(5);
  solu = zeros(leng, PowPowM);
  expense = solu;
  code = zeros(leng, N);
  loc_tmp = 6;
  expense(:) = x(loc_tmp : leng*PowPowM + loc_tmp - 1);
  solu(:) = x(loc_tmp + leng * PowPowM : 2 * leng * PowPowM + loc_tmp - 1);
  code(:) = x(loc_tmp + 2 * leng * PowPowM : loc_tmp + 2 * leng * PowPowM -1 + leng * N);

  %%%%% begin plot related %%%%%
  fig_position = x(2) + 1;
  if (x(1) > 0) & (rem(fig_position-leng, plot_flag - leng) == 0) & (fig_position >= plot_flag) & plot_flag_test
    set(handles(1), 'XLim', [fig_position-leng, fig_position + plot_flag - leng]);
    plot_curv_x = get(handles(2), 'XData');
    plot_curv_y = get(handles(2), 'YData');
    index = find(plot_curv_x == -(fig_position - leng));
%    index = find(plot_curv_x <= (fig_position - leng -1));
    leng_x = length(plot_curv_x);
    if isempty(index)
      % no move
    else
%      index = index(length(index));
      plot_curv_x = plot_curv_x(index(1)+2: leng_x);
      plot_curv_y = plot_curv_y(index(1)+2: leng_x);
    end;
    set(handles(2), 'XData', plot_curv_x, 'YData', plot_curv_y);
    for i = 0:n_std_sta - 1
      tmp = get(handles(5+n_std_sta+i), 'Position');
      tmp(1) = fig_position-leng-plot_flag/8;
      set(handles(5+n_std_sta+i), 'Position', tmp);
    end;
    tmp = get(handles(5+2*n_std_sta), 'Position');
    tmp(1) = fig_position-leng+plot_flag/2.1;
    set(handles(5+2*n_std_sta), 'Position', tmp);
    plot_curv_x = get(handles(4), 'XData');
    plot_curv_y = get(handles(4), 'YData');
    index = find(plot_curv_x <= fig_position-leng);
    index = max([index(length(index)) 0]) + 1;
    leng_x = length(plot_curv_x);
    plot_curv_x = plot_curv_x(index : leng_x);
    plot_curv_y = plot_curv_y(index : leng_x);
    set(handles(4), 'XData', plot_curv_x, 'YData', plot_curv_y);
%    drawnow;
  end;

  trace_num = x(3) + 1;
  trace_flag = x(4);
  code(trace_num, :) = u(1:length(u)-1);

  if (trace_flag == 0) & (trace_num == leng)
    trace_flag = 1;
  end;
  %%%%% end plot related %%%%%

  trace_pre = rem(trace_num - 2 + leng, leng) + 1;

  % fill up trace and solut
  if (trace_flag == 0) & (trace_num == 1)
    pre_state = starter + 1;
  else
    pre_state = [];
    for j2 = 1 : n_std_sta
      if max(~isnan(expense(trace_pre, [1-n_std_sta : 0] + j2*n_std_sta)))
          pre_state = [pre_state, j2];
      end;
    end;
  end;

  %%%%% begin plot related %%%%%
  if (x(1) > 0) & plot_flag_test
    plot_curv_x = [get(handles(2), 'XData'), NaN, -fig_position, NaN];
    plot_curv_y = [get(handles(2), 'YData'), NaN, -1000, NaN];
    dot_x = get(handles(4), 'XData');
    dot_y = get(handles(4), 'YData');
  end;
  %%%%% end plot related %%%%%

  expense(trace_num, :) = expense(trace_num,:) * NaN;
  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(trace_num, num_N)));
      end;
      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;
        % 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;