sviterba.c

Proakis《contemporarycommunication systems using matlab》matlab源代码

         *        out = B(tran_indx, :)';
         *    else
         *        cur_sta = cur_sta_pre(starter+1, :)';
         *        out = rem(C * cur_sta + D * inp,2);
         *    end;
         */
	/* here, expen_flag != 0 */ 
	
	for(i=0; i < PowPowM; i++){
	  if( sol[i] <= 0 ){
	    loc_exp = sol[i];
	    i = PowPowM;
	  }
	}
	for(i=0; i < PowPowM; i++){
	  if( sol[i] <= 0 && loc_exp < sol[i])
	    loc_exp = sol[i];
	}
	for(i=0; i < PowPowM; i++){
	  if( sol[i] == loc_exp ){
	    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];
	}
        /*    if expen_flag
         *      % find the expense by the transfer probability
         *      expen_0 = find(out' <= 0.5);
         *      expen_1 = find(out' > 0.5);
         *      loc_exp = sum([tran_prob(2,expen_work(expen_0)) 0])...
         *               +sum([tran_prob(3,expen_work(expen_1)) 0]);
         *    else
         *      loc_exp = sum(rem(code(trace_eli, :) + out', 2));
         *    end
         */
	lenIndx0= 0;
	for(i=0; i < N; i++){
	  if( out[i] <= 0.5 ){
	    expenOut[lenIndx0] = i;
	    lenIndx0++;
	  }
	}
	lenIndx1 = 0;
	for(i=0; i < N; i++){
	  if( out[i] > 0.5 ){
	    expenOut[lenIndx1+lenIndx0] = i;
	    lenIndx1++;
	  }
	}
	loc_exp = 0;
	for(i=0; i < lenIndx0; i++)
	  loc_exp = loc_exp + tran_prob_tmp[1+n_tran_prob*(expen_work[expenOut[i]]-1)];
	for(i=0; i < lenIndx1; i++)
	  loc_exp = loc_exp + tran_prob_tmp[2+n_tran_prob*(expen_work[expenOut[i+lenIndx0]]-1)];
	

        /*    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;
      for(i=0; i < leng*PowPowM; i++)
	x[i+5] = expense[i];
      for(i=0; i < leng*PowPowM; i++)
	x[i+5+leng*PowPowM] = (double)solu[i];
      for(i=0; i < leng*N; i++)
	x[i+5+2*leng*PowPowM] = (double)code[i];
      for(i=0; i < K+1; i++)
	x[i+5+2*leng*PowPowM+leng*N] = Y[i];
      /*  the end of "if (expen_flag != 0 ) */
    }else{ /* tran_prob is not 3-row matrix */
      /* In this kind of TRAN_PROB, the type of all of variables is integer */
      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;
      
      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];
      loc_tmp = 6;
      
      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];
      
      fig_position = x[1] + 1;
      if( x[0] > 0 && ((fig_position-leng)%(plot_flag - leng) == 0) && fig_position >= plot_flag && plot_flag_test != 0 ){
	
      }
      
      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++;
	}
      }

      /*  %%%%% begin plot related %%%%% */
      
      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;
	
	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;
      for(i=0; i < leng*PowPowM; i++)
	x[i+5] = (double)expense1[i];
      for(i=0; i < leng*PowPowM; i++)
	x[i+5+leng*PowPowM] = (double)solu[i];
      for(i=0; i < leng*N; i++)
	x[i+5+2*leng*PowPowM] = (double)code[i];
      for(i=0; i < K+1; i++)
	x[i+5+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