RLMS Algorithm - Periodic Noise Cancellation

#define taps 8
#define Hz 0.763
#define mu 8e-3
#define Pz 0.817
#define Fz 0.137

void main()
{
     float c[taps], a[taps], b[taps], y[taps], xn[10], x1[taps], x[taps],y1[taps];
     float yn=0, rn, en, e_n, j, val;
     int reg[14] = { 0,1, 0,0,0,0, 0,1,0,0, 0,1,0,1 };
     int i, k, l=0, shift, value=0;

     clrscr();

     for( i=0; i<taps; i++ )
          y1[i] = x1[i] = x[i] = c[i] = y[i] = a[i] = b[i] = 0;

     for( i=0; i<200; i++ )
     {
          shift = (reg[0] ^ reg[2]) ^ (reg[9] ^ reg[11]);

          for( k=13; k>0; k-- )
               reg[k] = reg[k-1];
          reg[0] = shift;

          for( k=13; k>=0; k-- )
                value += reg[k] * pow( 2, k );
 
          value = value - 8192;
          yn = 5.0/8192 * value;

          value = 0;

          en = Hz * yn;

          for( k=taps-1; k>0; k-- )
               y[k] = y[k-1];

          y[0] = yn;

          rn = 0;
          for( k=0; k<taps; k++ )
               rn += c[k] * y[k];
 
          e_n = en - rn;

          for( k=0; k<taps; k++ )                     // NO FEEDBACK PATH
               c[k] = c[k] + ( mu * e_n * y[k] );    // MODELLING
     }

     for( i=0; i<taps; i++ )
          y[i] = 0;

     for( i=0,j=0; j<=5e-3; j+=5e-4,i++ )
     {
          val = 5 * sin( 2*3.142*200*j );
          xn[i] = val;
     }

     for( i=0; i<300; i++ )
     {
          shift = (reg[1] ^ reg[3]) ^ (reg[10] ^ reg[12]);

          for( k=13; k>0; k-- )
               reg[k] = reg[k-1];
          reg[0] = shift;

          for( k=13; k>=0; k-- )
                 value += reg[k] * pow( 2, k );
 
          value = value - 8192;

          for( k=taps-1; k>0; k-- )
               x[k] = x[k-1];
          x[0] = xn[l] + (yn * Fz);

          en = Hz*Pz*x[0];

          yn = 0;
          for( k=0; k<taps; k++ )
               yn += (a[k] * x[k]) + (b[k] * y[k]);
          rn = yn*Hz;

          for( k=taps-1; k>0; k-- )
               y[k] = y[k-1];
          y[0] = yn;

          for( k=taps-1; k>0; k-- )
               x1[k] = x1[k-1];
 
          x1[0] = 0;
          for( k=0; k<taps; k++ )
               x1[0] += c[k] * x[k];

          for( k=taps-1; k>0; k-- )
               y1[k] = y1[k-1];
 
          y1[0] = 0;
          for( k=0; k<taps; k++ )
               y1[0] += c[k] * y[k+1];

          e_n = en + rn;

          for( k=0; k<taps; k++ )
               a[k] = a[k] - ( mu * e_n * x1[k] );

          for( k=0; k<taps; k++ )
               b[k] = b[k] - ( mu * e_n * y1[k] );

          printf( "%.2f\n", e_n );

          l++;
          l = l>9 ? 0 : l;
     }
}