vitcore.c

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

/*  P.S: This is a data processing function called by viterbi.m     */

/*-----------------------------------------------------------------------------
 *       Original Designed by Wes Wang
 *       Jun Wu        The MAthworks, Inc.
 *       Mar-3rd, 1996          
 *
 *       Copyright (c) 1996 by The MathWorks, Inc.
 *       All Rights Reserved
 *       $Revision: 1.1 $  $Date: 1996/04/01 18:14:52 $
 *-----------------------------------------------------------------------------
 */
#include <stdio.h>
#include <math.h>
#include "mex.h"

#include "vitlib.c"
#define Inf     mexGetInf()
#define one 1
/*
 * the main function
 */
void mexFunction(int nlhs, Matrix *plhs[], int nrhs, Matrix *prhs[])
{
  int    i, j, ii, jj, j2, j_k, j_pre, l, indx_j, num_N, num_K, NaN;
  int    n_code, m_code, colFunc, rowFunc, n_tran_prob, m_tran_prob, n_codd, n_msg, len_expen;
  int    N, M, K, num_state, n_std_sta, PowPowM, K2; 
  int    leng, len_C, len_pre_state, len_aft_state, len_plot_flag;
  int    expen_flag, trace_num, trace_eli, trace_pre, tran_indx, nex_sta_de, cur_sta_num;
  int    starter, dec, loc_tmp, numnotnan, loc_exp1, loca_exp1, tmp, lenIndx0, lenIndx1;
  double max, min, loca_exp, loc_exp;

  int    *A, *B, *C, *D, *expense1, *solution1, *solu, *sol1, *TRAN_A, *TRAN_B;  
  int    *inp_pre, *cur_sta_pre, *expen_work, *expen_tmp1, *pre_state, *cur_sta, *inp;
  int    *nex_sta, *out, *expenOut, *aft_state, *tmpIwork, *IWork; 
  double *code, *tran_func, *tran_prob, *Tran_prob, *plot_flag;
  double *msg, *expen, *codd;
  double *expense, *expen_tmp, *sol, *solution, *tmpRwork, *RWork;

  /*% routine check.
   *if nargin < 2
   *    error('Not enough input variable for VITERBI.');
   *elseif nargin < 3
   *    leng = 0;
   *end;
   *
   *% the case of unlimited delay and memory.
   *[n_code, m_code] = size(code);
   *if (leng <= 0) | (n_code <= leng)
   *    leng = n_code;
   *end;
   *if leng <= 1
   *    leng = 2;
   *end;
   *
   *if nargin < 4
   *    expen_flag = 0;     % expense_flag == 0 for BSC code;  == 1, with TRAN_PROB.
   *else
   *    [N, K] = size(tran_prob);
   *    if N < 3
   *        expen_flag = 0;
   *    else
   *        expen_flag = 1;
   *        expen_work = zeros(1, N); % work space.
   *        tran_prob(2:3, :) = log10(tran_prob(2:3, :));
   *    end;
   *end;
   *
   *if nargin < 5
   *    plot_flag(1) = 0;
   *else
   *    if (plot_flag(1) <= 0)
   *        plot_flag(1) = 0;
   *    elseif floor(plot_flag(1)) ~= plot_flag(1)
   *        plot_flag(1) = 0;
   *    end;
   *end;
   */

  /* deal properly with all kinds of input. */
  /* get the input(code, tran_func) of function */
  n_code = mxGetM(prhs[0]);
  m_code = mxGetN(prhs[0]);
  code = mxGetPr(prhs[0]);

  rowFunc = mxGetM(prhs[1]);
  colFunc = mxGetN(prhs[1]);
  tran_func = mxGetPr(prhs[1]);
  for(i=0; i < rowFunc*colFunc; i++){
    if( mxGetPr(prhs[1])[i] < 0 )
      tran_func[i] = -Inf;
  }

  if( tran_func[rowFunc*colFunc-1] == -Inf ){
    /* this kind of tran_func has A, B, C, D */ 
    N = (int)tran_func[ rowFunc*(colFunc-1) ];
    K = (int)tran_func[ rowFunc*(colFunc-1)+1 ];
    M = (int)tran_func[ rowFunc*(colFunc-1)+2 ];
    len_C = M*N;
  }else{
    /* this kind of tran_func is a two columns matrix */ 
    N = (int)tran_func[rowFunc];
    K = (int)tran_func[rowFunc+1];
    M = (int)tran_func[1];
    len_C = 0;
  }

  if(nrhs<3){
    leng = 0;
    expen_flag = 0;
    NaN = -1;
    len_plot_flag = 0;
  } else {
    leng = (int)mxGetPr(prhs[2])[0];  
    if( nrhs < 4 ){
      expen_flag = 0;
      NaN = -1;
      len_plot_flag = 0;
    } else {
      /* deal properly with all kinds of 'tran_prob' */
      n_tran_prob = mxGetM(prhs[3]);
      m_tran_prob = mxGetN(prhs[3]);
      Tran_prob = mxGetPr(prhs[3]);
      if(n_tran_prob < 3){
        expen_flag = 0;
        NaN = -1;
      }else{
        expen_flag = 1;
        NaN = 1;
        tran_prob = (double *)mxCalloc(m_tran_prob*n_tran_prob, sizeof(double));
        for(i=0; i < m_tran_prob; i++){
          tran_prob[i*n_tran_prob] = Tran_prob[i*n_tran_prob];
          tran_prob[i*n_tran_prob+1] = log10(Tran_prob[i*n_tran_prob+1]);
          tran_prob[i*n_tran_prob+2] = log10(Tran_prob[i*n_tran_prob+2]);
        }
      }

      /* deal properly with all kinds of 'plot_flag' */
      if( nrhs < 5 ){
        len_plot_flag = 0;
      }else{
        len_plot_flag = mxGetM(prhs[4])*mxGetN(prhs[4]);
        plot_flag = mxGetPr(prhs[4]);
        if( plot_flag[0] <= 0 )
          len_plot_flag = 0;
        else if ( (int)plot_flag[0] != plot_flag[0] )
          len_plot_flag = 0;
      }
    }
  }

  if( leng <= 0 || n_code <= leng )
    leng = n_code;
  if( leng <= 1 )
    leng = 2;
    
  /*if isstr(tran_func)
   *    tran_func = sim2tran(tran_func);
   *end;
   *% initial parameters.
   *[A, B, C, D, N, K, M] = gen2abcd(tran_func);
   */
  num_state = M;
  n_std_sta = 1;
  for(i=0; i < num_state; 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( expen_flag == 1 ){
    RWork =(double *)mxCalloc( (leng+4)*PowPowM+2*n_std_sta, sizeof(double));
                /* the number of real working space
                 * expense -------- leng*PowPowM
                 * solution ------- PowPowM
                 * sol ------------ PowPowM
                 * expen_tmp ------ PowPowM
                 * tmpRwork ------- 2*n_std_sta+PowPowM
                 */
    if( len_C != 0 )
      IWork =(int *)mxCalloc((M+N)*(M+K)+leng*PowPowM+K*(K2+1)+(4+M)*n_std_sta+3*N+2*M, sizeof(int));
                    /* A -------------- M*M
                     * B -------------- M*K
                     * C -------------- N*M
                     * D -------------- N*K
                     * solu ----------- leng*PowPowM
                     * inp_pre -------- K*K2
                     * cur_sta_pre ---- M*n_std_sta
                     * 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
      IWork = (int *)mxCalloc((rowFunc-2)*(M+N+2)+leng*PowPowM+K*(K2+1)+(4+M)*n_std_sta+3*N+2*M, sizeof(int));
                    /* TRAN_A --------- rowFunc-2
                     * TRAN_B --------- rowFunc-2
                     * A -------------- (rowFunc-2)*M
                     * B -------------- (rowFunc-2)*N
                     * same as above from *solu
                     */
  }else{
    /* no real working space */
    if( len_C != 0 )
      IWork =(int *)mxCalloc((M+N)*(M+K)+(2*leng+4)*PowPowM+K*(K2+1)+(6+M)*n_std_sta+2*N+2*M, sizeof(int));
                    /* A -------------- M*M
                     * B -------------- M*K
                     * C -------------- N*M
                     * D -------------- N*K
                     * expense1 ------- leng*PowPowM
                     * solution1 ------ PowPowM
                     * solu ----------- leng*PowPowM
                     * 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
      IWork = (int *)mxCalloc((rowFunc-2)*(M+N+2)+(2*leng+4)*PowPowM+K*(K2+1)+(6+M)*n_std_sta+2*N+2*M, sizeof(int));
                    /* TRAN_A --------- rowFunc-2
                     * TRAN_B --------- rowFunc-2
                     * A -------------- (rowFunc-2)*M
                     * B -------------- (rowFunc-2)*N
                     * same as above from *expense1
                     */
  }

  if( len_C != 0 ){
    A = IWork;
    B = A + M*M;
    C = B + M*K;
    D = C + M*N;
        
    /* Get the input Matrix A, B, C, D */
    for( i=0; i < M+N; i++ ){
      for( j=0; j < M+K; j++ ){
        if( i<M   && j<M )
          A[i+j*M] = (int)tran_func[i+j*(M+N)];
        if( i>M-1 && j<M )
          C[i+j*N-M] = (int)tran_func[i+j*(M+N)];
        if( i<M   && j>M-1 )
          B[i+j*M-M*M] = (int)tran_func[i+j*(M+N)];
        if( i>M-1 && j>M-1 )
          D[i+j*N-M*(N+1)] = (int)tran_func[i+j*(M+N)];
      }
    }
  }else{
    /* Assignment */
    TRAN_A = IWork;
    TRAN_B = TRAN_A + rowFunc-2;
    A = TRAN_B + rowFunc-2;
    B = A + M*(rowFunc-2);

    /* Get the input Matrix A, B */
    for(i=0; i < rowFunc-2; i++){
      TRAN_A[i] = (int)tran_func[i+2];
      TRAN_B[i] = (int)tran_func[rowFunc+i+2];
    }
    de2bi(TRAN_A, M, rowFunc-2, A);
    de2bi(TRAN_B, N, rowFunc-2, B);
  }
  /* end of Input Segment */

  if( expen_flag != 0 ){
    expense     = RWork;
    solution    = expense + leng*PowPowM;
    sol         = solution + PowPowM;
    expen_tmp   = sol + PowPowM;
    tmpRwork    = expen_tmp + PowPowM;

    if( len_C != 0 )
      solu = D + N*K;
    else
      solu = B + N*(rowFunc-2);

    inp_pre = solu + leng*PowPowM;        
    cur_sta_pre = inp_pre + K2*K;
    expen_work = cur_sta_pre + M*n_std_sta;
    pre_state = expen_work + N;
    cur_sta = pre_state + n_std_sta;
    inp = cur_sta + M;
    nex_sta = inp + K;
    out = nex_sta + M;
    expenOut = out + N;
    aft_state = expenOut + N;
    tmpIwork = aft_state + n_std_sta;
  }else{ /* tran_prob is not 3-row matrix */
    if( len_C != 0 )
      expense1 = D + N*K;
    else
      expense1 = B + N*(rowFunc-2);

    solu = expense1 + leng*PowPowM;    
    solution1 = solu + leng*PowPowM;
    inp_pre = solution1 + PowPowM;
    cur_sta_pre = inp_pre+ K2*K;
    pre_state =  cur_sta_pre + M*n_std_sta;
    cur_sta = pre_state + n_std_sta;
    inp = cur_sta + M;
    nex_sta = inp + K;
    out = nex_sta + M;
    expenOut = out + N;
    aft_state = expenOut + N;
    sol1 = aft_state + n_std_sta;
    expen_tmp1 = sol1 + PowPowM;
    tmpIwork = expen_tmp1 + PowPowM;
  }  

  /*if m_code ~= N
   *    if n_code == 1
   *        code = code';
   *        [n_code, m_code] = size(code);
   *    else
   *        error('CODE length does not match the TRAN_FUNC dimension in VITERBI.')
   *    end;
   *end;
   */

  if( m_code != N ){
    if(n_code == 1 ){
      n_code = m_code;
      m_code = 1;
	} else {
      mexErrMsgTxt("Code length does not match the TRAN_FUNC dimension in VITERBI.");
    }
  }

  /*% state size
   *if isempty(C)
   *    states = zeros(tran_func(2,1), 1);
   *else
   *    states = zeros(length(A), 1);
   *end;
   *num_state = length(states);
   */
  /* num_state has been set by the different TRAN_FUNC */

  /*% The STD state number reachiable is
   *n_std_sta = 2^num_state;
   *
   *PowPowM = n_std_sta^2;
   *
   *% at the very begining, the current state number is 1, only the complete zeros.
   *cur_sta_num = 1;
   *
   *% the variable expense is divided into n_std_sta section of size n_std_sta
   *% length vector for each row. expense(k, (i-1)*n_std_sta+j) is the expense of
   *% state j transfering to be state i
   *expense = ones(leng, PowPowM) * NaN;
   *expense(leng, [1:n_std_sta]) = zeros(1, n_std_sta);
   *solu = zeros(leng, PowPowM);
   *solution = expense(1,:);
   *
   *K2 = 2^K;
   */
  cur_sta_num = 1;