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

void main()
{
     float c[taps], d[taps], y[taps], xn[10], x1[taps], x[taps], w[taps];
     float x_n, un, fn, yn=0, rn, sn, 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          x1[i] = w[i] = x[i] = c[i] = y[i] = d[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;
          x_n = Fz * yn;

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

          y[0] = yn;

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

          for( k=0; k               c[k] = c[k] + ( mu * e_n * y[k] );

          sn=0;
          for( k=0; k               sn += d[k] * y[k];       // filter D(z)

          fn = x_n - sn;

          for( k=0; k               d[k] = d[k] + ( mu * fn * y[k] );
     }

     for( i=0; 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<30000; i++ )
     {
          for( k=taps-1; k>0; k-- )
               x[k] = x[k-1];
 
          un = xn[l] + (yn*Fz);          // u(n) = Actual i/p noise + Feedback

          sn=0;
          for( k=0; k               sn += d[k] * y[k];       // feedback estimate

          x[0] = un - sn;                // Actual noise recovered
          en = Hz*Pz*x[0];

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

          yn = 0;
          for( k=0; k               yn += w[k] * x[k];
          rn = yn*Hz;

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

          e_n = en + rn;

          for( k=0; k               w[k] = w[k] - ( mu * e_n * x1[k] );

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

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