sviternd.m

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

function [sys, x0, str, ts] = sviternd(t, x, u, flag, code_param, tran_func, leng, tran_prob, plot_flag);
%SVITERND 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.

% expense_flag == 0 when tran_prob IS NOT a three row matrix.
% In this case, all variable in this functions are integers.
% expense_flag == 0 when tran_prob IS a three row matrix.
% In this case, tran_prob, expense and some expense related
% variables are float, such as local_expense, right_exp, tmp etc.
% All others are interfer.
% The calculation C * cur_sta + D * inp are binary computation.

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

    % otherwise, there is a signal, have to process.
    %x = [0; 0; 0; 0; starter; trace(:); solut(:); store_code; expense(:); y];
    %     |  |  |  |      |       |        |           |         |         \- output
    %     \  \  \  \      \       \        \           \         \- weight, expense
    %      \  \  \  \      \start  \-state  \-decode    \- code_storage
    %       \  \  \  \-trace_flag
    %        \  \  \-trace_num
    %         \  \-fig_position, keep for plot m-files.
    %          \-figure handle, keep for plot m-files.


    % pre-process, assign parameters.
    % the last input is the trigger signal. The others are the codeword input
    u = u(:)';
    % codeword length
    N = code_param(1);
    % message length
    K = code_param(2);
    % transfer function memory length
    M = code_param(3);
    % state number
    num_state = code_param(4);
    % STD state number
    num_std_sta = code_param(5);

    % pre-process, assign parameters.
    % the current version uses only C and D. A and B will be used in the future
    A = tran_func(1:num_state, 1:num_state);
    B = tran_func(1:num_state, num_state+1:num_state+K);
    C = tran_func(num_state+1:num_state+N, 1:num_state);
    D = tran_func(num_state+1:num_state+N, num_state+1:num_state+K);

    % pre-process, determine the location
    [n_tran_prob, m_tran_prob] = size(tran_prob);
    if n_tran_prob == 3
      expense_flag = 1;     % expense flag == 0; BSC code. == 1, with TRAN_PROB.
    else
      expense_flag = 0;
    end;

    % varaibles 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.

    % pre-process, determine the location
    nn_code = num_std_sta * N * 1000;
    starter = x(5);
    trace = [];
    solut = [];
    store_code = [];

    % make the point of meaningful variabe to the storage variable x.
    for i = 1 : leng
      trace = [trace; x((i-1)*num_std_sta+6 : i*num_std_sta+5)'];
      solut = [solut; x((i+leng-1)*num_std_sta+6 : (i+leng) * num_std_sta + 5)'];
      store_code = [store_code; x(2*leng*num_std_sta+(i-1)*N+6 : 2*leng*num_std_sta+i*N+5)'];
    end;
    expense = [x(2*leng*num_std_sta+leng*N+6:(2*leng+1)*num_std_sta+leng*N+5)';
              x((2*leng+1)*num_std_sta+leng*N+6:(2*leng+2)*num_std_sta+leng*N+5)'];

    K2 = 2^K;

    trace_num = x(3) + 1;
    trace_flag = x(4);

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

    %clean-up spaces.
    trace_pre = rem(trace_num - 2 + leng, leng) + 1;
    trace(trace_num, :) = -ones(1, num_std_sta);
    solut(trace_num, :) = zeros(1, num_std_sta);
    store_code(trace_num,:) = u(1:N);

    % fill up trace and solut
    if (trace_flag == 0) & (trace_num == 1)
      pre_state = starter + 1;
    else
      pre_state = find(trace(trace_pre, :) >= 0);
    end;

    % loop in calculating the expense 
    for j = 1 : length(pre_state)
      jj = pre_state(j) - 1;           % index number - 1 is the state.
      cur_sta = de2bi(jj, num_state)';

      % find the expense region
      if expense_flag
        for num_N = 1 : N
          expen_work(num_N) = max(find(tran_prob(1,:) <= u(num_N)));
        end;
      end;

      % update the state and calculate the expense      
      for num_k = 1 : K2
        inp = de2bi(num_k - 1, K)';
        out = rem(C * cur_sta + D * inp,2);
        % X = A x + Bu to be added here.
        nex_sta = [inp; cur_sta(1 : num_state - K)];
        nex_sta_de = bi2de(nex_sta') + 1;
        if expense_flag
          % find the expense by the transfer probability
          expen_0 = find(out' <= 0.5);
          expen_1 = find(out' > 0.5);
          local_expense = sum([tran_prob(2,expen_work(expen_0)) 0])...
                      +sum([tran_prob(3,expen_work(expen_1)) 0]);
        else
          local_expense = sum(rem(u(1:N) + out', 2));
        end;

        if trace(trace_num, nex_sta_de) < 0
          % nothing happened yet. Fill all information.
          % previous state.
          trace(trace_num, nex_sta_de) = jj;
          % input information.
          solut(trace_num, nex_sta_de) = num_k - 1;
          % weighting/expense
          expense(2, nex_sta_de) = local_expense;
        else
          % something happended before. Take a comparison.
          replace_flag = 0;
          if expense_flag
            if (local_expense + expense(1, jj+1)) > (expense(2, nex_sta_de) + expense(1, trace(trace_num, nex_sta_de)+1))
              replace_flag = 1;
            end;
          else
            if (local_expense + expense(1, jj+1)) < (expense(2, nex_sta_de) + expense(1, trace(trace_num, nex_sta_de)+1))
              replace_flag = 1;
            end;
          end;
          if replace_flag
            % make a replacement other wise do nothing
            expense(2, nex_sta_de) = local_expense;
            solut(trace_num, nex_sta_de) = num_k - 1;
            trace(trace_num, nex_sta_de) = jj;
          end;
        end;
      end;
    end;

    % rearrange the first line of expense.
    tmp = expense(1,:);
    for j = 1 : num_std_sta
      if trace(trace_num, j) >= 0
        expense(1,j) = expense(2,j) + tmp(trace(trace_num, j)+1);
      end;
    end;

    % decision making.
    if trace_flag
      % strike out the unnecessary.
      for j_k = 1 : leng - 1
        j = rem(trace_num - j_k + leng, leng) + 1;
        j_pre = rem(trace_num - j_k - 1 + leng, leng) + 1;
        pre_state = find(trace(j_pre, :) >= 0);
        for num_k = 1 : length(pre_state)
          jj = pre_state(num_k) - 1;
          % in case of a state has not been continued, it will be strike out.
          if isempty(find(trace(j, :) == jj))
            trace(j_pre, pre_state(num_k)) = -1;
          end;
        end;
      end;

      % detect the converged result. eliminate a memory space.
      trace_eli = rem(trace_num, leng) + 1;
      trace_eli_non_z = find(trace(trace_eli,:) >= 0);
      if length(trace_eli_non_z) > 1
        % need to make a decision here.
        % take the one with the complete path has the best expense.
        if expense_flag
          right_exp = find(expense(1, :) == max(expense(1,:)));
        else
          right_exp = find(expense(1, :) == min(expense(1,:)));
        end;
        right_exp = right_exp(1) - 1;
        for j_k = 1 : leng - 1
          % trace from trace_num to trace_eli
          tmp = rem(trace_num - j_k + leng, leng) + 1;
          right_exp = trace(tmp, right_exp + 1);
        end;
        trace_eli_non_z = right_exp + 1;

        tmp = find(trace(trace_eli, :) ~= -1);
        tmp1 = find(tmp == trace_eli_non_z);
        tmp(tmp1) = [];
        for j_k = 1 : leng-1
          % clean the path which is derived from trace_non_z
          tmp1 = tmp;
          tmp = [];
          index_tmp = rem(trace_eli + j_k - 1, leng) + 1;
          for j_j = 1 : length(tmp1)
            tmp = [tmp find(trace(index_tmp, :) == tmp1(j_j)-1)];
          end;
          trace(index_tmp, tmp) = -ones(1, length(tmp));
        end;

        if expense_flag
          expense(1, tmp) = -ones(1, length(tmp)) * nn_code;
        else
           expense(1, tmp) = ones(1, length(tmp)) * nn_code;
        end;
      end;

      inp = de2bi(solut(trace_eli, trace_eli_non_z), K);
      cur_sta = de2bi(starter, num_state);
      out = rem(C * cur_sta' + D * inp', 2);

      if expense_flag
        % find the expense by the transfer probability
        for num_N = 1 : N
          expen_work(num_N) = max(find(tran_prob(1,:) <= store_code(trace_eli,num_N)));
        end;
        expen_0 = find(out' <= 0.5);
        expen_1 = find(out' > 0.5);
        local_expense = sum([tran_prob(2,expen_work(expen_0)) 0])...
                       +sum([tran_prob(3,expen_work(expen_1)) 0]);
      else
        local_expense = sum(rem(store_code(trace_eli,:) + out', 2));
      end;
      ind_exp = find(expense(1,:) ~= nn_code);
      starter = trace_eli_non_z - 1;
      output = [inp(:); local_expense];
    else %(if trace_flag)
      output = zeros(K+1, 1);
    end; %(if trace_flag)

    trace = trace';
    solut = solut';
    expense = expense';
    store_code = store_code';
    trace_num = rem(trace_num, leng);
    sys = [0; 0; trace_num; trace_flag; starter; trace(:); solut(:); store_code(:); expense(:); output];
elseif flag == 3 % output
    % if there is a trigger signal, output the new calculation
    % otherwise keep the old one.
    if u(length(u)) > 0.2
        K = code_param(2);
        len_x = length(x);
        sys = x(len_x - K: len_x);
    end;
elseif flag == 0
    N = code_param(1);
    K = code_param(2);
    M = code_param(3);
    num_state = code_param(4);
    num_std_sta = code_param(5);

    [n_tran_prob, m_tran_prob] = size(tran_prob);
    if n_tran_prob == 3
        expense_flag = 1;     % expense flag == 0; BSC code. == 1, with TRAN_PROB.
    else
        expense_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.
    nn_code = num_std_sta * N * 1000;
    if expense_flag 
        expense = [-ones(1, num_std_sta) * nn_code; zeros(1, num_std_sta)];
    else
        expense = [ones(1, num_std_sta) * nn_code; zeros(1, num_std_sta)];
    end;
    expense(1,1) = 0;

    % the column number is the state + 1.
    % the contents is the state it transfered from.
    % the starter is always single number.
    starter = 0;
    % a fixed length
    trace   = zeros(leng, num_std_sta);

    % 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.
    solut   = zeros(leng, num_std_sta);

    if plot_flag > 0
        x0 = -Inf;
    else
        x0 = 0;
    end;
    y = zeros(K+1, 1);
    store_code = zeros(leng*N, 1);
    expense = expense';
    % add output storage.
    x0 = [x0; 0; 0; 0; starter; trace(:); solut(:); store_code; expense(:); y];
%          |  |  |  |      |       |        |           |         |         \- output
%          \  \  \  \      \       \        \           \         \- weight, expense
%           \  \  \  \      \start  \-state  \-decode    \- code_storage
%            \  \  \  \-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
    sys = [];
end;