uep_qpsk.c

一些纠错编码源代码

/*
   file:         multi1.c
   author:       Robert H. Morelos-Zaragoza
   revised:      October 5, 1994
  
   ===========================================================================
  
            SIMULATION OF A BLOCK QPSK MODULATION CODE OF LENGTH 7
                   WITH TWO LEVELS OF ERROR PROTECTION
  
                           AWGN CHANNEL MODEL
  
          SINGLE-STAGE MAXIMUM-LIKELIHOOD DECODING (VITERBI ALGORITHM)
  
  COPYRIGHT NOTICE: This computer program is free for non-commercial purposes.
  You may implement this program for any non-commercial application. You may 
  also implement this program for commercial purposes, provided that you
  obtain my written permission. Any modification of this program is covered
  by this copyright.

     Copyright 1994. Robert Morelos-Zaragoza. All rights reserved
  
     Portions of this program were written by Dojun Rhee, University of
     Hawaii, 1990.
   
     See also:
  
     R.H. Morelos-Zaragoza and S. Lin, ``QPSK Block Modulation Codes for
     Unequal Error Protection,'' IEEE Transactions on Information Theory,
     Vol. 41, No. 2, pp. 576-581, March 1995. 
  
   ===========================================================================
  
   For this block modulation code, a trellis diagram is as follows:
  
                        0     0     0     0     0
                     @-----@-----@-----@-----@-----@
                    / \   / \   / \   / \   / \   / \
                   /  2\ /  2\ /  2\ /  2\ /  2\ /   \
                 0/     x     x     x     x     x     \0
                 /    2/ \  2/ \  2/ \  2/ \  2/ \     \
                /     /   \ /   \ /   \ /   \ /   \     \
               / ----@-----@-----@-----@-----@-----@---- \
              / /2      0     0     0     0     0      2\ \
             @-|                                         |-@
              \ \1      1     1     1     1     1      1/ /
               \ ----@-----@-----@-----@-----@-----@---- /
                \     \   / \   / \   / \   / \   /     /
                 \    3\ /  3\ /  3\ /  3\ /  3\ /     /
                 3\     x     x     x     x     x     /3
                   \  3/ \  3/ \  3/ \  3/ \  3/ \   /
                    \ /   \ /   \ /   \ /   \ /   \ /
                     @-----@-----@-----@-----@-----@
                        1     1     1     1     1
  
   The structure of this diagram is read into the program from a file 
   called 'data1' which shall consist of the following 11 lines:
  
   0 0 0 0 1
   0 0 1 2 3
   0 7 7 7 7
   1 0 0 0 1
   1 0 1 2 3
   1 1 0 2 3
   1 1 1 0 1
   1 7 7 7 7
   2 0 0 0 1
   2 1 0 2 3
   2 7 7 7 7
  
   Each row in the file is of the form:
  
   +----------------------------------------------------------+
   | intime | initial | final | output1[0][0] | output1[0][1] |
   +----------------------------------------------------------+
  
   intime:        Type of section, initial (0), intermediate (1) and final (2)
   initial:       Initial state of a branch
   final:         Final state of a branch
   output1[0][0]: Branch label of upper subtrellis (coset 0)
   output1[0][1]: Branch label of lower subtrellis (coset 1)
  
   For example, the first row in the file indicates that the initial section
   has a branch, from the "0" state to the "0" state, with label "0" for the
   upper subtrellis and label "1" for the lower subtrellis. 
  
   To indicate the end of a section, the program checks whether the initial 
   state (second entry in a row) is equal to 7 (This number can be changed
   to suit).
 */

# include <math.h>
# include <stdio.h>
# include <limits.h>

/* 	#define MAXRAND 0x7fffffff              /* */
#define MAXRAND LONG_MAX                    /* for random number generation */
#define NUMSIM 5000                         /* no. simulations per 200 blocks */
#define DEBUG2                              /* For debugging only */

main()
  {
/* variables to index cosets of subtrellis */
   struct out1 { int coset[2][2];
               } ;
   int cosetvalue; 

/* variables for a trellis section */
   struct  trellis2 {                       
                    int time;               /* section type index */
                    int initial_state[5];   /* initial states */
                    int count;              /* number of branches in section */
                    int final_state[5];     /* final states */
                    struct out1 outbit1[5]; /* branch labels */
                  } ;

/* array for the trellis diagram */
   struct trellis2 trellis3[9];

/* miscelaneous variables */
   int test;
   int ifinal;
   int j;
   int i,i1,iii,ii;
   int jj,kk;
   int k;
   int l;
   int pathnum;                             /* branch counter in a section */
   int intime;                              /* section type */
   int initial;                             /* initial state */
   int final;                               /* final state */
   int output1[2][2];                       /* branch labels */
   int m;
   int n;
   int fflush();

/* variables for surviving paths */
   struct surv1 {
                double surv_metric[2][2];
                int surv_ini[2][2];
                int surv_out[2][2];
                int est[2][2][2];
               }; 

   struct sur_final1 {
                     double metric[2];
                     int ini[2];
                     int final[2];
                     int est[2][2];
                    }; 

   struct surv1 suv1[8];
   struct surv1 suv[8];
   struct sur_final1 surv_final1[8]; 
   struct sur_final1 surv_final[8];
   int suv_count1[2];
   int suv_count[2];
   int min_state;

/* variables for distance metrics */
   struct dis1 {
               double distance[2];
               double final;
               double f[2];
               int c12est[3];
              };
   struct metric1 {
                  struct dis1 branch[5];
                 };
   struct metric1 trellis_metric1[8];
   struct metric1 trellis_metric[8];
   double final0,final1;

/* variables for signal points */
   double q[2][4];

/* variables for estimates */
   struct est {
                int in[2];
                int out;
                int coset[8];
                int ini[3]; 
                int est[3];
              }; 
   struct est x_estimate[201];
   struct est x_estimate1[201];
   struct est y_estimate[201];
   long nn;

/* variables for block encoding and signal mapping */
   struct data {
                 int first;
                 int second[7];
               };             
   struct data code1[201];

/* variables for transmitted signal */
   int x[7];
   double rx[201][7];
   double ry[201][7];
   double pi,amp; 

/* variables for random numbers */
   double u1,u2,randmum;
   double y1,y2,v1,v2,s,w1,w2,x1,x2;
   int seed;

/* variables for bit error rate computations */
   double snr;
   double error = 0.0, error1 = 0.0, error2 = 0.0;
   double first_error, second_error;
   double biterror, biterr1, biterr2;

/* trellis data file pointer */
   FILE *fp1; 

/* =======================  READ TRELLIS STRUCTURE ======================== */
   seed = 1;
   printf("Enter amplitude, amp ? ");
   scanf("%lf", &amp);
   printf("amp = %lf\n", amp);
   error = 0;
   fp1=fopen("data1","r");
   for(i=1;i<=3;i++)
      {
      fscanf(fp1, "%d %d %d %d %d", &intime,&initial,&final,
						&output1[0][0],&output1[0][1]);
      pathnum = 1;  

/* read trellis diagram data */
      do{ 
            trellis3[i].time = intime;
            trellis3[i].initial_state[pathnum] = initial;
            trellis3[i].final_state[pathnum] = final;
            trellis3[i].outbit1[pathnum].coset[0][0] = output1[0][0];
            trellis3[i].outbit1[pathnum].coset[0][1] = output1[0][1];
            fscanf(fp1, "%d %d %d %d %d", &intime,&initial,&final,
						&output1[0][0],&output1[0][1]);
            pathnum++;
            } while (initial !=  7); 
         trellis3[i].count = pathnum - 1;
       }
   fclose(fp1);

/* construct last section of the trellis diagram */
   for(i=7;i<=7;i++)
      {
      for(j=1;j<=trellis3[3].count;j++)
         {
        trellis3[i].time = trellis3[3].time;
        trellis3[i].initial_state[j] = trellis3[3].initial_state[j];
        trellis3[i].final_state[j] = trellis3[3].final_state[j];
        trellis3[i].outbit1[j].coset[0][0] = trellis3[3].outbit1[j].coset[0][0];
        trellis3[i].outbit1[j].coset[0][1] = trellis3[3].outbit1[j].coset[0][1];
         }
         trellis3[i].count = trellis3[3].count;
       }

/* construct intermediate sections of trellis diagram */
   for(i=2;i<=6;i++)
      {
      for(j=1;j<=trellis3[2].count;j++)
         {
        trellis3[i].time = trellis3[2].time;
        trellis3[i].initial_state[j] = trellis3[2].initial_state[j];
        trellis3[i].final_state[j] = trellis3[2].final_state[j];
        trellis3[i].outbit1[j].coset[0][0] = trellis3[2].outbit1[j].coset[0][0];
        trellis3[i].outbit1[j].coset[0][1] = trellis3[2].outbit1[j].coset[0][1];
         }
         trellis3[i].count = trellis3[2].count;
       }

/* in-phase and quadrature components of the QPSK signal points */
   q[0][0] = amp;	q[1][0] = 0.0;
   q[0][1] = 0.0;	q[1][1] = amp;
   q[0][2] = -amp;	q[1][2] = 0.0;
   q[0][3] = 0.0;	q[1][3] = -amp;




/* ===========  MAIN LOOP OF SIMULATION (per 200 blocks)  ============ */

   for(ifinal=0; ifinal<NUMSIM; ifinal++)
   {
     seed++;  
     pi = 3.141592654;
     srandom(seed); 

     first_error = second_error = 0.0;

/* ===========  SECONDARY LOOP OF SIMULATION (200 blocks)  ============ */
     for(i=1;i<=200;i++)
     { 




/***************************************************************************
 *                                                                         *
 *                            ENCODING                                     *
 *                                                                         *
 ***************************************************************************/

/* most important information bit:   (7,1,7) code */
 
   code1[i].first = (random()&1024) >> 10 ; 

/* least important information bits: (7,6,2) code */
 
   code1[i].second[0] = (random()&1024) >> 10 ;
   code1[i].second[1] = (random()&1024) >> 10 ;
   code1[i].second[2] = (random()&1024) >> 10 ;
   code1[i].second[3] = (random()&1024) >> 10 ;
   code1[i].second[4] = (random()&1024) >> 10 ;
   code1[i].second[5] = (random()&1024) >> 10 ;
   code1[i].second[6] = code1[i].second[0] + code1[i].second[1] 
                        +code1[i].second[2] +code1[i].second[3] 
                        +code1[i].second[4] +code1[i].second[5]; 
   code1[i].second[6] %= 2;

/* Transmitted sequence, ontained by Gray mapping between two-bit symbols */
/* and QPSK signal points */
 
   x[0] = code1[i].first + 2*code1[i].second[0];
   x[1] = code1[i].first + 2*code1[i].second[1];
   x[2] = code1[i].first + 2*code1[i].second[2];
   x[3] = code1[i].first + 2*code1[i].second[3];
   x[4] = code1[i].first + 2*code1[i].second[4];
   x[5] = code1[i].first + 2*code1[i].second[5];
   x[6] = code1[i].first + 2*code1[i].second[6];





/***************************************************************************
 *                                                                         *
 *                         AWGN CHANNEL                                    *
 *                                                                         *
 ***************************************************************************/
 
   for(j=0;j<=6;j++)
     {
      /* Generate pair of iid Gaussian rv's */
      do
      {
       randmum = (double)(random())/MAXRAND;
       u1 = randmum*2 - 1.0;
       randmum = (double)(random())/MAXRAND;
       u2 = randmum*2 - 1.0;
       s = u1*u1 + u2*u2;