sviterba.c

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

   *  num_state = M;
   *  n_std_sta = 2^M;
   *  PowPowM = n_std_sta^2;
   *  K2 = 2^K;
   *  inp_pre = de2bi([0:K2-1]', K); 
   *  cur_sta_pre = de2bi([0:n_std_sta-1], M);
   */
  if ( n_tran_prob == 3 ){
    expen_flag = 1;
    NaN = 1; 
  }else{
    expen_flag = 0;
    NaN = -1;
  }
  
  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;
    
    A = ssGetIWork(S);
    B = A + M*M;
    C = B + M*K;
    D = C + N*M;
    
    /* 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)mxGetPr(TRAN_FUNC)[i+j*(M+N)];
    	if( i>M-1 && j<M )
          C[i+j*N-M] = (int)mxGetPr(TRAN_FUNC)[i+j*(M+N)];
    	if( i<M   && j>M-1 )
    	  B[i+j*M-M*M] = (int)mxGetPr(TRAN_FUNC)[i+j*(M+N)];
    	if( i>M-1 && j>M-1 )
    	  D[i+j*N-M*(N+1)] = (int)mxGetPr(TRAN_FUNC)[i+j*(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;
    
    /* Assignment */
    TRAN_A =    ssGetIWork(S);/*   !--- size of *TRAN_A */
    TRAN_B =    ssGetIWork(S) + (rowFunc-2);/* <--- size of *TRAN_B */
    A =         ssGetIWork(S) + 2*(rowFunc-2);/*    !----- size of *A */ 
    B =         ssGetIWork(S) + 2*(rowFunc-2) + (rowFunc-2)*M;/*    !----- size of *B */
    
    /* Get the input Matrix A, B */
    for(i=0; i < rowFunc-2; i++){
      TRAN_A[i] = (int)mxGetPr(TRAN_FUNC)[i+2];
      TRAN_B[i] = (int)mxGetPr(TRAN_FUNC)[rowFunc+i+2];
    }
    de2bi(TRAN_A, M, rowFunc-2, A);
    de2bi(TRAN_B, N, rowFunc-2, B);
  }
  
  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 ( u[N] >= 0.2 ){
    if( expen_flag != 0 ){  /* Tran_prob is a THREE rows Matrix */
      tran_prob_tmp    = ssGetRWork(S);
      expense     = ssGetRWork(S) + n_tran_prob*m_tran_prob;
      Y           = ssGetRWork(S) + n_tran_prob*m_tran_prob + leng*PowPowM;
      sol         = ssGetRWork(S) + n_tran_prob*m_tran_prob + leng*PowPowM + (K+1);
      expen_tmp   = ssGetRWork(S) + n_tran_prob*m_tran_prob + leng*PowPowM + (K+1) + PowPowM;
      tmpRwork    = expen_tmp + PowPowM;
      
      if( mxGetPr(TRAN_FUNC)[rowFunc*colFunc-1] < 0 )
	     solu = D + N*K;
      else
	     solu = B + N*(rowFunc-2);
      
      code        = solu + leng*PowPowM;        
      inp_pre     = code + leng*N;        /* allocate K*2^K for *inp_pre */
      cur_sta_pre = inp_pre + K2*K;   /*  M*2^M for *cur_sta_pre. */
      expen_work  = cur_sta_pre + M*n_std_sta;  /* allocate N for *expen_work */
      pre_state   = expen_work + N;     /* 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 */
      tmpIwork    = aft_state + 2*n_std_sta;
      
      inp_pre[0] = -1;
      cur_sta_pre[0] = -1;
      de2bi(inp_pre, K, K2, inp_pre);
      de2bi(cur_sta_pre, M, n_std_sta, cur_sta_pre);
      
      /*  [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.
       *    expen_work = zeros(1, N); % work space.
       *    tran_prob(2:3, :) = log10(tran_prob(2:3, :));
       *  else
       *    expen_flag = 0;
       *  end;
       */
      for(i=0; i < N; i++)
	     expen_work[i] = 0;
      
      for(i=0; i < m_tran_prob; i++){
	     tran_prob_tmp[i*n_tran_prob] = mxGetPr(TRAN_PROB)[i*n_tran_prob];
	     tran_prob_tmp[i*n_tran_prob+1] = log10(mxGetPr(TRAN_PROB)[i*n_tran_prob+1]);
	     tran_prob_tmp[i*n_tran_prob+2] = log10(mxGetPr(TRAN_PROB)[i*n_tran_prob+2]);
      }
      
      /*  % 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.
       *  starter = x(5);
       *  solu = zeros(leng, PowPowM);
       *  expense = solu;
       *  code = zeros(leng, N);
       *  loc_tmp  = 6;
       *  expense(:) = x(loc_tmp : leng*PowPowM + loc_tmp - 1); 
       *  solu(:) = x(loc_tmp + leng * PowPowM : 2 * leng * PowPowM + loc_tmp - 1);
       *  code(:) = x(loc_tmp + 2 * leng * PowPowM : loc_tmp + 2 * leng * PowPowM -1 + leng * N);
       */
      starter = (int)x[4];
      loc_tmp = 6;
      
      for(i=0; i < leng*PowPowM; i++)
    	expense[i] = 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(2) + 1;
       *  if (x(1) > 0) & (rem(fig_position-leng, plot_flag - leng) == 0)
       *             & (fig_position >= plot_flag) & plot_flag_test
       *  see also sviplot2.m
       */
      fig_position = x[1] + 1;
      
      /*  trace_num = x(3) + 1;
       *  trace_flag = x(4);
       *  code(trace_num, :) = u(1:length(u)-1);
       *
       *  if (trace_flag == 0) & (trace_num == leng)
       *      trace_flag = 1;
       *  end;
       *  trace_pre = rem(trace_num - 2 + leng, leng) + 1;
       */
      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;
      
      /*  % fill up trace and solut
       *  if (trace_flag == 0) & (trace_num == 1)
       *      pre_state = starter + 1;
       *  else
       *      pre_state = [];
       *      for j2 = 1 : n_std_sta
       *          if max(~isnan(expense(trace_pre, [1-n_std_sta : 0] + j2*n_std_sta)))
       *              pre_state = [pre_state, j2];
       *          end;
       *      end;
       *  end;
       */
      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( expense[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++;
	       }
	     }
      }
      
      /*  expense(trace_num, :) = expense(trace_num,:) * NaN;
       *  if expen_flag
       *      for j = 1 : length(pre_state)
       *          jj = pre_state(j) - 1;           % index number - 1 is the state.
       *          cur_sta = cur_sta_pre(pre_state(j),:)';
       *          indx_j = (pre_state(j) - 1) * n_std_sta;
       *          for num_N = 1 : N
       *              expen_work(num_N) = max(find(tran_prob(1,:) <= code(trace_num, num_N)));
       *          end;
       *          for num_K = 1 : K2
       *              inp = inp_pre(num_K, :)';
       *              if isempty(C)
       *                  tran_indx = pre_state(j) + (num_K -1) * n_std_sta;
       *                  nex_sta = A(tran_indx, :)';
       *                  out = B(tran_indx, :)';
       *              else
       *                  out = rem(C * cur_sta + D * inp,2);
       *                  nex_sta = rem(A * cur_sta + B * inp, 2);
       *              end;
       */
      for(i=0; i < PowPowM; i++)
    	expense[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_N=0; num_N < N; num_N++){
          max = 0;
          for(i=0; i < m_tran_prob; i++){
            if( tran_prob_tmp[i*n_tran_prob] <= code[trace_num-1+num_N*leng] )
              max = i;
          }
	      expen_work[num_N] = max + 1;
	    }
        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;
	        }
	     }
	  /*              nex_sta_de = bi2de(nex_sta') + 1;
	   *              % find the expense by the transfer probability
	   *              expen_0 = find(out' <= 0.5);
	   *              expen_1 = find(out' > 0.5);
	   *              loca_exp = sum([tran_prob(2,expen_work(expen_0)) 0])...
	   *                  +sum([tran_prob(3,expen_work(expen_1)) 0]);
	   *              tmp = (nex_sta_de-1)*n_std_sta + pre_state(j);
	   *              if isnan(expense(trace_num, tmp))
	   *                  expense(trace_num, tmp) = loca_exp;
	   *                  solu(trace_num, nex_sta_de + indx_j) = num_K;   
	   *              elseif expense(trace_num, tmp) < loca_exp
	   *                  expense(trace_num, tmp) = loca_exp;
	   *                  solu(trace_num, nex_sta_de + indx_j) = num_K;
	   *              end;
	   *          end;
	   *      end;
	   */
	  bi2de(nex_sta,1, M, &nex_sta_de);
	  nex_sta_de = nex_sta_de + 1;
	  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++;
	    }
	  }
	  loca_exp = 0;
	  for(i=0; i < lenIndx0; i++)
	    loca_exp = loca_exp + tran_prob_tmp[1 + n_tran_prob*(expen_work[ expenOut[i] ]-1) ];
	  for(i=0; i < lenIndx1; i++)
	    loca_exp = loca_exp + tran_prob_tmp[2 + n_tran_prob*(expen_work[expenOut[i+lenIndx0]]-1)];
	  tmp = (nex_sta_de - 1) * n_std_sta + pre_state[j] - 1;
	  if( expense[trace_num - 1 + tmp*leng] > 0 ){
	    expense[trace_num - 1 + tmp*leng] = loca_exp;
	    solu[trace_num - 1 + leng*(nex_sta_de+indx_j-1)] = num_K + 1;
	  }else if( expense[trace_num - 1 + tmp*leng] < loca_exp ){
	    expense[trace_num - 1 + tmp*leng] = loca_exp;
	    solu[trace_num - 1 + leng*(nex_sta_de+indx_j-1)] = num_K + 1;
	  }
	}
      }
      /*  aft_state = [];
       *  for j2 = 1 : n_std_sta
       *      if max(~isnan(expense(trace_num, [1-n_std_sta : 0] + j2*n_std_sta)))
       *          aft_state = [aft_state, j2];
       *      end
       *  end;
       */
      len_aft_state = 0;
      for(j2=0; j2 < n_std_sta; j2++){
    	numnotnan = 0;
    	for(i=0; i < n_std_sta; i++){
          if( expense[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 %%%%% */
      
      /*  % decision making.
       *  if trace_flag
       *    sol = expense(trace_num,:);
       *    trace_eli = rem(trace_num, leng) + 1;
       *    % strike out the unnecessary.
       *    for j_k = 1 : leng - 2
       *      j_pre = rem(trace_num - j_k - 1 + leng, leng) + 1;
       *      sol = vitshort(expense(j_pre, :), sol, n_std_sta, expen_flag);
       *    end;
       *
       *    tmp = (ones(n_std_sta,1) * expense(trace_eli, [starter+1:n_std_sta:PowPowM]))';
       *    sol = sol + tmp(:)';
       */
      if( trace_flag != 0 ){
        trace_eli = (trace_num % leng) + 1;
	
        for( i=0; i < PowPowM; i++)
          sol[i] = expense[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_tmp[i] = expense[j_pre-1 + i*leng];
          shortdbl(expen_tmp, sol, n_std_sta, tmpRwork, tmpIwork);
        }
        for(j=0; j < n_std_sta; j++){
          for(i=0; i < n_std_sta; i++){
            if( expense[trace_eli-1+(starter+i*n_std_sta)*leng] > 0 )
              sol[i+j*n_std_sta] = 1;
	        else
	          sol[i+j*n_std_sta] = sol[i+j*n_std_sta] + expense[trace_eli-1+(starter+i*n_std_sta)*leng];
          }
        }

        /*    if expen_flag
         *      loc_exp =  max(sol(find(~isnan(sol))));   
         *    else
         *      loc_exp =  min(sol(find(~isnan(sol))));
         *    end
         *    dec = find(sol == loc_exp);
         *    dec = dec(1);
         *    dec = rem((dec - 1), n_std_sta);
         *    num_K = solu(trace_eli, starter*n_std_sta+dec+1);
         *    inp = inp_pre(num _K, :)';
         *    if isempty(C)
         *        tran_indx =  starter+1 + (num_K -1) * n_std_sta;