gf.c

Kasami序列产生器 用于扩频通信

/*
 * $Log: gf.c,v $
 * Revision 1.1  2000/05/03 14:30:04  bjc97r
 * Initial revision
 *
 */
char *_gf_id = "$Id: gf.c,v 1.1 2000/05/03 14:30:04 bjc97r Exp $";

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "gf.h"

/*
 * gf_table() calculates all the GF(2^M) elements and store them
 * in a table in ascending order. The element for a^n can be accessed
 * at (n+1)th index of the table, where a is the root of Poly.
 * The table index starts from 0.
 */
unsigned long * gf_table( unsigned M, unsigned long Poly ) 
{
  unsigned long *table;      /* the table containing the GF element */
  unsigned long  poly_mask;  /* the polynomial mask */
  unsigned long  abs_mask;   /* the absolute bit mask */
  unsigned long  i;          /* loop index */

  {  /* abs_mask has '1's in the last M bits. */
    unsigned long pos;
    for( abs_mask = 1, pos=1, i=0; i < M; i++, pos <<= 1 ) 
      abs_mask |= pos;
  }

  /* poly_mask represents the feedback part of Galois LFSR */
  poly_mask = Poly & abs_mask;

  /* Now, generate all the elements in the Galois Field. */
  {
    unsigned long size = 1 << M;  /* the size of the GF */
    unsigned long element;        /* the current element in the GF */
    unsigned long carry;          /* carry flag */
    unsigned long c_mask = size >> 1; /* mask for the most significant
					 bit of an elem. */

    table = (unsigned long*) malloc( size * sizeof(unsigned long) );
    element = 0; table[0] = element;
    element = 1; table[1] = element;
    for ( i=2; i < size; i++ ) {
      carry = element & c_mask;
      element <<= 1;
      element &= abs_mask;
      if ( carry ) element ^= poly_mask;
      table[i] = element;
    }
  }

  return table;
}

/*
 * gf_exp_table() calculates all the GF(2^M) elements and store only
 * a^(R*n), where a is the root of the polynomial 'Poly', R is the exponent,
 * and n is from 0 to N-1. 
 */
unsigned long * gf_exp_table( unsigned M, unsigned long Poly, 
			      unsigned R, unsigned N ) 
{
  unsigned long *table;     /* the table containing the GF element */
  unsigned long  poly_mask; /* the polynomial mask */
  unsigned long  abs_mask;  /* the absolute bit mask */
  unsigned long  n, r;      /* loop index */

  { /* abs_mask has '1's in the last r bits. */
    unsigned long pos; int i;
    for( abs_mask = 1, pos=1, i=0; i < M; i++, pos <<= 1 ) 
      abs_mask |= pos;
  }

  /* poly_mask represents the feedback part of Galois LFSR */
  poly_mask = Poly & abs_mask;

  /* Now, generate all the elements in the Galois Field. */
  {
    unsigned long element;       /* the current element in the GF */
    unsigned long carry;         /* carry flag */
    unsigned long c_mask = 1 << (M-1); /* mask for the most significant
					  bit of an elem. */

    table = (unsigned long*) malloc( N * sizeof(unsigned long) );
    element = 1; table[0] = element;
    for ( n=1; n < N; n++ ) {
      for ( r=0; r < R; r++ ) {
	carry = element & c_mask;
	element <<= 1;
	element &= abs_mask;
	if ( carry ) element ^= poly_mask;
      }
      table[n] = element;
    }
  }

  return table;
}
/*
 * gfl_exp_table() is a 64-bit version of gf_exp_table().
 */
u2long * gfl_exp_table( unsigned M, u2long Poly, unsigned R, unsigned N ) 
{
  u2long *table;       /* the table containing the GF element */
  u2long  poly_mask;   /* the polynomial mask */
  u2long  abs_mask;    /* the absolute bit mask */
  unsigned Mh, Ml;
  unsigned long  n, r; /* loop index */

  if ( M > 32 ) { Ml = 32; Mh = M - Ml; }
  else          { Ml =  M; Mh = 0; }

  { /* abs_mask has '1's in the last M bits. */
    unsigned long pos; int i;
    for( abs_mask.l=0, pos=1, i=0; i < Ml; i++, pos <<= 1 ) abs_mask.l |= pos;
    for( abs_mask.h=0, pos=1, i=0; i < Mh; i++, pos <<= 1 ) abs_mask.h |= pos;
  }

  /* poly_mask represents the feedback part of Galois LFSR */
  poly_mask.h = Poly.h & abs_mask.h;
  poly_mask.l = Poly.l & abs_mask.l;

  /* Now, generate all the elements in the Galois Field. */
  {
    u2long element; /* the current element in the GF */
    u2long c_mask;  /* mask for the most significant bit of an elem. */
    unsigned long  carry;  /* carry flag */
    unsigned long  overflow;
    unsigned long  ovf_mask=0;

    if ( Mh ) {  /* More than 32 bits */
      c_mask.h = 1 << ( Mh-1 );
      c_mask.l = 0;
      ovf_mask = 1 << ( sizeof(unsigned long)*8 - 1 );
    }
    else {  /* Less than or equal to 32 bits */
      c_mask.h = 0;
      c_mask.l = 1 << ( Ml-1 ); 
    }

    table = (u2long*) malloc( N * sizeof(u2long) );
    element.l = 1; element.h = 0;
    table[0]  = element;
    for ( n=1; n < N; n++ ) {
      for ( r=0; r < R; r++ ) {
	if ( Mh ) {
	  carry       = element.h & c_mask.h;
	  overflow    = element.l & ovf_mask;
	  element.l <<= 1;
	  element.h <<= 1;
	  if ( overflow ) element.h |= 1;
	  element.h  &= abs_mask.h;
	  if ( carry ) {
	    element.l ^= poly_mask.l;
	    element.h ^= poly_mask.h;
	  }
	}
	else {
	  carry       = element.l & c_mask.l;
	  element.l <<= 1;
	  element.l  &= abs_mask.l;
	  if ( carry ) element.l ^= poly_mask.l;
	}
      }
      table[n] = element;
    }
  }

  return table;
}

/*
 * Find the generator polynomial for the sub m-sequence based on 
 * the main polynomial to generate a set of Kasami sequences.
 * It returns the generator polynomail for the sub m-sequence.
 */
unsigned long find_subpoly( unsigned L, unsigned long P )
{
  unsigned long l1;  /* the order of the sub m-sequence */
  unsigned long p1;  /* the polynomial for the sub m-sequence */
  unsigned long r;   /* the decimation exponent */
  unsigned long *gf; /* GF(2^L) element table for a^(r*n),
			where a is the root. */
  unsigned long sum0;/* a^(r*l1)+1, where a is the root of P */
  unsigned long sum; /* p1(a^r) */
  unsigned long np;  /* the number of subpoly candidates */
  unsigned long mask;/* a bit mask */

  int i, j; /* loop indices */

  if ( L & 1 ) {
    fprintf( stderr, "find_subpoly(): L=%d is odd!!\n", L);
    return 0;
  }

  l1 = L >> 1;          /* the order of the subpoly */
  r  = ( 1 << l1 ) + 1; /* decimation exponent */
  gf = gf_exp_table( L, P, r, l1+1 );
  p1 = ( 1 << l1 ) + 3; /* a candidate for the subpoly */
  sum0 = gf[l1] ^ 1;    /* a^(r*l1)+1 */
  np = 1 << (l1-1);     /* ..because the last digit of poly is always 1. */

  for ( i=1; i < np; i++, p1 += 2 ) {
    sum  = sum0; /* initial value */
    mask = 2;
    for ( j=1; j < l1; j++, mask <<= 1 ) {
      if ( p1 & mask ) 
	sum ^= gf[j];
    }
    if ( sum == 0 )
      break;
  }
  free(gf);

  if ( i == np ) {
    fprintf( stderr, "find_subpoly(): no subpoly found!\n");    
    return 0L;
  }
  
  return p1;
}

/*
 * This is 64bit version of find_subpoly() function.
 * It can handle polynomials with high orders uptp 62.
 * The polynomial is represented as a octal string
 * without leading 0 in 'P'. 'L' is the order of the polynomail 'P'.
 * It returns the generator polynomail for the sub m-sequence.
 */
unsigned long find_longsubpoly( unsigned L, char *P )
{
  unsigned long l1;    /* the order of the sub m-sequence */
  unsigned long p1;    /* the polynomial for the sub m-sequence */
  unsigned long r;     /* the decimation exponent */
  u2long *gf;  /* GF(2^L) element table for a^(r*n), where a is the root. */
  u2long sum0; /* a^(r*l1)+1, where a is the root of P */
  u2long sum;  /* p1(a^r) */
  unsigned long np;    /* the number of subpoly candidates */
  unsigned long mask;  /* a bit mask */
  u2long poly;

  int i, j;  /* loop indices */
  

  /*
   * Bounday check first.
   */
  if ( L & 1 ) {
    fprintf( stderr, "find_longsubpoly(): L(%d) is odd!!\n", L);
    return 0;
  }
  if ( L > 62 ) {
    fprintf(stderr, "find_longsubpoly: L(%d) is too big.\n", L );
    fprintf(stderr, "find_longsubpoly: ..cannot handle deg > 62\n");
    return 0;
  }
  if ( strlen(P) > 21 ) {
    fprintf(stderr, "find_longsubpoly: P(%s) is too long.\n", P);
    fprintf(stderr, "find_longsubpoly: P should be "
	    "an octal representation without leading 0.\n");
    return 0;
  }

  poly = strtou2long( P );

  l1 = L >> 1;           /* the order of the subpoly */
  r  = ( 1 << l1 ) + 1;  /* decimation exponent */
  gf = gfl_exp_table( L, poly, r, l1+1 );
  p1 = ( 1 << l1 ) + 3;  /* a candidate for the subpoly */
  sum0 = gf[l1]; sum0.l ^= 1;  /* a^(r*l1)+1 */
  np = 1 << (l1-1);      /* ..because the last digit of poly is always 1. */

  for ( i=1; i < np; i++, p1 += 2 ) {
    sum  = sum0;  /* initial value */
    mask = 2;
    for ( j=1; j < l1; j++, mask <<= 1 ) {
      if ( p1 & mask ) {
	sum.h ^= gf[j].h;
	sum.l ^= gf[j].l;
      }
    }
    if ( (sum.h == 0) && (sum.l == 0) )
      break;
  }
  free(gf);

  if ( i == np ) {
    fprintf( stderr, "find_longsubpoly(): no subpoly found!\n");    
    return 0L;
  }
  
  return p1;
}