sviterbi.c

通信系统的matlab仿真程序

          mxGetPr(prhs4[0])[4] = dec;
          mxGetPr(prhs4[0])[5] = plot_flag;
          mxGetPr(prhs4[0])[6] = x[0];

          mexCallMATLAB(nlhs,plhs4,nrhs,prhs4,"sviplot4");

          mxFreeMatrix(plhs4[0]);
        }    
#endif

        /*    starter = dec;
         *    output = [inp(:); loca_exp];
         *  else %(if trace_flag)
         *    output = zeros(K+1, 1);
         *    starter = 0;
         *  end; %(if trace_flag)
         *
         *  trace_num = rem(trace_num, leng);
         *  sys = [x(1); fig_position; trace_num; trace_flag; starter; expense(:); solu(:); code(:); output(:)];
         */
        starter = dec;
        for(i=0; i < K; i++)
          Y[i] = (double)inp[i];
        Y[K] = loca_exp;    
      }else{  /* if (trace_flag != 0 ) */
        for(i=0; i < K+1; i++)
          Y[i] = 0;
        starter = 0;
      }  /* the end of "if (trace_flag != 0 )" */
      
      trace_num = trace_num % leng;
      
      x[1] = (double)fig_position;
      x[2] = (double)trace_num;
      x[3] = (double)trace_flag;
      x[4] = (double)starter;
      x[5] = (double)plot_flag_test;
      x[6] = (double)initial_flag;
      for(i=0; i < leng*PowPowM; i++)
    	x[i+7] = expense[i];
      for(i=0; i < leng*PowPowM; i++)
    	x[i+7+leng*PowPowM] = (double)solu[i];
      for(i=0; i < leng*N; i++)
    	x[i+7+2*leng*PowPowM] = (double)code[i];
      for(i=0; i < K+1; i++)
    	x[i+7+2*leng*PowPowM+leng*N] = Y[i];
    /*  the end of "if (expen_flag != 0 ) */
    }else{ /* tran_prob is not 3-row matrix */
      if( mxGetPr(TRAN_FUNC)[rowFunc*colFunc-1] < 0 )
        expense1 = D + N*K;
      else
        expense1 = B + N*(rowFunc-2);

      solu = expense1 + leng*PowPowM;    
      code = solu + leng*PowPowM;
      Y1 = code + leng*N;              /* size of *Y1 is (K+1) */
      inp_pre = Y1 + (K+1);            /* allocate K*2^K for *inp_pre */
      cur_sta_pre = inp_pre + K2*K;   /*  M*2^M for *cur_sta_pre. */
      pre_state = cur_sta_pre + M*n_std_sta;  /* allocate n_std_sta for *pre_state */
      cur_sta = pre_state + n_std_sta;/* allocate M for *cur_sta */
      inp = cur_sta + M;              /* allocate K for *inp */
      nex_sta = inp + K;              /* allocate M for *nex_sta */
      out = nex_sta + M;              /* allocate N for *out */
      expenOut = out + N;             /* allocate N for *expenOut */
      aft_state = expenOut + N;       /* allocate n_std_sta for *aft_state */
      sol1 = aft_state + n_std_sta;    /* allocate PowPowM for *sol */
      expen_tmp1 = sol1 + PowPowM;
      tmpIwork = expen_tmp1 + PowPowM;

      if ( x[6] != 12345 ){
        for(i=0; i < leng*PowPowM; i++){
          expense1[i] = NaN;
          solu[i] = 0;
        }        
        for(i=0; i < n_std_sta; i++)
          expense1[leng+i*leng-1] = 0;
        starter = 0;
        x[0] = 0;
        if(plot_flag > 0)
          x[5] = 1;
        else
          x[5] = 0;
        x[6] = 12345;
        for(i=0; i < leng*N; i++)
          code[i] = 0;
        for(i=0; i < K+1; i++)
          Y1[i] = 0;
        fig_position = 0;
        trace_num = 0;
        trace_flag = 0;
        
        x[1] = (double)fig_position;
        x[2] = (double)trace_num;
        x[3] = (double)trace_flag;
        x[4] = (double)starter;
      }

      inp_pre[0] = -1;
      cur_sta_pre[0] = -1;
      de2bi(cur_sta_pre, M, n_std_sta, cur_sta_pre);
      de2bi(inp_pre, K, K2, inp_pre);
      
      starter = (int)x[4];
      plot_flag_test = (int)x[5];
      initial_flag = (int)x[6];
      loc_tmp = 8;
#ifdef MATLAB_MEX_FILE
      if(plot_flag_test > 0){
        int nlhs, nrhs;
        Matrix *plhs1[1], *prhs1[1];

        nlhs = 1;
        nrhs = 1;
        prhs1[0] = V1;
        mxGetPr(prhs1[0])[0] = x[0];
        mxGetPr(prhs1[0])[1] = (double)n_std_sta;
        mxGetPr(prhs1[0])[2] = (double)num_state;
        mxGetPr(prhs1[0])[3] = (double)plot_flag;
        mxGetPr(prhs1[0])[4] = (double)initial_flag;
        
        mexCallMATLAB(nlhs, plhs1, nrhs, prhs1, "sviplot1");

        x[0] = mxGetPr(plhs1[0])[0];
        plot_flag_test = (int)mxGetPr(plhs1[0])[1];
        initial_flag = (int)mxGetPr(plhs1[0])[2];

        mxFreeMatrix(plhs1[0]);
      }
#endif

      for(i=0; i < leng*PowPowM; i++)
        expense1[i] = (int)x[loc_tmp-1+i];
      for(i=0; i < leng*PowPowM; i++)
        solu[i] = (int)x[loc_tmp+leng*PowPowM-1+i];
      for(i=0; i < leng*N; i++)
        code[i] = (int)x[loc_tmp+2*leng*PowPowM-1+i];

#ifdef MATLAB_MEX_FILE
      fig_position = x[1] + 1;
      if( x[0] > 0 && ((fig_position-leng)%(plot_flag - leng) == 0) && fig_position >= plot_flag && plot_flag_test != 0 ){
        int nrhs, nlhs;
        Matrix *plhs2[1], *prhs2[1];

        nlhs = 1;
        nrhs = 1;
        prhs2[0] = V1;
        mxGetPr(prhs2[0])[0] = (double)fig_position;
        mxGetPr(prhs2[0])[1] = (double)leng;
        mxGetPr(prhs2[0])[2] = (double)plot_flag;
        mxGetPr(prhs2[0])[3] = (double)n_std_sta;
        mxGetPr(prhs2[0])[4] = (double)x[0];

        mexCallMATLAB(nlhs,plhs2,nrhs,prhs2,"sviplot2");

        mxFreeMatrix(plhs2[0]);
      }
#endif

      trace_num = (int)x[2] + 1;
      trace_flag = (int)x[3];
      
      for(i=0; i < N; i++)
        code[trace_num-1+i*leng] = (int)u[i];
      
      if(trace_flag == 0 && trace_num == leng)
        trace_flag = 1;
      
      trace_pre = (trace_num - 2 + leng) % leng + 1;
      
      len_pre_state = 0;
      if( trace_flag == 0 && trace_num == 1 ){
        pre_state[0] = starter + 1;
        len_pre_state = 1;
      }else{
        for(j2=0; j2 < n_std_sta; j2++){
          numnotnan = 0;
          for(i=0; i < n_std_sta; i++){
            if( expense1[trace_pre-1 + i*leng+j2*leng*n_std_sta] >= 0 )
              numnotnan ++;
          }
          if(numnotnan != 0){
            pre_state[len_pre_state] = j2 + 1;
            len_pre_state++;
          }
        }
      }
      
      for(i=0; i < PowPowM; i++)
        expense1[trace_num-1+i*leng] = NaN;
      
      for(j=0; j < len_pre_state; j++){
        jj = pre_state[j] - 1;
        for(i=0; i < M; i++)
          cur_sta[i] = cur_sta_pre[jj + i*n_std_sta];
        indx_j = jj * n_std_sta;
        for(num_K=0; num_K < K2; num_K++){
          for(i=0; i < K; i++)
            inp[i] = inp_pre[num_K+i*K2];
          if( len_C == 0 ){
            tran_indx = pre_state[j] + num_K*n_std_sta;
            for(i=0; i < M; i++)
              nex_sta[i] = A[tran_indx-1+i*(rowFunc-2)];
            for(i=0; i < N; i++)
              out[i] = B[tran_indx-1+i*(rowFunc-2)];
          }else{
            for(i=0; i < N; i++){
              out[i] = 0;
              for(l=0; l < M; l++)
                out[i] = out[i] + C[i+l*N]*cur_sta[l];
              for(l=0; l < K; l++)    
                out[i] = out[i] + D[i+l*N]*inp[l];
              out[i] = out[i] % 2;
            }    
            for(i=0; i < M; i++){
              nex_sta[i] = 0;
              for(l=0; l < M; l++)
                nex_sta[i] = nex_sta[i] + A[i+l*M]*cur_sta[l];
              for(l=0; l < K; l++)
                nex_sta[i] = nex_sta[i] + B[i+l*M]*inp[l];
              nex_sta[i] = nex_sta[i] % 2;
            }
          }
          bi2de(nex_sta, 1, M, &nex_sta_de);
          nex_sta_de = nex_sta_de + 1;
          loca_exp1 = 0;
          for(i=0; i < N; i++)
            loca_exp1 = loca_exp1 + (code[trace_num-1+leng*i] + out[i]) % 2;
          tmp = (nex_sta_de - 1) * n_std_sta + pre_state[j] - 1;
          if( expense1[trace_num - 1 + tmp*leng] < 0 || expense1[trace_num - 1 + tmp*leng] > loca_exp1 ){
            expense1[trace_num - 1 + tmp*leng] = loca_exp1;
            solu[trace_num - 1 + leng*(nex_sta_de+indx_j-1)] = num_K + 1;
          }else if( expense1[trace_num - 1 + tmp*leng] > loca_exp1 ){
            expense1[trace_num - 1 + tmp*leng] = loca_exp1;
            solu[trace_num - 1 + leng*(nex_sta_de+indx_j-1)] = num_K + 1;
          }
        }
      }
      
      len_aft_state = 0;
      for(j2=0; j2 < n_std_sta; j2++){
        numnotnan = 0;
        for(i=0; i < n_std_sta; i++){
          if( expense1[trace_num-1+i*leng+j2*leng*n_std_sta] >= 0 )
            numnotnan ++;
        }
        if(numnotnan != 0){
          aft_state[len_aft_state] = j2 + 1;
          len_aft_state++;
        }
      }

#ifdef MATLAB_MEX_FILE
      /*  %%%%% begin plot related %%%%% */
      if (x[0] > 0 && plot_flag_test != 0 ){
        int nlhs, nrhs, aaa;
        Matrix *plhs3[1], *prhs3[4];

        nlhs = 1;
        nrhs = 4;
        prhs3[0] = V1;
        prhs3[1] = V2;
        prhs3[2] = V3;
        prhs3[3] = V4;
        mxGetPr(prhs3[0])[0] = (double)M;
        mxGetPr(prhs3[0])[1] = (double)trace_num;
        mxGetPr(prhs3[0])[2] = (double)expen_flag;
        mxGetPr(prhs3[0])[3] = (double)fig_position;
        mxGetPr(prhs3[0])[4] = (double)trace_flag;
        mxGetPr(prhs3[0])[5] = (double)leng;
        mxGetPr(prhs3[0])[6] = (double)x[0];
        for(i=0; i < len_pre_state; i++)
          mxGetPr(prhs3[1])[i] = (double)pre_state[i];
        for(i=0; i < leng*PowPowM; i++){
          if( expense1[i] == NaN )
            mxGetPr(prhs3[2])[i] = mexGetNaN();
          else
            mxGetPr(prhs3[2])[i] = (double)expense1[i];
        } 
        for(i=0; i < len_aft_state; i++)
          mxGetPr(prhs3[3])[i] = (double)aft_state[i];

        mexCallMATLAB(nlhs, plhs3,nrhs,prhs3,"sviplot3");
	
        mxFreeMatrix(plhs3[0]);
      }
#endif

      if( trace_flag != 0 ){
        trace_eli = (trace_num % leng) + 1;

        for( i=0; i < PowPowM; i++)
          sol1[i] = expense1[trace_num-1 + i*leng];
        for( j_k=1; j_k <= leng-2; j_k++){
          j_pre =(trace_num - j_k -1 + leng) % leng + 1;
          for( i=0; i < PowPowM; i++)
            expen_tmp1[i] = expense1[j_pre-1 + i*leng];
          shortint(expen_tmp1, sol1, n_std_sta, tmpIwork);
        }
        for(j=0; j < n_std_sta; j++){
          for(i=0; i < n_std_sta; i++){
            if( expense1[trace_eli-1+(starter+i*n_std_sta)*leng] < 0 )
              sol1[i+j*n_std_sta] = -1;
            else
              sol1[i+j*n_std_sta] = sol1[i+j*n_std_sta] + expense1[trace_eli-1+(starter+i*n_std_sta)*leng];
          }
        }

        for(i=0; i < PowPowM; i++){
          if( sol1[i] >= 0 ){
            loc_exp1 = sol1[i];
            i = PowPowM;
          }
        }
        for(i=0; i < PowPowM; i++){
          if( sol1[i] >= 0 && loc_exp1 > sol1[i])
            loc_exp1 = sol1[i];
        }
        for(i=0; i < PowPowM; i++){
          if( sol1[i] == loc_exp1 ){
            dec = i;
            i = PowPowM;
          }
        }
        dec = dec % n_std_sta;
        num_K = solu[trace_eli-1+leng*(starter*n_std_sta+dec)];
        for(i=0; i < K; i++)
          inp[i] = inp_pre[num_K-1+i*K2];
	
        if( len_C == 0 ){
          tran_indx = starter + 1 + (num_K-1)*n_std_sta;
          for(i=0; i < N; i++)
            out[i] = B[tran_indx-1+i*(rowFunc-2)];
        }else{
          for(i=0; i < N; i++){
            out[i] = 0;
            for(l=0; l < M; l++)
              out[i] = out[i] + C[i+l*N]*cur_sta[l];
            for(l=0; l < K; l++)
              out[i] = out[i] + D[i+l*N]*inp[l];
            out[i] = out[i] % 2;
          }
          for(i=0; i < M; i++)
            cur_sta[i] = cur_sta_pre[starter + i*n_std_sta];
        }

        /*  loc_exp = sum(rem(code(trace_eli, :) + out', 2));   */
        loc_exp1 = 0;
        for(i=0; i < N; i++)
          loc_exp1 = loc_exp1 + (code[trace_eli-1+i*leng] + out[i]) % 2;
	
#ifdef MATLAB_MEX_FILE
        if( plot_flag != 0 && plot_flag_test !=0 ){
          int nlhs, nrhs, aaa;
          Matrix *plhs4[1], *prhs4[1];

          nlhs = 1;
          nrhs = 1;
          prhs4[0] = V1;
          mxGetPr(prhs4[0])[0] = fig_position;
          mxGetPr(prhs4[0])[1] = leng;
          mxGetPr(prhs4[0])[2] = starter;
          mxGetPr(prhs4[0])[3] = num_state;
          mxGetPr(prhs4[0])[4] = dec;
          mxGetPr(prhs4[0])[5] = plot_flag;
          mxGetPr(prhs4[0])[6] = x[0];

          mexCallMATLAB(nlhs,plhs4,nrhs,prhs4,"sviplot4");

          mxFreeMatrix(plhs4[0]);
        }    
#endif

        starter = dec;
        for(i=0; i < K; i++)
          Y1[i] = inp[i];
        Y1[K] = loca_exp1;    
      }else{  /* if (trace_flag != 0 ) */
        for(i=0; i < K+1; i++)
          Y1[i] = 0;
        starter = 0;
      }  /* the end of "if (trace_flag != 0 )" */
      
      trace_num = trace_num % leng;

      x[1] = (double)fig_position;
      x[2] = (double)trace_num;
      x[3] = (double)trace_flag;
      x[4] = (double)starter;
      x[5] = (double)plot_flag_test;
      x[6] = (double)initial_flag;
      for(i=0; i < leng*PowPowM; i++)
    	x[i+7] = (double)expense1[i];
      for(i=0; i < leng*PowPowM; i++)
    	x[i+7+leng*PowPowM] = (double)solu[i];
      for(i=0; i < leng*N; i++)
    	x[i+7+2*leng*PowPowM] = (double)code[i];
      for(i=0; i < K+1; i++)
    	x[i+7+2*leng*PowPowM+leng*N] = (double)Y1[i];
   	}/* the end of " if( expen_flag != 0 )" */
  } /*  the end of "if ( u[N] >= 0.2 ) */
}

/*
 * mdlDerivatives - compute the derivatives
 *
 * In this function, you compute the S-function block's derivatives.
 * The derivatives are placed in the dx variable.
 */
static void mdlDerivatives(dx, x, u, S, tid)
     double *dx, *x, *u; 
     SimStruct *S;
     int tid;
{
}

/*
 * mdlTerminate - called when the simulation is terminated.
 *
 * In this function, you should perform any actions that are necessary
 * at the termination of a simulation.  For example, if memory was allocated
 * in mdlInitializeConditions, this is the place to free it.
 */
static void mdlTerminate(S)
     SimStruct *S;
{
}

#ifdef      MATLAB_MEX_FILE    /* Is this file being compiled as a MEX-file? */
#include "simulink.c"      /* MEX-file interface mechanism */
#else
#include "cg_sfun.h"       /* code generation registration function */
#endif