ldpc_2.cpp

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

     if (i == rightbound_index)	continue;

     Aux *= GetEdge(i).LeftBoundMessage;
  }

  Aux.Normalize();     
  
  return Aux;
}




void variable_node::CalcFinalMessage()
{
   FinalEstimate = CalcRightboundMessage( /* NO rightbound index */ -1 );
   Symbol = FinalEstimate.Decision();
}



void variable_node::CalcAllRightboundMessages()
{
   for (int i = 0; i < GetDegree(); i++)
      GetEdge(i).RightBoundMessage = CalcRightboundMessage( /* rightbound index */ i);

   CalcFinalMessage();
}




message &GenerateChannelMessage(GFq v,
                                channel &TransmitChannel, 
                                mapping &MapInUse,
                                double ChannelOut )
{
   static message InitialMessage;
   int q = MapInUse.GetQ();
   double CandidateIn, ChannelIn;
  
   ChannelIn = MapInUse.map(v);        // mapping of (0 + v) % q;
  
   // Generate InitialMessage
   InitialMessage.Set_q(q);
   for (GFq i(0); i.val < q; i.val++)
   {
      CandidateIn  = MapInUse.map(i + v);
      InitialMessage[i.val] = TransmitChannel.CalcProbForInput(
		                       ChannelOut, CandidateIn );
   }

   InitialMessage.Normalize();    // Make valid probability vector 

   return InitialMessage;
}






void variable_node::Initialize( channel &TransmitChannel, double ChannelOut )
{
   InitialMessage = GenerateChannelMessage(v, TransmitChannel, *MapInUse, ChannelOut );
   FinalEstimate = InitialMessage;
  
   // Generate Rightbound messages
   for (int j = 0; j < degree; j++)
      GetEdge(j).RightBoundMessage = InitialMessage;

}





BOOLEAN variable_node::IsRightConnectedTo( node *n )
{
  for (int i = 0; i < degree; i++)
    if (&edges[i]->RightNode() == n) return TRUE;

  return FALSE;
}

/**********************************************************************************
 *
 * Graph
 *
 **********************************************************************************/




void bipartite_graph::Reset( 
	     int p_N,                            // number of variable nodes
	     int lambda_degs[],                  // left-degrees, terminated by -1
	     double lambda_wts[],                // weights corresponding to above degrees
	     int rho_degs[],       		  // right-degree, terminated by NULL
	     double rho_wts[],                  // weights corresponding to above degrees
		 mapping &MapInUse)
{
  //-----------------------------------------------------------------------
  // Clean up
  //-----------------------------------------------------------------------
  Clear();

  //-----------------------------------------------------------------------
  // Calc M, N and E
  //-----------------------------------------------------------------------
  double Ratio;     // M/N
  double Nominator, Denominator;

  Nominator = 0;
  for (int i = 0; rho_degs[i] != -1; i++)
    Nominator += rho_wts[i] / rho_degs[i];

  Denominator = 0;
  for (int j = 0; lambda_degs[j] != -1; j++)
    Denominator += lambda_wts[j] / lambda_degs[j];

  Ratio = Nominator / Denominator;

  N = p_N;
  M = (int)ceil(Ratio * N);
  E = (int)floor(N / Denominator);
  
  //----------------------------------------------------------------------------------------
  // Create variables, checks and edges
  //----------------------------------------------------------------------------------------
  edges = new edge[E];
  variable_nodes = new variable_node[N];
  check_nodes = new check_node[M];
  EdgeStack = new edge* [E * 2];   // Buffer used to manage memory allocation for nodes

  if ((edges == NULL) ||
      (variable_nodes == NULL) ||
      (check_nodes == NULL) ||
      (EdgeStack == NULL))
    {
      cout << "Error allocating memory\n";
      exit(1);
    }

  //--------------------------------------------------------------------------
  // Generate sockets
  //--------------------------------------------------------------------------
  variable_node **left_sockets = new variable_node*[E];
  check_node **right_sockets = new check_node*[E];
  edge **EdgeStackPointer = EdgeStack;         // Auxiliary pointer to stack
  int node_index, socket_index, count_nodes_of_degree;
  
  node_index = 0;
  socket_index = 0;
  for (int i = 0; lambda_degs[i] != -1; i++)                // Loop through all left-degrees
    {
      int count_nodes_of_degree = (int)floor(lambda_wts[i] * E / lambda_degs[i]); // No. nodes of required degree
      for (int j = 0; j < count_nodes_of_degree; j++)    // Number of nodes with required left-degree
	{
	  for (int k = 0; k < lambda_degs[i]; k++)              // Number of sockets for each degree
	    left_sockets[socket_index++] = &variable_nodes[node_index];

	  variable_nodes[node_index].AllocateEdges(EdgeStackPointer, lambda_degs[i]);
	  variable_nodes[node_index].SetMapInUse(MapInUse);
	  variable_nodes[node_index].SetID(node_index);

	  node_index++;
	}
    }

  // Record no. of left sockets, may be different than E
  // due to the fact that lambda * E may not be whole numbers
  int count_left_sockets = socket_index;      

  // Modify N for same reason
  N = node_index;

  // Generate right sockets
  node_index = 0;
  socket_index = 0;
  for (int i = 0; rho_degs[i] != -1; i++)                  // Loop through all right-degrees
    {
      int CurrentDegree  = rho_degs[i];
      count_nodes_of_degree = (int)floor(rho_wts[i] * E / 
					 CurrentDegree );  // No of nodes of required degree
      for (int j = 0; j < count_nodes_of_degree; j++)         
	{
	  for (int k = 0; k < CurrentDegree; k++)      // Number of sockets for each degree
	    right_sockets[socket_index++] = &check_nodes[node_index];

	  check_nodes[node_index].AllocateEdges( EdgeStackPointer, CurrentDegree );
	  check_nodes[node_index].SetID(node_index);
	  
	  node_index++;
	}
    }
  // Record no. of right sockets, may be different than E
  // due to the fact that lambda * E may not be whole numbers
  int count_right_sockets = socket_index;      

  // Modify M,E for same reason
  M = node_index;
  E = min(count_left_sockets, count_right_sockets);

  //----------------------------------------------------------------------------------------
  // Generate permutations
  //----------------------------------------------------------------------------------------
  srand(time(NULL)); // Init random seed so that each call to function returns different set of values

  cout << "Starting bipartite graph...";

  int left_index, right_index;
  for (int i = 0; i < E; i++)
    {
      // Randomly select socket from first E - left_index (last left_index sockets represent
      // sockets that have already been selected)
      //      cout << " i = " << i << " E = " << E;
      int attempts = 0;

      do
	{
	  // It is important to select left_index randomly and not only right_index,
	  // because there is a significance to the order of bits within a code, sometimes with 
	  // first or last bits better protected.  If left_index is not selected randomly,
	  // the result would be a random tree, but in which lower degree left-nodes are of
	  // lower index within each constituent code, contrary to complete randomness.
	  left_index = uniform_random(E - i);
	  //left_index = i;
	  right_index = uniform_random(E - i);
	  if ((attempts++) > 100)
	    {
	      cout << "Warning: cycles\n";
	      break;
	    }
	}
    while (left_sockets[left_index]->IsRightConnectedTo(right_sockets[right_index]));

      if (right_index >= (E-i))
	{
	  cout << "right index or left_index exceed range\n"
	       << "right_index = " << right_index
	       << " left_index = " << left_index
               << " E = " << E 
               << " left_index = " << left_index << "\n";
	  exit(1);
	}

      edges[i].set_nodes(left_sockets[left_index], right_sockets[right_index]);

      // Overwrite current sockets with last sockets, so that they are not selected again
      right_sockets[right_index] = right_sockets[E - i - 1];
      left_sockets[left_index] = left_sockets[E - i - 1];
    }

  cout << "Done\n";

  // Clean-up
  delete left_sockets;
  delete right_sockets;
}




/************************************************************************
 *
 * LDPC code
 *
 ************************************************************************/



void LDPC_Code::Init_Messages( vector &ChannelOutput )
{
  if (ChannelOutput.GetSize() != Variables.GetLength())
  {
     cout << "LDPC_Code::Init_Messages: Incompatible sizes\n";
     exit(1);
  }

  for (int i = 0; i < Variables.GetLength(); i++)
     Variables[i].Initialize( *Channel, ChannelOutput[i] );
}





void LDPC_Code::GetZeroCodeword( vector &Codeword )
{
   Codeword.Allocate(Variables.GetLength());

   for (int i = 0; i < Variables.GetLength(); i++)
      Codeword[i] = Variables[i].GetZeroSignal();
}





void LDPC_Code::Leftbound_Iteration()
{
  for (int i = 0; i < Checks.GetLength(); i++)
    Checks[i].CalcAllLeftboundMessages();
}




void LDPC_Code::Rightbound_Iteration()
{
  for (int i = 0; i < Variables.GetLength(); i++)
     Variables[i].CalcAllRightboundMessages();
}




void LDPC_Code::FinalIteration()
{
  for (int i = 0; i < Variables.GetLength(); i++)
     Variables[i].CalcFinalMessage();
}




double LDPC_Code::Calc_Symbol_Error_Rate()
{
  double acc_correct = 0;

  for (int i = 0; i < Variables.GetLength(); i++)
     acc_correct += Variables[i].FinalEstimate.Decision().IsZero();

  return 1.0 - acc_correct / (double)Variables.GetLength();
}




double LDPC_Code::Calc_Rightbound_Symbol_Error_Rate()
{
  double acc_correct = 0;

  for (int i = 0; i < Graph.E; i++)
    acc_correct += Graph.edges[i].RightBoundMessage.Decision().IsZero();

  return 1.0 - acc_correct / Graph.E;
}







double LDPC_Code::Belief_Propagation_Decoder(int Count_Iterations)
{
  static char buffer[500];
  double Func_RC;
  double LastMin = INF;
  int CountIncreaseIterations = 0;

  cout << "Initial symbol error rate = " << Calc_Symbol_Error_Rate() << "\n";

  for (int i = 0; i < Count_Iterations; i++)
  {

      Leftbound_Iteration( );
      Rightbound_Iteration();

      Func_RC = Calc_Rightbound_Symbol_Error_Rate();
      sprintf(buffer, "%d: Rightbound SER = %1.10f, %s", i+1, Func_RC, CharTime());
      cout << buffer;

      // Stop when Func_RC doesn't fall for some consecutive iterations
      if (Func_RC < LastMin)
      {
         LastMin = Func_RC;
         CountIncreaseIterations = 0;
      }
      else
      {
         CountIncreaseIterations++;
         if (CountIncreaseIterations > 50)
            break;
      }

      if (Func_RC < 1e-7) break;
    }

   return Func_RC;
}