sviterba.c

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

/*============================================================================
 *  Syntax: [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.
 *
 *=============================================================================
 * 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:12 $
 *=============================================================================
 */
#define S_FUNCTION_NAME sviterba

#ifdef MATLAB_MEX_FILE
#include "mex.h"
#include <stdio.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            4           /* four 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 Prim 2
static void de2bi(pde, dim, pow_dim, pbi)
     int *pde, dim, pow_dim, *pbi;
{
  int     i,j, tmp;

  if( pde[0] < 0 ){
    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{
    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];
    }
  }
}

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;
  
  leng = (int)mxGetPr(LENG)[0];
  plot_flag = (int)mxGetPr(PLOT_FLAG)[0];
  rowFunc = mxGetM(TRAN_FUNC);
  colFunc = mxGetN(TRAN_FUNC);
  n_tran_prob = mxGetM(TRAN_PROB);
  m_tran_prob = mxGetN(TRAN_PROB);
  
  if ( n_tran_prob == 3 )
    expen_flag = 1;
  else
    expen_flag = 0;
  
  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];
  }else{
    M = (int)mxGetPr(TRAN_FUNC)[1];
    N = (int)mxGetPr(TRAN_FUNC)[rowFunc];
    K = (int)mxGetPr(TRAN_FUNC)[rowFunc+1];
  }
  
  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;
  
  ssSetNumContStates(     S, 0);          /* number of continuous states */
  ssSetNumDiscStates(     S, 6+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
     * handles -------- 5+2*n_std_sta
     * 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);/* number of real working space */
    /* handles -------- 5+2*n_std_sta */
    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)