ldpc_2.cpp

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

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <iostream.h>
#include <fstream.h>
#include <ctype.h>
#include "Utils_1.h"
#include "LDPC_1.h"
#include "LDPC_2.h"


/*********************************************************************************
 *
 * MESSAGE
 *
 *********************************************************************************/


void message::DFT2()          // A real-valued DFT - also IDFT
{
	static message Aux;
      GFq mask, n0_index, n1_index;
	BYTE j_bit;

      if (GFq::IsPrimeQ)
      {
         cout << "Error, message::DFT2:  Only q that is a power of 2 is supported.\n";
         exit(1);
      }

      for (int i  = 0; i < GFq::log_2_q; i++)
      {
	   Aux = *this;			// Initialize
	   mask.val = 1 << i;

         for (GFq j(0); j.val < q; j.val++)
         {
            j_bit = (j.val & mask.val) >> i;	// obtain value of j
		n0_index.val = j.val & (~mask.val);	// turn bit off
		n1_index.val = j.val | mask.val;    // turn bit on

		//cout << (int)j_bit;

            if (j_bit == 0)
               Probs[j.val] = Aux[n0_index] + Aux[n1_index];
            else
               Probs[j.val] = Aux[n0_index] - Aux[n1_index];
         }

      }
}




/*********************************************************************************
 *
 * CHANNEL
 *
 *********************************************************************************/


void channel::SimulateOutputVector(vector &InVector, vector &OutVector)
{
   OutVector.Allocate(InVector.GetSize());

   for (int i = 0; i < InVector.GetSize(); i++)   // Add noise to each component
      OutVector[i] = this->SimulateOutput( /* channel in */ InVector[i] );    
}


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


void BSC_Channel::ProcessMapping(LDPC_Code &Code)
{
   mapping &MapInUse = Code.MapInUse;

   for (GFq i(0); i.val < MapInUse.q; i.val++)
      if ((MapInUse.map(i) != 0) && (MapInUse.map(i) != 1))
      {
         cout << "mapping levels must be (0,1)\n";
         exit(1);
      }
}




double H2(double p)
{
   return -p*log2(p) - (1-p)*log2(1-p);
}




double H2Reverse(double H)
{
   double pmin = 0, 
          pmax = 0.5, 
          tol = 0.0000001,
          pmid;

   if (H == 0)
      return 0;
   else if (H == 1.)
      return 0.5;

   while ((pmax - pmin) > tol)
   {
      pmid = (pmin + pmax) / 2.;
      if (H2(pmid) > H)
         pmax = pmid;
      else
         pmin = pmid;
   }

   return pmid;
}





double BSC_Channel::CapacityInBits()
{
   return 1 - H2(channel_p);
}




void BSC_Channel::PrintChannelData(LDPC_Code &Code)
{
   double BitRate;

   BitRate = Code.Calc_Bit_Rate();
   cout  << "channel_p = " << channel_p
         << "\nCapacity (bits per channel use) = " << 1 - H2(channel_p) 
         << " Max channel_p for Bit rate: " << H2Reverse(1. - BitRate);
}


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


double GaussGenerate(double sigma)
  // Simulate the result of passing the zero vector through the AWGN
{
  static const double pi = 3.141592653;
  double normal_random_number, x1, x2;

  x1 = my_rand();
  x2 = my_rand();

  if (x1 < EPSILON)
   x1 = EPSILON;

  normal_random_number = sqrt(-2 * log(x1)) * cos(2*pi * x2);
  
  clip(normal_random_number);
  
  return sigma * normal_random_number;
}



void AWGN_Channel::ProcessMapping(LDPC_Code &Code)
{
   Code.MapInUse.Normalize();
}

  



double AWGN_Channel::CapacityInBits()
{
   double No = pow(noise_sigma,2.);
   double SNR = 1./No;
   return 0.5 * log(1. + SNR) / log(2.);
}



void AWGN_Channel::PrintChannelData(LDPC_Code &Code)
{
   double BitRate, No, SNR, SNR_dB;
   
   BitRate = Code.Calc_Bit_Rate();
   No = pow(noise_sigma,2.);
   SNR = 1./No;
   SNR_dB = 10. * log10(SNR);

   cout << "SNR(dB) = " << SNR_dB
        << " SNR = " << SNR
        << " Noise Sigma = " << noise_sigma
        << "\nCapacity at SNR (symbols per channel use) = " 
        << 0.5 * log(1. + SNR) / log((double)GFq::q)
        << "\nCapacity at SNR (bits per channel use) = " 
        << 0.5 * log(1. + SNR) / log(2.)
        << "\nMinimum SNR for rate (dB) = " 
        << 10. * log(pow(2., 2.*BitRate) - 1.) / log(10.)
        << " (absolute value) = " << pow(2., 2.*BitRate) - 1;
}



// override virtual functions
double AWGN_Channel::CalcProbForInput( double ChannelOutput, 
				       double ChannelInput )
{
   static const double sqrt_2pi = sqrt(2*3.141592653);

   double noise_prob = (1/(sqrt_2pi * noise_sigma)*
	    exp(-pow(ChannelOutput-ChannelInput,2.)/(2.*NoiseVariance())));

   return noise_prob;
}



double AWGN_Channel::SimulateOutput(double ChannelInput)
  // Simulate the result of passing the zero vector through the AWGN
{
  return ChannelInput + GaussGenerate(noise_sigma);
}



double AWGN_Channel::NoiseVariance() 
{ 
  return pow(noise_sigma, 2); 
}




double AWGN_Channel::NoiseStdDev() 
{ 
  return noise_sigma; 
}








/****************************************************************************
 *
 * Check node
 *
 ****************************************************************************/



GFq &check_node::Value()
{
   static GFq Aux;

   Aux = 0;
   for (int i = 0; i < GetDegree(); i++)
      Aux += GetEdge(i).label * GetEdge(i).LeftNode().Symbol;

   return Aux;
}




node &check_node::AdjacentNode(int index)
{
   return GetEdge(index).LeftNode();
}



message &FastCalcLeftboundMessage(message AuxLeft[], 
                                  message AuxRight[],
			          int left_index,
                                  int degree)
{
  static message Aux;

  // Init to delta pulse 
  Aux.Set_q(GFq::q);

  Aux = AuxLeft[left_index];
  Aux *= AuxRight[left_index];
  Aux.DFT2();

  Aux.Reverse();	// estimate of symbol value = (- sum of others)

  return Aux;
}



void check_node::CalcAllLeftboundMessages( )
{
  static message Vectors[MAX_DEGREE];
  static message AuxLeft[MAX_DEGREE];
  static message AuxRight[MAX_DEGREE];
  static GFq ZERO(0);

  // Calc auxiliary DFTs
  for (int i = 0; i < GetDegree(); i++)
  {
     Vectors[i] = GetEdge(i).RightBoundMessage;
     Vectors[i].PermuteTimes(GetEdge(i).label.Inverse());
  }

  //-------------------------------------------------------------
  // If power of two - use DFT2
  //-------------------------------------------------------------
  if (!GFq::IsPrimeQ)
  {
     for (int i = 0; i < GetDegree(); i++)
	 {
         Vectors[i].DFT2();
	 }

      // Calc auxiliary values
     AuxLeft[0].Set_q(GFq::q);
     AuxLeft[0] = 1.;
     for (int i = 1; i < GetDegree(); i++)
     {
        AuxLeft[i] = AuxLeft[i - 1];
        AuxLeft[i] *= Vectors[i - 1];
     }

     AuxRight[GetDegree() - 1].Set_q(GFq::q);
     AuxRight[GetDegree() - 1] = 1.;
     for (int i = GetDegree() - 2; i >= 0; i--)
     {
        AuxRight[i] = AuxRight[i + 1];
        AuxRight[i] *= Vectors[i + 1];
     }

     // Calc leftbound messages
     for (int i = 0; i < GetDegree(); i++)
     {
       GetEdge(i).LeftBoundMessage = FastCalcLeftboundMessage(AuxLeft, AuxRight, i, GetDegree());
       GetEdge(i).LeftBoundMessage.PermuteTimes(GetEdge(i).label);
     }
  }
}



/*************************************************************************
 *
 * Node
 *
 *************************************************************************/


void node::Disconnect()
{
   // while edges remain, disconnect the first one 
   while (degree > 0) edges[0]->Disconnect();
}



/*************************************************************************
 *
 * Variable node
 *
 *************************************************************************/


node &variable_node::AdjacentNode(int index)
{
   return GetEdge(index).RightNode();
}




message &variable_node::CalcRightboundMessage( int rightbound_index )
// rightbound_index is -1 if calculation is meant for final estimate
{
  static message Aux;

  Aux = InitialMessage;

  for (int i = 0; i < GetDegree(); i++)
  {
     // ignore intrinsic information