vitcore.c

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

  /* considering the different type of variable 'expense' */
  if(expen_flag != 0){
    for(i=0; i < leng*PowPowM; i++){
      expense[i] = NaN;
      solu[i] = 0;
    }
    for(i=0; i < n_std_sta; i++)
      expense[leng-1+i*leng] = 0;
    for(i=0; i < PowPowM; i++)
      solution[i] = expense[i*leng];

  } else {
    for(i=0; i < leng*PowPowM; i++){
      expense1[i] = NaN;
      solu[i] = 0;
    }
    for(i=0; i < n_std_sta; i++)
      expense1[leng-1+i*leng] = 0;
    for(i=0; i < PowPowM; i++)
      solution1[i] = expense1[i*leng];
  }

  /*pre_state = 1;
   *inp_pre = de2bi([0:K2-1]', K);
   *cur_sta_pre = de2bi([0:n_std_sta-1], M);
   *starter = 0;
   *msg = [];
   *expen = [];
   *codd = [];
   */
  len_pre_state = 1;
  pre_state[0] = 1;

  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);

  starter = 0;
  msg = mxGetPr(plhs[0]=mxCreateFull(n_code, K, 0));
  expen = mxGetPr(plhs[1]=mxCreateFull(n_code, 1, 0));
  codd = mxGetPr(plhs[2]=mxCreateFull(n_code, N, 0));
  n_msg = 0;
  len_expen = 0;
  n_codd = 0;
  /*for i = 1 : n_code
   *    % make room for one more storage.
   *    trace_pre = rem(i-2+leng, leng) + 1;  % previous line of the trace.
   *    trace_num = rem(i-1, leng) + 1;       % current line of the trace.
   *    expense(trace_num,:) = solution;
   */
  for(ii=1; ii <= n_code; ii++){
    trace_pre = (ii-2+leng) % leng + 1; /* previous line of the trace. */
    trace_num = (ii-1) % leng + 1;      /* current line of the trace. */
    for( i=0; i < PowPowM; i++){
      if ( expen_flag != 0 )
        expense[trace_num-1+i*leng] = solution[i];
      else
        expense1[trace_num-1+i*leng] = solution1[i];
    }
    /*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(i, num_N)));
     *        end;
     *        for num_k = 1 : K2
     *            inp = inp_pre(num_k, :)';
     *            if isempty(C)
     *                tran_indx = pre_state(j) + (num_k -1) * K2;
     *                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;
     *            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;
     */
    if( expen_flag != 0){
      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[i*n_tran_prob] <= code[ii-1 + num_N*n_code] )
              max = i+1;
          }
	      expen_work[num_N] = max;
	    }
        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;
          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[1 + n_tran_prob*(expen_work[ expenOut[i] ]-1) ];
          for(i=0; i < lenIndx1; i++)
            loca_exp = loca_exp + tran_prob[2 + n_tran_prob*(expen_work[expenOut[i+lenIndx0]]-1)];
          tmp = (nex_sta_de - 1) * n_std_sta + pre_state[j] - 1; /* minus 1 for index shift */
          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;
          }
        }
      }

    /*    else
     *        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_k = 1 : K2
     *                inp = inp_pre(num_k, :)';
     *                if isempty(C)
     *                    tran_indx = pre_state(j) + (num_k -1) * K2;
     *                    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;
     *                nex_sta_de = bi2de(nex_sta') + 1;
     *                loca_exp = sum(rem(code(i, :) + out', 2));
     *                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;
     *    end;
     */
    } else {
      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 + ((int)code[ii-1+i*n_code] + 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;
            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;
          }
        }
      }
	}
	
    /*    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;
     *    % go back one step to re-arrange the lines.
     */
    len_aft_state = 0;
    for(j2=0; j2 < n_std_sta; j2++){
      numnotnan = 0;
      for(i=0; i < n_std_sta; i++){
        if ( expen_flag != 0 ){
          if( expense[trace_num-1+i*leng+j2*leng*n_std_sta] <=0 )
            numnotnan ++;
        } else {
          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++;
      }
    }

    if( len_plot_flag > 0 ){
      Matrix *lhs1[1], *rhs1[4];

      rhs1[0] = mxCreateFull(4, 1, 0);
      mxGetPr(rhs1[0])[0] = num_state;
      mxGetPr(rhs1[0])[1] = ii;
      mxGetPr(rhs1[0])[2] = expen_flag;
      mxGetPr(rhs1[0])[3] = n_code;
      rhs1[1] = mxCreateFull(1, PowPowM, 0);
      for(i=0; i < PowPowM; i++){
        if(expen_flag != 0){
          if ( expense[trace_num-1 + i*leng] > 0 )
            mxGetPr(rhs1[1])[i] = mexGetNaN();
          else
            mxGetPr(rhs1[1])[i] = expense[trace_num-1 + i*leng];
        }else{
          if ( expense1[trace_num-1 + i*leng] < 0 )
            mxGetPr(rhs1[1])[i] = mexGetNaN();
          else
            mxGetPr(rhs1[1])[i] = expense1[trace_num-1 + i*leng];
        }
      }
      rhs1[2] = mxCreateFull(1, len_aft_state, 0);
      for(i=0; i < len_aft_state; i++)
        mxGetPr(rhs1[2])[i] = aft_state[i];
      rhs1[3] = mxCreateFull(1, len_plot_flag, 0);
      for(i=0; i<len_plot_flag; i++)
        mxGetPr(rhs1[3])[i] = plot_flag[i];
        
      mexCallMATLAB(1, lhs1, 4, rhs1, "vitplot1");

      mxFreeMatrix(lhs1[0]);
    }
    /*    sol = expense(trace_num,:);
     *    % decision making.
     *    if i >= leng
     *        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( expen_flag != 0 ){
      for( i=0; i < PowPowM; i++)
        sol[i] = expense[trace_num-1 + i*leng];
    } else {
      for( i=0; i < PowPowM; i++)
        sol1[i] = expense1[trace_num-1 + i*leng];
    }
        
    if( ii >= leng ){
      trace_eli = trace_num % leng + 1;
	  if( expen_flag != 0 ){
        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];
          }
        }
      } else {
        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];
          }
        }
      }

      /*        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);