vitcore.c

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

       *        dec = rem((dec - 1), n_std_sta);
       *
       *        inp = de2bi(solu(trace_eli, starter*n_std_sta+dec+1)-1, K);
       *        if isempty(C)
       *            tran_indx = starter + 1 + (num_k -1) * n_std_sta;
       *            out = B(tran_indx, :)';
       *        else
       *            cur_sta = cur_sta_pre(starter+1, :)';
       *            out = rem(C * cur_sta + D * inp(:),2);
       *        end;
       */
      if ( expen_flag != 0 ){  /* 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;
      } else {
        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)] - 1;
      de2bi(&num_K, K, one, inp);

      if( len_C == 0 ){
        tran_indx = starter + 1 + num_K*n_std_sta;
        for(i=0; i < N; i++)
          out[i] = B[tran_indx-1+i*(rowFunc-2)];
      } else {
        for(i=0; i < M; i++)
          cur_sta[i] = cur_sta_pre[starter + i*n_std_sta];
        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;
        }
      }

      /*        msg = [msg; inp];
       *        codd = [codd; out'];
       *        [n_msg, m_msg] = size(msg);
       *        if nargout > 1
       *            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(n_msg, :) + out', 2));
       *            end
       *            expen = [expen; loc_exp];
       *        end;
       */
      for( i=0; i < K; i++)
        msg[n_msg+i*n_code] = inp[i];
      n_msg++ ;
      for( i=0; i < N; i++)
        codd[n_codd+i*n_code] = out[i];
      n_codd++ ;

      /* calculate the second output 'expen' */
      if( expen_flag != 0 ){
        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[1+n_tran_prob*(expen_work[expenOut[i]]-1)];
        for(i=0; i < lenIndx1; i++)
          loc_exp = loc_exp + tran_prob[2+n_tran_prob*(expen_work[expenOut[i+lenIndx0]]-1)];
        expen[len_expen] = loc_exp;
      }else{
        loca_exp1 = 0;
        for(i=0; i < N; i++)
          loca_exp1 = loca_exp1 + ((int)code[n_msg-1+i*n_code] + out[i]) % 2;
        expen[len_expen] = loca_exp1;
      }
      len_expen++;

      if( len_plot_flag > 0 ){
        Matrix *lhs2[1], *rhs2[1];

        rhs2[0] = mxCreateFull(5, 1, 0);
        mxGetPr(rhs2[0])[0] = n_msg;
        mxGetPr(rhs2[0])[1] = starter;
        mxGetPr(rhs2[0])[2] = dec;
        mxGetPr(rhs2[0])[3] = plot_flag[0];
        mxGetPr(rhs2[0])[4] = M;

        mexCallMATLAB(1, lhs2, 1, rhs2, "vitplot2");

        mxFreeMatrix(lhs2[0]);
      }
      /*        starter = dec;
       *    end; %(if i >= leng)
       *    pre_state = aft_state;
       *    aft_state = [];
       *end;
       */
      starter = dec;
    } /* the end of  (if i >= leng) */
    for( i=0; i < len_aft_state; i++ )
      pre_state[i] = aft_state[i];
    len_pre_state = len_aft_state;  
    len_aft_state = 0;
  }

  /*for i =  1 : leng-1
   *    sol = solution;
   *    trace_eli = rem(trace_num + i, leng) + 1;
   *    if i < leng-1
   *        sol(1:n_std_sta) = expense(trace_num, 1:n_std_sta);
   *        for j_k = 1 : leng - 2 - i
   *            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 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);
   *    else
   *        dec = 0;
   *    end;
   */
  for(ii=1; ii <= leng-1; ii++){
    if( expen_flag != 0 ){        /* here, expen_flag != 0 */ 
      for(i=0; i < PowPowM; i++)
        sol[i] = solution[i];
    } else {    
      for(i=0; i < PowPowM; i++)
        sol1[i] = solution1[i];
    }
    trace_eli = (trace_num + ii) % leng + 1;

    if( ii < leng-1 ){
      if( expen_flag != 0 ){        /* here, expen_flag != 0 */ 
        for(i=0; i < n_std_sta; i++)
          sol[i] = expense[trace_num-1 + i*leng];
        for(j_k=1; j_k <= leng-2-ii; 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] > 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];
          }
        }
        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;
          }
        }
      } else {
        for(i=0; i < n_std_sta; i++)
          sol1[i] = expense1[trace_num-1 + i*leng];
        for(j_k=1; j_k <= leng-2-ii; 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] < 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;
    } else {
      dec = 0;
	} /* end -- if( ii < leng-1 ) */

    /*    inp = de2bi(solu(trace_eli, starter*n_std_sta+dec+1)-1, K);
     *    cur_sta = de2bi(starter, num_state);
     *    out = rem(C*cur_sta' + D * inp', 2);
     *    msg = [msg; inp];
     *    codd = [codd; out'];
     *    [n_msg, m_msg] = size(msg);
     */
    num_K = solu[trace_eli-1+leng*(starter*n_std_sta+dec)] - 1;
    de2bi(&num_K, K, one, inp);
    de2bi(&starter, num_state, one, cur_sta);

    if( len_C == 0 ){
      tran_indx = starter + 1 + num_K*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 < K; i++)
      msg[n_msg+i*n_code] = inp[i];
    n_msg++ ;
    for( i=0; i < N; i++)
      codd[n_codd+i*n_code] = out[i];
    n_codd++ ;

    /*    if nargout > 1
     *        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(n_msg, :) + out', 2));
     *        end
     *        expen = [expen; loc_exp];
     *    end;
     */
    if(expen_flag != 0){
      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[1+n_tran_prob*(expen_work[expenOut[i]]-1)];
      for(i=0; i < lenIndx1; i++)
        loc_exp = loc_exp + tran_prob[2+n_tran_prob*(expen_work[expenOut[i+lenIndx0]]-1)];
      expen[len_expen] = loc_exp;
      len_expen ++;
    } else {
      loca_exp1 = 0;
      for(i=0; i < N; i++)
        loca_exp1 = loca_exp1 + ((int)code[n_msg-1 + i*n_code] + out[i]) % 2;
      expen[len_expen] = loca_exp1;
      len_expen ++;
    }

    if( len_plot_flag > 0 ){
      Matrix *lhs2[1], *rhs2[1];

      rhs2[0] = mxCreateFull(5, 1, 0);
      mxGetPr(rhs2[0])[0] = n_msg;
      mxGetPr(rhs2[0])[1] = starter;
      mxGetPr(rhs2[0])[2] = dec;
      mxGetPr(rhs2[0])[3] = plot_flag[0];
      mxGetPr(rhs2[0])[4] = M;
        
      mexCallMATLAB(1, lhs2, 1, rhs2, "vitplot2");

      mxFreeMatrix(lhs2[0]);
    }
    /*    starter = dec;
     *end;
     */
    starter = dec;
  }

  /*% cut the extra message length
   *[n_msg, m_msg] = size(msg);
   *msg(n_msg-M+1:n_msg, :) = [];
   *if plot_flag(1)
   *    set(xx,'Color',[0 1 1]);
   * hold off
   *end;
   *% end of VITERBI.M
   */
  mxSetM(plhs[0], n_msg-M);
  mxSetN(plhs[0], K);
  for(i=0; i < n_msg-M; i++){
    for(j=0; j < K; j++)
      mxGetPr(plhs[0])[i + j*(n_msg-M)] = msg[i + j*n_code];
  }
  return;
}
/*--end of VITERBI.C --*/