rbds.c

error correction code

// ------------------------------------------------------------------------
// An example of error trapping decoding of a cyclic code.
//
// This is a shortened cyclic (26,16) code capable of correcting 
// any random burst of up to 5 errors. This code was invented by 
// Professor Tadao Kasami in 1963:
// T. Kasami, ``Optimum Shortened Cyclic Codes for Burst-Error Correction,''
// IEEE Trans. Info. Theory, vol. IT-9, no. 4, pp. 105-109, July 1963.
// (See Table I, b=5, first column, generator = 2671 in octal, n = 27.)
//
// This code is used in the U.S. Radio Broadcast Data System (RBDS) standard,
// National Radio Systems Committee, April 9, 1998
//
// The main idea behind this program is to emulate the operation of the
// encoder and decoder of a binary cyclic codes, using bitwise shifts
// and xor for modulo g(x) operation
//
// Compiling (in Linux/UNIX) with
//                 gcc -DPRINT_GEN -o RBDS.exe RBDS.c -lm
// gives the generator polynomial matrix of the code, as columns 1 and 2 of
// the output, the received word syndrome in the third column and the decoded
// message in the fourth column, as shown below:
//
//                       M    S        Sr        Mhat
//                     8000  77 --- > 3d2  --- > 8000
//                     4000 2e7 --- > 3d2  --- > 4000
//                     2000 3af --- > 3d2  --- > 2000
//                     1000 30b --- > 3d2  --- > 1000
//                      800 359 --- > 3d2  --- >  800
//                      400 370 --- > 3d2  --- >  400
//                      200 1b8 --- > 3d2  --- >  200
//                      100  dc --- > 3d2  --- >  100
//                       80  6e --- > 3d2  --- >   80
//                       40  37 --- > 3d2  --- >   40
//                       20 2c7 --- > 3d2  --- >   20
//                       10 3bf --- > 3d2  --- >   10
//                        8 303 --- > 3d2  --- >    8
//                        4 35d --- > 3d2  --- >    4
//                        2 372 --- > 3d2  --- >    2
//                        1 1b9 --- > 3d2  --- >    1
// ------------------------------------------------------------------------
// 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 <stdio.h>
#include <math.h>
#include <float.h>
#include <limits.h>

#define MAXRAND LONG_MAX    // for random number generation
#define G     0x5b9         // x^{10} + x^8 + x^7 + x^5 + x^4 + x^3 + 1
#define P     0x31b         // x^9 + x^8 + x^4 + x^3 + x + 1
#define K21   32767         // 2^{k-1} - 1
#define K2    32768         // 2^{k-1} 
#define FLAG   1024         // 2^{n-k}
#define NK2     512         // 2^{n-k-1}
#define NK    0x7ff         // n-k+1 ones
#define TRAP  0x01f         // trapping pattern = 0000011111
#define RED2      9         // n-k-2

main()
{
int M;                      // message
int St;                     // syndrome of transmitted message
int S;                      // syndrome register
int Mr;                     // received message
int Mhat;                   // estimated message
int Sr;                     // receiver syndrome register
int bit;                    // information bit
int i, j, k, ell;           // auxiliary variables
int error;

long seed;

  printf("\nSimulation of encoding and decoding for burst error correction with a\n");
  printf("shortened (26,16) cyclic code with generator polynomial (hexa) G = %x\n\n", 
         G);
  time(&seed);
  srandom(seed);

#ifdef PRINT_GEN
  printf("  M    S        Sr        Mhat\n");
#endif

  for (M=32768; M>=1; M=M/2) {

  // ---------------------------------------------------------------
  //                          ENCODER
  // Shift register implementation: Message is divided by G. In the
  // reminder, stored in the syndrome register S, are the redundant
  // bits.
  // ---------------------------------------------------------------

  S = 0;
  j = K2;
  for (i=15; i>=0; i--) {
    bit = ( (j & M) >> i ) & 0x01;
    j = j>>1;
    if (bit)
      S = ((S<<1)^FLAG) & NK;
    else
      S = (S<<1) & NK;
    if (S & FLAG)
      S = (S^G) & NK;
    }

#ifdef PRINT_GEN
  // Write the generator matrix of the code:
  printf("%4x %3x", M, S);
#endif

  // ---------------------------------------------------------------
  //                          CHANNEL
  // ---------------------------------------------------------------

  Mr = M ^ 0x001f; // A burst of errors of length 5
#ifndef PRINT_GEN
  printf("Transmitted = %4x (S=%3x)\n", M, S);
  printf("Received    = %4x \n", Mr);
#endif

  // ---------------------------------------------------------------
  //                          DECODER
  // ---------------------------------------------------------------
  // STAGE 1: Feed the information bits
  // Multiply by P and divide by G: Same as encoder, except for multiply

  Sr = 0;
  j = K2;
  for (i=15; i>=0; i--) {
    bit = ( (j & Mr) >> i ) & 0x01;
    j = j>>1;
    if (bit)
      Sr = ((Sr<<1)^P) & NK;
    else
      Sr = (Sr<<1) & NK;
    if (Sr & FLAG)
      Sr = (Sr^G) & NK;
    }

  // STAGE 2: Feed the redundant bits into the receiver shift register
  j = NK2;
  for (i=9; i>=0; i--) {
    bit = ( (j & S) >> i ) & 0x01;
    j = j>>1;
    if (bit)
      Sr = ((Sr<<1)^P) & NK;
    else
      Sr = (Sr<<1) & NK;
    if (Sr & FLAG)
      Sr = (Sr^G) & NK;
    }

#ifdef PRINT_GEN
   printf(" --- > %3x ", Sr);
#endif

  // STAGE 3: Syndrome checking and error correction of information
  j = K2;
  Mhat = 0;
  for (i=15; i>=0; i--) {
    bit = ( (j & Mr) >> i ) & 0x01;
    // If errors are trapped correct them
    if (!(Sr & TRAP)) {
      error = ((Sr & NK2) >> RED2) & 0x01;
      bit ^= error;
      Sr = (Sr<<1) & NK;
      }
    else {
      Sr = (Sr<<1) & NK;
      if (Sr & FLAG)
        Sr = (Sr^G) & NK;
      }
    Mhat += (bit*j);
    j = j>>1;
    }

#ifdef PRINT_GEN
  printf(" --- > %4x\n", Mhat);
#endif

#ifndef PRINT_GEN
  printf("Corrected   = %4x\n", Mhat);
#endif

  }  // end for M

}