sviterbi.c

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

/*============================================================================
 *  Syntax: [sys, x0, str, ts] =
 *      sviterbi(t,x,u,flag,tran_func,leng,tran_prob,plot_flag,V1,V2,V3,V4);
 *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.
 *
 * P.S. The most of comments are written in Matlab file.
 *=============================================================================
 * Original designed by Wes Wang,
 * Jun Wu,  The Mathworks, Inc.
 * Feb-07, 1996
 *
 * Copyright (c) 1996 by The MAthWorks, Inc.
 * All Rights Reserved
 * $Revision: 1.1 $  $Date: 1996/04/01 19:06:23 $
 *=============================================================================
 */
#define S_FUNCTION_NAME sviterbi

/* M-files sviplot1.m sviplot2.m sviplot3.m and sviplot4.m
 * are needed if you want plot the figures.
 */
#ifdef MATLAB_MEX_FILE
#include "mex.h"
#endif

#include <math.h>

/* need to include simstruc.h for the definition of the SimStruct and
 * its associated macro definitions.
 */ 
#include "simstruc.h"

#define NUM_ARGS            8           /* eight additional input arguments */
#define TRAN_FUNC       ssGetArg(S,0)   /* the length of  output vector */
#define LENG            ssGetArg(S,1)   /* the memory length */
#define TRAN_PROB       ssGetArg(S,2)   /* the transfer probability */
#define PLOT_FLAG       ssGetArg(S,3)   /* the flag of plot */
#define V1              ssGetArg(S,4)
#define V2              ssGetArg(S,5)  
#define V3              ssGetArg(S,6)  
#define V4              ssGetArg(S,7)  

#define Prim 2
void de2bi(pde, dim, pow_dim, pbi)
    int *pde, dim, pow_dim, *pbi;
{
  int     i,j, tmp;

  if( pde[0] < 0 ){
    /* the first part is for converting the decimal numbers(from 0 to pow_dim)
     * to binary.
     */  
    for(i=0; i < pow_dim; i++){
      tmp = i;
      for(j=0; j < dim; j++){
        pbi[i+j*pow_dim] = tmp % Prim;
        if(j < dim-1)
          tmp = (int)(tmp/Prim);
      }
    }
  }else{
    /* the second part is for converting the decimal numbers(setting by user)
     * to binary.
     */  
    for(i=0; i < pow_dim; i++){
      tmp = pde[i];
      for(j=0; j < dim; j++){
        pbi[i+j*pow_dim] = tmp % Prim;
        if(j < dim-1)
          tmp = (int)(tmp/Prim);
      }
    }
  }
}
static void bi2de(pbi, pow_dim, dim, pde)
     int *pbi, pow_dim, dim, *pde;
{
  int     i, j, k;
  
  for(i=0; i<pow_dim; i++)
    pde[i] = 0;
  
  for (i=0; i < pow_dim; i++){
    for (j=0; j < dim; j++){
      if (pbi[i+j*pow_dim] != 0){
        if(j > 0){
          for (k=0; k < j; k++)
            pbi[i+j*pow_dim] = pbi[i+j*pow_dim]*Prim;
        }
      }
      pde[i] = pde[i] + (int)pbi[i+j*pow_dim];
    }
  }
}
/*
 * See also vitshort.c and vitshort.m for the two functions following.
 */
static void shortdbl(expense, sol, n_std_sta, Rwork, Iwork)
     double *expense, *sol, *Rwork; 
     int *Iwork, n_std_sta;
{        
  int     i, j, k, PowPowM, aft_indx, len_indx, len, indx;
  int     *wei_indx, *sub_indx;
  double  max;
  double  *wei_mat_exp, *wei_sum, *sol_tmp;
  double  wei_mat_sol;    
  
  wei_mat_exp = Rwork;
  wei_sum     = Rwork + n_std_sta;
  sol_tmp     = Rwork + 2*n_std_sta;
  wei_indx    = Iwork;
  sub_indx    = Iwork + n_std_sta;
  
  PowPowM = n_std_sta * n_std_sta;
  for(i=0; i < PowPowM; i++)
    sol_tmp[i] = sol[i];
  for(i=1; i <= n_std_sta; i++){
    aft_indx = i * n_std_sta - n_std_sta;
    for(j=1; j <= n_std_sta; j++){
      for(k=0; k < n_std_sta; k++)
        wei_mat_exp[k] = expense[k * n_std_sta + j - 1];
      len_indx = 0;
      for(k=0; k < n_std_sta; k++){
        wei_mat_sol = sol_tmp[aft_indx + k];
        if ( wei_mat_exp[k] > 0 || wei_mat_sol > 0 ) {
          wei_sum[k] = 1;
	    }else{
          wei_sum[k] = wei_mat_exp[k] + wei_mat_sol;
          wei_indx[len_indx] = k;
          len_indx++;
	    }
      }
      
      if (len_indx > 0) {
        len = 0;
        max = wei_sum[wei_indx[0]];
        sub_indx[0] = wei_indx[0];                    
        k = 1;
        while (k < len_indx) {
          if ( max < wei_sum[wei_indx[k]] ) {
            max = wei_sum[wei_indx[k]];
            sub_indx[0] = wei_indx[k];
          }
          k++;
        }
        for(k=0; k < len_indx; k++){
          if (wei_sum[wei_indx[k]] == max ){
            sub_indx[len] = wei_indx[k];
            len++;
          }
        }
        indx = sub_indx[0];                        
        if (len > 1){
          max = wei_mat_exp[sub_indx[0]];
          k = 1;
          while(k < len){
            if( max < wei_mat_exp[sub_indx[k]] ) {
              max = wei_mat_exp[sub_indx[k]];
              indx = sub_indx[k];
            }
            k++;
          }
        }
        sol[aft_indx + j - 1] = wei_sum[indx];
      }else{
        sol[aft_indx + j - 1] = 1;
      }
    }
  }
}
static void shortint(expense, sol, n_std_sta, Iwork)
     int     *expense, *sol, *Iwork;
     int     n_std_sta;
{
  int     i, j, k, PowPowM, aft_indx, len_indx, len, indx;
  int     min;
  int     wei_mat_sol;    
  int     *wei_mat_exp, *wei_sum, *sol_tmp, *wei_indx, *sub_indx;
  
  wei_mat_exp = Iwork; 
  wei_sum     = Iwork + n_std_sta;
  wei_indx    = Iwork + 2*n_std_sta;
  sub_indx    = Iwork + 3*n_std_sta;
  sol_tmp     = Iwork + 4*n_std_sta;
  
  PowPowM = n_std_sta * n_std_sta;
  for(i=0; i < PowPowM; i++)
    sol_tmp[i] = sol[i];
  for(i=1; i <= n_std_sta; i++){
    aft_indx = i * n_std_sta - n_std_sta;
    for(j=1; j <= n_std_sta; j++){
      for(k=0; k < n_std_sta; k++)
	    wei_mat_exp[k] =expense[k * n_std_sta + j - 1];
      len_indx = 0;
      for(k=0; k < n_std_sta; k++){                
	    wei_mat_sol = sol_tmp[aft_indx + k];
	    if ( wei_mat_exp[k] < 0 || wei_mat_sol < 0 ) {
	       wei_sum[k] = -1;
	    }else{
	       wei_sum[k] = wei_mat_exp[k] + wei_mat_sol;
	       wei_indx[len_indx] = k;
	       len_indx++;
	    }
      }
      
      if (len_indx > 0) {
	     len = 0;
	     min = wei_sum[wei_indx[0]];
	     sub_indx[0] = wei_indx[0];                    
	     k = 1;
	     while (k < len_indx) {
	       if ( min > wei_sum[wei_indx[k]] ) {
	         min = wei_sum[wei_indx[k]];
	         sub_indx[0] = wei_indx[k];
	       }
	       k++;
	     }
	     for(k=0; k < len_indx; k++){
	       if (wei_sum[wei_indx[k]] == min ){
	         sub_indx[len] = wei_indx[k];
	         len++;
	       }
	     }
	     indx = sub_indx[0];                        
	     if (len > 1){
	       min = wei_mat_exp[sub_indx[0]];
	       k = 1;
	       while(k < len){
	         if( min > wei_mat_exp[sub_indx[k]] ) {
	           min = wei_mat_exp[sub_indx[k]];
	           indx = sub_indx[k];
	         }
	         k++;
	       }
	     }
	     sol[aft_indx + j - 1] = wei_sum[indx];
      }else{
	     sol[aft_indx + j - 1] = -1;
      }
    }
  }
}
static void mdlInitializeSizes(S)
    SimStruct *S;
{
  int     i;
  int     N, M, K, K2, colFunc, rowFunc, n_tran_prob, m_tran_prob;
  int     num_state, n_std_sta, PowPowM;
  int     leng, plot_flag, expen_flag, len_C;

  leng = (int)mxGetPr(LENG)[0];
  plot_flag = (int)mxGetPr(PLOT_FLAG)[0];
  colFunc = mxGetN(TRAN_FUNC);
  rowFunc = mxGetM(TRAN_FUNC);
  n_tran_prob = mxGetM(TRAN_PROB);
  m_tran_prob = mxGetN(TRAN_PROB);
  
  if( mxGetPr(TRAN_FUNC)[rowFunc*colFunc-1] < 0 ){
    N = (int)mxGetPr(TRAN_FUNC)[ rowFunc*(colFunc-1) ];
    K = (int)mxGetPr(TRAN_FUNC)[ rowFunc*(colFunc-1)+1 ];
    M = (int)mxGetPr(TRAN_FUNC)[ rowFunc*(colFunc-1)+2 ];
    len_C = M*N;

  }else{
    N = (int)mxGetPr(TRAN_FUNC)[rowFunc];
    K = (int)mxGetPr(TRAN_FUNC)[rowFunc+1];
    M = (int)mxGetPr(TRAN_FUNC)[1];
    len_C = 0;
  }

  num_state = M;
  n_std_sta = 1;
  for(i=0; i < M; i++)
    n_std_sta = n_std_sta * 2;
  PowPowM = n_std_sta * n_std_sta;
  K2 = 1;
  for(i=0; i < K; i++)
    K2 = K2 * 2;

  if ( n_tran_prob == 3 )
    expen_flag = 1;
  else
    expen_flag = 0;
	
  ssSetNumContStates(     S, 0);          /* number of continuous states */
  ssSetNumDiscStates(     S, 8+2*leng*PowPowM+leng*N+K); /* number of discrete states */
  ssSetNumOutputs(        S, K+1);        /* number of outputs */
  ssSetNumInputs(         S, N+1);        /* number of inputs */
  ssSetDirectFeedThrough( S, 0);          /* direct feedthrough flag */
  ssSetNumSampleTimes(    S, 1);          /* number of sample times */
  ssSetNumInputArgs(      S, NUM_ARGS);   /* number of input arguments */
  if( expen_flag == 1 ){
    ssSetNumRWork(      S, n_tran_prob*m_tran_prob+(leng+3)*PowPowM+K+1+2*n_std_sta);/* number of real working space */
                /* tran_prob_tmp -------- n_tran_prob*m_tran_prob
                 * expense -------- leng*PowPowM
                 * Y -------------- K+1
                 * sol ------------ PowPowM
                 * expen_tmp ------ PowPowM
                 * tmpRwork ------- 2*n_std_sta+PowPowM
                 */
    if( mxGetPr(TRAN_FUNC)[rowFunc*colFunc-1] < 0 ){
      ssSetNumIWork( S,  (M+N)*(M+K)+leng*(PowPowM+N)+K*(K2+1)+(4+M)*n_std_sta+3*N+2*M);
                    /* A -------------- M*M
                     * B -------------- M*K
                     * C -------------- N*M
                     * D -------------- N*K
                     * solu ----------- leng*PowPowM
                     * code ----------- leng*N
                     * inp_pre -------- K*2^K
                     * cur_sta_pre ---- M*2^M
                     * expen_work ----- N
                     * pre_state ------	n_std_sta
                     * cur_sta -------- M
                     * inp ------------ K
                     * nex_sta -------- M
                     * out ------------ N
                     * expenOut ------- N
                     * aft_state ------	n_std_sta
                     * tmpIwork ------- 2*n_std_sta
                     */
    }else{
      ssSetNumIWork( S, (rowFunc-2)*(M+N+2)+leng*(N+PowPowM)+K*(K2+1)+(4+M)*n_std_sta+3*N+2*M);
                    /* TRAN_A --------- rowFunc-2
                     * TRAN_B --------- rowFunc-2
                     * A -------------- (rowFunc-2)*M
                     * B -------------- (rowFunc-2)*N
                     * same as above from *solu
                     */
    }
  }else{
    ssSetNumRWork(          S, 0);  /* no real working space */

    if( mxGetPr(TRAN_FUNC)[rowFunc*colFunc-1] < 0 ){
      ssSetNumIWork( S,  (M+N)*(M+K)+(2*leng+3)*PowPowM+K*(K2+2)+(6+M)*n_std_sta+(leng+2)*N+2*M+1);
                    /* A -------------- M*M
                     * B -------------- M*K
                     * C -------------- N*M
                     * D -------------- N*K
                     * expense1 ------- leng*PowPowM
                     * solu ----------- leng*PowPowM
                     * code ----------- leng*N
                     * Y1 ------------- K+1
                     * inp_pre -------- K*2^K
                     * cur_sta_pre ---- M*2^M
                     * pre_state ------	n_std_sta
                     * cur_sta -------- M
                     * inp ------------ K
                     * nex_sta -------- M
                     * out ------------ N
                     * expenOut ------- N
                     * aft_state ------	n_std_sta
                     * sol1 ----------- PowPowM
                     * expen_tmp1 ----- PowPowM
                     * tmpIwork ------- 4*n_std_sta+PowPowM
                     */
    }else{
      ssSetNumIWork( S, (rowFunc-2)*(M+N+2)+(2*leng+3)*PowPowM+K*(K2+2)+(6+M)*n_std_sta+(leng+2)*N+2*M+1);
                    /* TRAN_A --------- rowFunc-2
                     * TRAN_B --------- rowFunc-2
                     * A -------------- (rowFunc-2)*M
                     * B -------------- (rowFunc-2)*N
                     * same as above from *expense1
                     */
    }
  }
  ssSetNumPWork(          S, 0);
}

/*
 * mdlInitializeConditions - initialize the states
 * Initialize the states, Integers and real-numbers
 */
static void mdlInitializeConditions(x0, S)
    double *x0;
    SimStruct *S;
{
    /*  [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.