sviterbi.c

通信系统的matlab仿真程序

        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

        /*  [n_tran_prob_tmp, m_tran_prob_tmp] = size(tran_prob_tmp);
         *  if n_tran_prob_tmp == 3
         *    expen_flag = 1;     % expense flag == 0; BSC code. == 1, with tran_prob_tmp.
         *    expen_work = zeros(1, N); % work space.
         *    tran_prob_tmp(2:3, :) = log10(tran_prob_tmp(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);
         */
      
      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
         */
#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 = 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_tmp(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_tmp(2,expen_work(expen_0)) 0])...
         *                  +sum([tran_prob_tmp(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 %%%%% */
#ifdef MATLAB_MEX_FILE
      if (x[0] > 0 && plot_flag_test != 0 ){
        int nlhs, nrhs;
        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( expense[i] == NaN )
            mxGetPr(prhs3[2])[i] = mexGetNaN();
          else
            mxGetPr(prhs3[2])[i] = (double)expense[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
        /*  % 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;
         *        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_tmp(2,expen_work(expen_0)) 0])...
         *               +sum([tran_prob_tmp(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)];
	
#ifdef MATLAB_MEX_FILE
        if( plot_flag != 0 && plot_flag_test !=0 ){
          int nlhs, nrhs;
          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;