uep_qpsk.c

一些纠错编码的c程序实现

       } while( s >= 1);
       x1 = u1*sqrt((-2.0*log(s))/s);
       x2 = u2*sqrt((-2.0*log(s))/s);
       /* Received (complex) signal */
       rx[i][j] = amp*cos(2*((double)x[j])*pi/4.0) + x1 ;
       ry[i][j] = amp*sin(2*((double)x[j])*pi/4.0) + x2 ; 
      }




/***************************************************************************
 *                                                                         *
 *                             DECODING                                    *
 *                                                                         *
 ***************************************************************************/



/* PART (I): Compare with (1111111) --> decode using the lower subtrellis */
 


/* initialize path metrics */
   for(i1=1;i1<=7;i1++)
      {
      for(jj=1;jj<=trellis3[i1].count;jj++)
         {
         trellis_metric1[i1].branch[jj].final = 0.0;   
         trellis_metric1[i1].branch[jj].f[0] = 0.0;   
         trellis_metric1[i1].branch[jj].f[1] = 0.0;   
         } 
      }
 
     for(j=1;j<=7;j++)
        { 
         jj = j - 1;
         for(k=1;k<=trellis3[j].count;k++)
           {
 
           /* compute metrics for coset[0][1] */
           cosetvalue = trellis3[j].outbit1[k].coset[0][1];
 
           /* squared Euclidean distance to received signal */
           trellis_metric1[j].branch[k].distance[0] = 
               (rx[i][jj] - q[0][cosetvalue]) *
                    (rx[i][jj] - q[0][cosetvalue])
             + (ry[i][jj] - q[1][cosetvalue]) *
                    (ry[i][jj] - q[1][cosetvalue]);
 
           trellis_metric1[j].branch[k].final 
                    = trellis_metric1[j].branch[k].distance[0];
		   trellis_metric1[j].branch[k].c12est[0]
                    = trellis3[j].outbit1[k].coset[0][1];
		   trellis_metric1[j].branch[k].c12est[1]
                    = trellis3[j].outbit1[k].coset[0][1];
 
           }
         } 

/* Decode a sequence in the lower subtrellis by the VITERBI ALGORITHM */

   surv_final1[0].metric[0] = 0.0;
   for(j=1;j<=7;j++)
       {
       for(iii =0;iii<=1; iii++)
         { suv_count1[iii] = 0; }
       i1 = j -1;
       for(k=1;k<=trellis3[j].count;k++)
          {
           n = trellis3[j].final_state[k] ; 
           suv_count1[n] = suv_count1[n] + 1;
           m = trellis3[j].initial_state[k];
           if(suv_count1[n] == 1)
              { 
              suv1[j].surv_metric[n][0]
                       = surv_final1[i1].metric[m] + 
                             trellis_metric1[j].branch[k].final;
              suv1[j].surv_ini[n][0] = trellis3[j].initial_state[k];
              suv1[j].surv_out[n][0] = trellis3[j].final_state[k];
              suv1[j].est[n][0][0] = trellis_metric1[j].branch[k].c12est[0];
              suv1[j].est[n][0][1] = trellis_metric1[j].branch[k].c12est[1];

              surv_final1[j].metric[n] = suv1[j].surv_metric[n][0];
              surv_final1[j].ini[n] = suv1[j].surv_ini[n][0];
              surv_final1[j].final[n] = suv1[j].surv_out[n][0];
              surv_final1[j].est[n][0] = suv1[j].est[n][0][0]; 
              surv_final1[j].est[n][1] = suv1[j].est[n][0][1]; 
              } 
          else 
              { 
              suv1[j].surv_metric[n][1]
                 		= surv_final1[i1].metric[m] + 
                     		 trellis_metric1[j].branch[k].final;
              suv1[j].surv_ini[n][1] = trellis3[j].initial_state[k];
              suv1[j].surv_out[n][1] = trellis3[j].final_state[k];
              suv1[j].est[n][1][0] = trellis_metric1[j].branch[k].c12est[0];
              suv1[j].est[n][1][1] = trellis_metric1[j].branch[k].c12est[1];

              surv_final1[j].metric[n] = suv1[j].surv_metric[n][1];
              surv_final1[j].ini[n] = suv1[j].surv_ini[n][1];
              surv_final1[j].final[n] = suv1[j].surv_out[n][1];
              surv_final1[j].est[n][0] = suv1[j].est[n][1][0]; 
              surv_final1[j].est[n][1] = suv1[j].est[n][1][1]; 
              } 
          if(suv_count1[n] == 2)
              { 
              if( suv1[j].surv_metric[n][1] >= suv1[j].surv_metric[n][0])
                 {   
                 surv_final1[j].metric[n] = suv1[j].surv_metric[n][0];
                 surv_final1[j].ini[n] = suv1[j].surv_ini[n][0];
                 surv_final1[j].final[n] = suv1[j].surv_out[n][0];
                 surv_final1[j].est[n][0] = suv1[j].est[n][0][0]; 
                 surv_final1[j].est[n][1] = suv1[j].est[n][0][1]; 
                 }
              else 
                 {   
                 surv_final1[j].metric[n] = suv1[j].surv_metric[n][1];
                 surv_final1[j].ini[n] = suv1[j].surv_ini[n][1];
                 surv_final1[j].final[n] = suv1[j].surv_out[n][1];
                 surv_final1[j].est[n][0] = suv1[j].est[n][1][0]; 
                 surv_final1[j].est[n][1] = suv1[j].est[n][1][1]; 
                 }
             } 
          }
    }

    min_state = n;
    ii = 7;

/* Total metric of the surviving sequence in lower subtrellis (coset[0][1]) */
    final1 = surv_final1[7].metric[min_state];

    while (ii >= 1 )
       { 
       x_estimate[ii].est[0] = surv_final1[ii].est[min_state][0];
       x_estimate[ii].est[1] = surv_final1[ii].est[min_state][1];
       x_estimate[ii].est[2] = 1;

/* INVERSE MAPPING: Recover bits from QPSK signal */
 
				test = x_estimate[ii].est[1];
				x_estimate[ii].est[1] = 0;
				if ( (test == 2) || (test == 3) )
					x_estimate[ii].est[1] = 1;
 
       nn = min_state;
       min_state = surv_final1[ii].ini[nn];
       ii--;
       }




/* PART (II): Compare with (0000000) --> decode using the upper subtrellis */


 
/* initialize path metrics */
   for(i1=1;i1<=7;i1++)
      {
      for(jj=1;jj<=trellis3[i1].count;jj++)
         {
         trellis_metric[i1].branch[jj].final = 0.0;   
         trellis_metric[i1].branch[jj].f[0] = 0.0;   
         trellis_metric[i1].branch[jj].f[1] = 0.0;   
         } 
      }
 
     for(j=1;j<=7;j++)
        { 
        jj = j - 1;
        for(k=1;k<=trellis3[j].count;k++)
           {
 
           /* compute metrics for coset[0][0] */
           cosetvalue = trellis3[j].outbit1[k].coset[0][0];
 
           /* squared Euclidean distance to received signal */
           trellis_metric[j].branch[k].distance[0] =
               (rx[i][jj] - q[0][cosetvalue]) *
                   (rx[i][jj] - q[0][cosetvalue])
             + (ry[i][jj] - q[1][cosetvalue]) *
                   (ry[i][jj] - q[1][cosetvalue]);
 
     	   trellis_metric[j].branch[k].final 
                 = trellis_metric[j].branch[k].distance[0];
     	   trellis_metric[j].branch[k].c12est[0] 
                 				= trellis3[j].outbit1[k].coset[0][0];
           trellis_metric[j].branch[k].c12est[1]
                                = trellis3[j].outbit1[k].coset[0][0];
    	   }
        } 

/* Decode a sequence in the upper subtrellis by the VITERBI ALGORITHM */

   surv_final[0].metric[0] = 0.0;
   for(j=1;j<=7;j++)
       {
      for(iii =0;iii<=1; iii++)
         { suv_count[iii] = 0; }
         i1 = j -1;
         for(k=1;k<=trellis3[j].count;k++)
          {
           n = trellis3[j].final_state[k] ; 
           suv_count[n] = suv_count[n] + 1;
           m = trellis3[j].initial_state[k];
           if(suv_count[n] == 1)
           { 
            suv[j].surv_metric[n][0]
                 = surv_final[i1].metric[m] + 
                     trellis_metric[j].branch[k].final;
            suv[j].surv_ini[n][0] = trellis3[j].initial_state[k];
            suv[j].surv_out[n][0] = trellis3[j].final_state[k];
            suv[j].est[n][0][0] = trellis_metric[j].branch[k].c12est[0];
            suv[j].est[n][0][1] = trellis_metric[j].branch[k].c12est[1];

            surv_final[j].metric[n] = suv[j].surv_metric[n][0];
            surv_final[j].ini[n] = suv[j].surv_ini[n][0];
            surv_final[j].final[n] = suv[j].surv_out[n][0];
            surv_final[j].est[n][0] = suv[j].est[n][0][0]; 
            surv_final[j].est[n][1] = suv[j].est[n][0][1]; 
           } 
          else 
           { 
            suv[j].surv_metric[n][1]
                 = surv_final[i1].metric[m] + 
                     trellis_metric[j].branch[k].final;
            suv[j].surv_ini[n][1] = trellis3[j].initial_state[k];
            suv[j].surv_out[n][1] = trellis3[j].final_state[k];
            suv[j].est[n][1][0] = trellis_metric[j].branch[k].c12est[0];
            suv[j].est[n][1][1] = trellis_metric[j].branch[k].c12est[1];

            surv_final[j].metric[n] = suv[j].surv_metric[n][1];
            surv_final[j].ini[n] = suv[j].surv_ini[n][1];
            surv_final[j].final[n] = suv[j].surv_out[n][1];
            surv_final[j].est[n][0] = suv[j].est[n][1][0]; 
            surv_final[j].est[n][1] = suv[j].est[n][1][1]; 
           } 
           if(suv_count[n] == 2)
             { 
               if( suv[j].surv_metric[n][1] >= suv[j].surv_metric[n][0])
                 {   
                 surv_final[j].metric[n] = suv[j].surv_metric[n][0];
                 surv_final[j].ini[n] = suv[j].surv_ini[n][0];
                 surv_final[j].final[n] = suv[j].surv_out[n][0];
                 surv_final[j].est[n][0] = suv[j].est[n][0][0]; 
                 surv_final[j].est[n][1] = suv[j].est[n][0][1]; 
                 }
              else 
                 {   
                 surv_final[j].metric[n] = suv[j].surv_metric[n][1];
                 surv_final[j].ini[n] = suv[j].surv_ini[n][1];
                 surv_final[j].final[n] = suv[j].surv_out[n][1];
                 surv_final[j].est[n][0] = suv[j].est[n][1][0]; 
                 surv_final[j].est[n][1] = suv[j].est[n][1][1]; 
                 }
             } 
         }
    }
    min_state = n;
    ii = 7;

/* Total metric of the surviving sequence in upper subtrellis (coset[0][0]) */
    final0 = surv_final[7].metric[min_state];

    while (ii >= 1 )
       { 
       x_estimate1[ii].est[0] = surv_final[ii].est[min_state][0];
       x_estimate1[ii].est[1] = surv_final[ii].est[min_state][1];
       x_estimate1[ii].est[2] = 0;

/* INVERSE MAPPING: Recover bits from QPSK signal */
 
       test = x_estimate1[ii].est[1];
       x_estimate1[ii].est[1] = 0;
       if ( (test == 2) || (test == 3) )
          x_estimate1[ii].est[1] = 1;
          nn = min_state;
          min_state = surv_final[ii].ini[nn];
          ii--;
          }

/* Decoding of a block finished. Count the number of decoding errors */

    if( final1 >= final0 )  /* if sequence in the upper subtrellis */
       {
       if(x_estimate1[1].est[2] != code1[i].first)
          {
          first_error += 1.0;
          }
       for(kk=1;kk<=7;kk++)
          {
          k = kk - 1;
          if(x_estimate1[kk].est[1] != code1[i].second[k])
             {
             second_error += 1.0;
             }
          }
        }
        else                   /* if sequence in the lower subtrellis */
        {
         if(x_estimate[1].est[2] != code1[i].first)
            {
            first_error += 1.0;
            }
         for(kk=1;kk<=7;kk++)
            {
            k = kk - 1;
            if(x_estimate[kk].est[1] != code1[i].second[k])
               {
               second_error += 1.0;
               }
            }
        }



     }

     error += (second_error + first_error);
     error1 += first_error;
     error2 += second_error;

#ifdef DEBUG2
printf("%8d %10.0f %10.0f %10.0f\n", ifinal, first_error, second_error, error);
fflush(stdout) ;
#endif



   }



   biterror = error/(NUMSIM*(7*200)); 
   biterr1 = error1/(NUMSIM*200); 
   biterr2 = error2/(NUMSIM*(6*200)); 

   first_error = NUMSIM*(7*200);
   printf("Done... A total of %e bits were transmitted\n", first_error);
   printf("error,error1,error2 = %20.0f %20.0f %20.0f\n",error, error1, error2);

   printf("Pe  = %20.17e \n", biterror);
   printf("Pe1 = %20.17e \n", biterr1);
   printf("Pe2 = %20.17e \n", biterr2);

   snr = 10*log10(amp*amp/2.0);
   printf("Eb/No = %20.10f \n", snr);
}