#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 functionsdouble 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