ldpc_2.h

关于Q元LDPC码的C++仿真程序

#ifndef LDPC_2
#define LDPC_2

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <iostream.h>
#include <fstream.h>
#include <ctype.h>
#include "Utils_1.h"
#include "Utils_2.h"
#include "Encoding.h"

typedef unsigned char BOOLEAN;
typedef unsigned char BYTE;

#define TRUE   1
#define FALSE  0


#define MAX_DEGREE	1000


/****************************************************************************
 *
 * Channel
 *
 ****************************************************************************/

class LDPC_Code;
class mapping;

// Channel types

class channel
{
public:
   virtual char *GetChannelName() = 0;

   // General functions ------------------------------------------
   GFq &MaxProbableForOutput( mapping &MapInUse );

   double ErrorUndefined(char *FuncName)
   {
      cout << GetChannelName() << "::" << FuncName << " undefined.\n";
      exit(1);
   }
   virtual void PrintChannelData(LDPC_Code &Code) = 0;
   virtual void ProcessMapping(LDPC_Code &Code) = 0;
   
   // Channel functions-------------------------------------------
   void SimulateOutputVector(vector &InVector, vector &OutVector);
   
   virtual double SimulateOutput(   double ChannelInput ) = 0;
   virtual double CalcProbForInput( double ChannelOutput, 
				    double ChannelInput ) = 0;

   // Statistical data -------------------------------------------
   virtual double CapacityInBits() = 0;
} ;


/****************************************************************************
 *
 * BSC Channel
 *
 ****************************************************************************/


class BSC_Channel : public channel
{
private:
   double channel_p;
   double source_p;
   double InterferenceProb;      // Dirty paper

public:
   // General functions ------------------------------------------

   BSC_Channel( )  : channel_p(0), InterferenceProb(0.5)
   { }
   void SetChannel_p( double p_channel_p )
   {
      channel_p = p_channel_p;
   }
   int Flip( double prob )
   {  return (my_rand() < prob); }
   virtual char *GetChannelName()
   { return "BSC_Channel"; }

   virtual void PrintChannelData(LDPC_Code &Code);
   virtual void ProcessMapping(LDPC_Code &Code);

   // Channel coding functions -----------------------------------
   virtual double SimulateOutput(   double ChannelInput )
   {  return (double)(int(ChannelInput) ^ Flip( channel_p ));  }
   virtual double CalcProbForInput( double ChannelOutput, double ChannelInput )
   {  
      if (ChannelOutput == ChannelInput)
         return 1.-channel_p;
      else
         return channel_p;
   }
   
   // Statistical functions -------------------------------------
   virtual double CapacityInBits();
   virtual double CalcProbForOutput( double ChannelOutput )
   {
      if (ChannelOutput == 1)
         return source_p;
      else if (ChannelOutput == 0)
         return 1 - source_p;
      else 
      {
         cout << "BSC_Channel.CalcProbForOutput: Invalid parameter\n";
         exit(1);
      }
   }
};




/****************************************************************************
 *
 * AWGN Channel
 *
 ****************************************************************************/

double GaussGenerate(double sigma);


class AWGN_Channel : public channel
{
private:
  double source_sigma;
  double noise_sigma;
public:
   // Specialized functions ----------------------------
   AWGN_Channel(double p_noise_sigma = -1) : 
     noise_sigma(p_noise_sigma) 
    { 
    }

   double NoiseVariance();
   double NoiseStdDev();
   void SetNoiseSigma(double p_noise_sigma) { noise_sigma = p_noise_sigma; }

   // General functions --------------------------------
   virtual char *GetChannelName()
   { return "AWGN_Channel"; }
   virtual void PrintChannelData(LDPC_Code &Code);
   virtual void ProcessMapping(LDPC_Code &Code);

   // Channel coding functions -------------------------
   virtual double SimulateOutput(double ChannelInput);
   virtual double CalcProbForInput(double ChannelOutput, double ChannelInput);

   // Statistical data ---------------------------------
   virtual double CapacityInBits();
} ;



/***************************************************************************
 *
 * mapping
 *
 ***************************************************************************/


class mapping
{
public:
   int q;
   double vals[MAX_Q];
public:
   mapping(int p_q = -1, double *p_vals = NULL) 
   {
        if (p_q != -1)
           Set_Q(p_q);

	if (p_vals != NULL)
	  bcopy( /*from*/ p_vals, /*to*/ vals, q * sizeof(double));
      }

   void Set_Q(int p_q)
   {
      if (p_q > MAX_Q)
        {
         cout << "Exceeded MAX_Q in mapping (simply increase)\n";
         exit(1);
        }

      q = p_q;
   }

   mapping(mapping &p_MapInUse) :
     q(p_MapInUse.q)
   {
      bcopy( /*from*/ p_MapInUse.vals, /*to*/ vals, q * sizeof(double));
   }

   BOOLEAN IsBinary()
   {
      if (q != 2)
         return FALSE;

      if (((vals[0] != 0) && (vals[0] != 1))   ||
          ((vals[1] != 0) && (vals[1] != 1)))
          return FALSE;
      
      return TRUE;
   }

   void GetFromFile(ifstream &file);

   double Average_E()
      {
         double sum_E = 0;
         for (int i = 0; i < q; i++)
            sum_E += pow(vals[i], 2);

         return sum_E / (double)q;
      }

   void Normalize()
      {
         // Normalize average energy to 1

         double factor = sqrt(Average_E());

         for (int i = 0; i < q; i++)
            vals[i] /= factor;
      }


   void operator/=(double d)
   {
      for (int i = 0; i < q; i++)
         vals[i] /= d;
   }

   void operator*=(double d)
   {
      for (int i = 0; i < q; i++)
         vals[i] *= d;
   }

   double map(int x)
     { return vals[x]; }

   double map(GFq &g)
     { return map(g.val); }
   
   int GetQ() {return q; }
} ;






inline ostream &operator<<( ostream &s, mapping &MapInUse )
{
  s << "[";
  for (int i = 0; i < MapInUse.GetQ(); i++)
  {
     if (i != 0) s << " ";
     s << MapInUse.vals[i];
  }
  
  s << "]";

  return s;
}


/***************************************************************************
 *
 * Messages
 *
 ***************************************************************************/

class message
{
public:
   signed int    q;
   double Probs[MAX_Q];
public:
   message(int p_q = -1) : q(p_q)
      { }

   message(message &M) :   q(M.q)
      { *this = M; }

   void HardMessage(GFq &g)
   {
      *this = 0;
      (*this)[g] = 1;
   }

   void Set_q(int p_q) 
      { q = p_q; }


   void DFT2();

   BOOLEAN operator<(message &m2)
   {
      message &m1 = *this;

      for (int i = 0; i <= m1.q; i++)
         if (m1[i] >= m2[i])
            return FALSE;

      return TRUE;
   }

   double Maximum()
   {
      double maximum = -INF;
      for (int i = 0; i < q; i++)
         if (Probs[i] > maximum)
            maximum = Probs[i];

      return maximum;
   }

   double &operator[](int i) { return Probs[i]; }
   double &operator[](GFq i) { return Probs[i.val]; }
   message &operator=(message &M)
      {
         if (q != M.q) Set_q(M.q);

         bcopy(/* from */ M.Probs, /* to */ Probs, sizeof(double)*q);
         return *this;
      }

   message &operator=(double d)
      {
         for (int i = 0; i < q; i++)
            Probs[i] = d;
         return *this;
      }

   message &operator*=(message &M)
      {
         for (int i = 0; i < q; i++)
            Probs[i] *= M.Probs[i];

         return *this;
      }

   message &operator*(message &M2)
   {
      static message Aux(q);
      message &M1 = *this;

      for (int i = 0; i < q; i++)
         Aux.Probs[i] = M1.Probs[i] * M2.Probs[i];

      return Aux;
   }

   message &operator*(double d)
      {
         static message Aux(q);

         for (int i = 0; i < q; i++)
            Aux[i] = Probs[i] * d;

         return Aux;
      }

   BOOLEAN operator==(message &m)
   {
      if (q != m.q)
         return FALSE;
      
      for (int i = 0; i < q; i++)
//         if (Probs[i] != m[i])
         if (fabs(Probs[i] - m[i]) > EPSILON)
            return FALSE;

      return TRUE;
   }

   BOOLEAN operator==(double d)
   {
      for (int i = 0; i < q; i++)
         //if (fabs(Probs[i] - d) > EPSILON)
         if (Probs[i] != d)
            return FALSE;

      return TRUE;
   }
   
   message &operator+=(message &M)
   {
      for (int i = 0; i < q; i++)
         Probs[i] += M.Probs[i];

      return *this;
   }
   
   message &operator+(message &M)
   {
      static message Aux(q);

      for (int i = 0; i < q; i++)
         Aux[i] = Probs[i] + M.Probs[i];

      return Aux;
   }

   message &operator/=(double d)
   {
      for (int i = 0; i < q; i++)
         Probs[i] /= d;

      return *this;
   }


   double sum()
   {
      double aux = 0;
      for (int i = 0; i < q; aux += Probs[i++]);
      return aux;
   }

   void Normalize()
   {
	   double aux = sum();

	   if (aux > 0)
	   {
	      for (int i = 0; i < q; i++)
		     Probs[i] /= aux;
	   }
	   else
	   {
		  for (int i = 0; i < q; i++)
		     Probs[i] = 1./(double)q;
	   }
   }

   void Clear()
   {
      bzero(Probs, sizeof(double)*q);
   }


   message &Convolve(message &M2)
   {
      message M1(*this);   // Auxiliary

      // Clear this
      Clear();

      for (GFq i(0); i.val < q; i.val++)
         for (GFq j(0); j.val < q; j.val++)
            Probs[i.val] += M1[j] * M2[i - j];

      return *this;
   }

   message &ApproxConvolve(message &M2);
   void Approximate();

   message &MaxConvolve(message &M2)
   {     // max-prod version
      message M1(*this);   // Auxiliary

      // Clear this
      Clear();
      double max;
      double candidate;

      for (GFq i(0); i.val < q; i.val++)
      {
         max = -1;
         for (GFq j(0); j.val < q; j.val++)
         {
            candidate = M1[j] * M2[i - j];
            if (candidate > max)
               max = candidate;
         }
         Probs[i.val] = max;
      }

      return *this;
   }



   void PermutePlus(GFq &g)
   {
      message Aux;

      Aux = *this;
      for (GFq i(0); i.val < q; i.val++)
         Probs[i.val] = Aux[i + g];
   }


   void PermuteTimes(GFq &g)
   {
      message Aux;

      Aux = *this;
      for (GFq i(0); i.val < q; i.val++)
         Probs[i.val] = Aux[i * g];
   }


   message &Reverse()
   {
      message Aux;

      Aux = *this;
      for (GFq i(0); i.val < q; i.val++)
         Probs[i.val] = Aux[i.Minus()];

      return *this;
   }


   GFq &Decision()
   {
     static GFq Candidate(0);
     double max = -1;
     int count_max = 0;

     for (GFq i(0); i.val < q; i.val++)
     {
       if (Probs[i.val] > max)
       {
          max = Probs[i.val];
          Candidate = i;
          count_max = 1;
       }
       else if (Probs[i.val] == max)
          count_max++;
     }

     if (count_max > 1)
     {     // If more than one maximum - randomly select
        int selection = uniform_random(count_max) + 1;      // Between 1 and count_max
        int found_so_far = 0;

        for (GFq i(0); i.val < q; i.val++)
           if (Probs[i.val] == max)
           {
              Candidate = i;
              found_so_far++;
              if (found_so_far == selection) break;
           }
     }

     return Candidate;
   }


   double ProbCorrect()
     {
       int count_zero = 1;

       for (int i = 1; i < q; i++)
       {
         if (Probs[i] > Probs[0])      
	 {
            return 0;
         }
         else if (Probs[i] == Probs[0])
         {
            count_zero++;
         }
       }

       return 1. / count_zero;
     }


   double Entropy()
   {
      double aux = 0;

      for (int i = 0; i < q; i++)
      {
         if (Probs[i] != 0)    
         {
            double aux2 = Probs[i] * log(Probs[i]);
            aux += - clip(aux2);
         }
      }

      return aux;
   }