order0.c

error correction code

/* This function "order0_dec" written by Marc Fossorier performs
   order-0 decoding as described in his paper.     "Gg_t" is the
   two-dimensional integer array containing the generator matrix of
   the code. "R" is the length "N" double-valued received vector
   containing the received real numbers. "out_D" is the result of
   order-0 decoding. */ 
// ------------------------------------------------------------------------
// This program is complementary material for the book:
//
// R.H. Morelos-Zaragoza, The Art of Error Correcting Coding, Wiley, 2002.
//
// ISBN 0471 49581 6
//
// This and other programs are available at http://the-art-of-ecc.com
//
// You may use this program for academic and personal purposes only. 
// If this program is used to perform simulations whose results are 
// published in a journal or book, please refer to the book above.
//
// The use of this program in a commercial product requires explicit
// written permission from the author. The author is not responsible or 
// liable for damage or loss that may be caused by the use of this program. 
//
// Copyright (c) 2002. Robert H. Morelos-Zaragoza. All rights reserved.
// ------------------------------------------------------------------------

# include "def.h"

void order0(Gg_t,R,out_D)
    int Gg_t[K][N];
    double R[N],out_D[N];
{
 
    double abs_R[N],C[N][2];
    int Dd[N],DD[N],GG[K][N],i,j,Gg[K][N];
    int permutation_final[N],permutation_I[K],permutation_R[N];
    double zero[K],fabs();
    void peterson_I(),quick_sort_track();
    void switch_column_I();
    void switch_vector(),switch_matrix();
 
   for (i=0;i < K;i++)
	for (j=0;j < N;j++)
	    Gg[i][j] = Gg_t[i][j]; 

   for (i=0;i < N;i++)
   {
       abs_R[i] = (double) fabs((double) R[i]);
       if (i < K)
       {  
           permutation_I[i] = i;
       }  
       permutation_R[i] = i;
   }

   quick_sort_track(abs_R,permutation_R,0,N-1);

/* Note that the permutations returned by qs are from small to big ! */

/* Move the K independent most informative symbols to the front */
/* Column operations => Changes the code into an equivalent code */
 
   switch_vector(R,permutation_R,N-1);
 
   switch_matrix(Gg,permutation_R,K-1,N-1);
 
/* Put into systematic code */
/* Row operations => keeps the code */
 
   peterson_I(Gg,zero);
                           /* Obtain the K columns of I */
 
   quick_sort_track(zero,permutation_I,0,K-1);

/* Switch the rows to obtain the K I_rows in the right order
   => keeps the code */
 
   for (i=0; i<K; i++)
   {
       for (j=0; j<N; j++)
       { 
           GG[i][j] = Gg[permutation_I[i]][j];
       }
   }
 
 
/* Switch the column to obtain the K I_column in a consecutive order
   => changes the code */
 
   switch_column_I(GG,Gg,R,permutation_R,permutation_final);

/* Compute 2N costs and K most probable bits */
 
   for (i=0; i<N; i++)
   {
       C[i][0] = (R[i] - 1.0) * (R[i] - 1.0);
       C[i][1] = (R[i] + 1.0) * (R[i] + 1.0);
       if (i < K)
       { 
           if (C[i][0] < C[i][1])
           {
              Dd[i] = 0;
           }
           else
           {
              Dd[i] = 1;
           }
       }
       else
       {
          Dd[i] = 0;
          for (j=0; j<K; j++)
          {
             Dd[i] = Dd[i] + Dd[j] * Gg[j][i];
          }
          Dd[i] = Dd[i] % 2;
       }
    }
 
   for (i=0; i<N; i++)
   {
       DD[permutation_final[i]] = Dd[i];
   }
 
   /*for (i=0; i<O_OUT_iovec_len; i++)*/ 
   for (i=0; i < N; i++)   /* Using a systematic code */
   {
      if (DD[i] == 0)
         out_D[i] = - 1.0;
      else
         out_D[i] = 1.0;
   }

}


void switch_column_I(GG,Gg,R,permutation_R,permutation_final)
  /* copy GG into Gg */   
   int GG[K][N],Gg[K][N];
   double R[];
   int permutation_R[],permutation_final[];
{
    int i,j,index,count,start;
    int record[N],record_I[N];
    double temp[N];
 
    index = 0;
    count = 0;
 
    for (i=0; i < K+index; i++)   /* record the column positions to switch */
    {
         if (GG[i-index][i] == 0)
         {
              record[index] = i;
              index++;
         }
         else
         {
              record_I[count] = i;
              count ++;
         }
    }
    start = K + index - 1;
    i = N-1;
    while (i > start)            /* Unchanged part */
    {
        for(j=0; j<K; j++)
        {
            Gg[j][i] = GG[j][i];
            permutation_final[i] = permutation_R[N-1-i];
        }
        i--;
    }    
 

    while (index > 0)            /* Copy from R to L dependent columns found
                                    in first positions */
    {
        temp[i] = R[i];
        for(j=0; j<K; j++)
        {
            Gg[j][i] = GG[j][record[index-1]];
            permutation_final[i] = permutation_R[N-1-record[index-1]];
        }
        if (i < record[index-1])
        {
            R[i] = temp[record[index-1]];
        }
        else
        {
            R[i] = R[record[index-1]];
        }
        i--;
        index--;
    }
    while (count > 0)    /* Copy I_columns */
    {
        temp[i] = R[i];
        for(j=0; j<K; j++)
        {
            Gg[j][i] = GG[j][record_I[count-1]];
            permutation_final[i] = permutation_R[N-1-record_I[count-1]];
        }
 
        if (i < record_I[count-1])
        {
            R[i] = temp[record_I[count-1]];
        }
        else
        {
            R[i] = R[record_I[count-1]];
        } 
        i--;
        count--;
    }
}

void switch_vector(R,permutation,length)
     double R[];
     int permutation[];
     int length;
{
  int i;
  double temp[N];
  
  for (i=0;i<=length;i++)
    {
      temp[i] = R[i];
      if (permutation[length-i] > i)   /* data not overwritten */
	{
	  R[i] = R[permutation[length-i]];
	}
      else
	{
	  R[i] = temp[permutation[length-i]];
	  /* data overwritten but stored in temp */
	}
    }
}

void switch_matrix(Gg,permutation,k,n)
     int Gg[K][N];
     int permutation[];
     int k,n;
{
  int i,j;
  double temp[K][N];
  
  for (j=0;j<=n;j++)
    {
      for (i=0; i<=k; i++)
	{
	  temp[i][j] = Gg[i][j];
	  if (permutation[n-j] > j)   /* column not overwritten */
	    {
	      Gg[i][j] = Gg[i][permutation[n-j]];
	    }
	  else
	    {
	      Gg[i][j] = temp[i][permutation[n-j]];
	      /* data overwritten but stored in temp */
	    }
	}
    }
}

void peterson_I(Gg,zero)
     int Gg[K][N];
     double zero[N];
{
  int i,j,l,m;
  
  for (i=0; i<K; i++)
    {
      j = 0;
      zero[i] = 0.0;
      while (Gg[i][j] == 0)
	{
	  j++;
	  zero[i] = zero[i] + 1.0;
	}
      for (l=0; l<K ; l++)
	{
          if ((l != i) && (Gg[l][j] == 1))
	    {
	      for (m=0; m<N; m++)
		{
                  Gg[l][m] = (Gg[l][m] + Gg [i][m])%2;
		}
	    }
	}
    }
}

/* Recursive quick_sort */
void quick_sort_track(vals,permutation,left,right) /* Tracks the permutations */        
     double vals[];
     int permutation[];
     int left,right;
{
  /* ascending = left wall, descending = right wall */
  
  int ascending = left-1, descending = right;
  double ref_val,temp;
  int temp2;
  
  /* select the comparison value */
  
  if (right > left)                 
    {
      ref_val = vals [right];
      
      for(;;)                 
        {
	  
	  /* while element smaller than reference value,
	     move left wall upwards
	     */
	  
	  while (vals [++ascending] < ref_val);
	  
	  /* while element larger than reference value,
	     move right wall downwards
	     */
	  
	  while (vals [--descending] > ref_val);
	  
	  if (ascending >= descending) break;    
	  
	  /* if the walls have not passed each other,
	     exchange values
	     (if not access to function need to save index)
	     */
	  
	  temp=vals[ascending];
	  temp2=permutation[ascending];
	  vals[ascending]=vals[descending];
	  vals[descending]=temp;
	  permutation[ascending]=permutation[descending];
	  permutation[descending]=temp2;
	  
	}                                 
      
      
      temp=vals[ascending];
      temp2=permutation[ascending];
      vals[ascending]=vals[right];
      vals[right]=temp;
      permutation[ascending]=permutation[right];
      permutation[right]=temp2;
      
      
      /* if descending wall has not yet reached left partition,
	 call quick-sort again with this smaller array
	 */
      
      quick_sort_track( vals, permutation, left, ascending-1);
      
      /* if ascending wall has not yet reached right partition,   
	 call quick-sort again with this smaller array
	 */
      
      quick_sort_track( vals, permutation, ascending+1, right);
    }     
  
}