/* adp_sys_ident.c Adaptive system identification
                   using the LMS algorithm          */

# include <stdlib.h>
# include <stdio.h>
# define N0 48                  /* Adaptive filter order */
# define N1 48                  /* Unknown system order  */
# define Ns 128              /* Number of input signal*/

double w[N0],        /* Adaptive filter coefficients */
      d_sys[N0],    /* Adaptive filter signal buffer */
      d_fir[N1],    /* Unknown system signal buffer */
      in[Ns],       /* Input sample buffer */
      d[Ns],        /* Unknow system output buffer */
      LP_coef[N0]= {-0.00042,
-0.00110,
-0.00132, 
-0.00069,
0.00090,             
0.00285,
0.00371, 
0.00198,
0.00251, 
-0.00762,                            
-0.00946,
-0.00482,
0.00590,
0.01731,
0.02098,                          
0.01054,
-0.01282,
-0.03796,
-0.04723,
-0.02494,
0.03316, 
0.11491, 
0.19471, 
0.24395, 
0.24395,
0.19471, 
0.11491, 
0.03316, 
-0.02494, 
-0.04723,
-0.03796,
-0.01282,
0.01054,
0.02098, 
0.01731,             
0.00590,
-0.00482,
-0.00946,
-0.00762,
-0.00251, 
0.00198,
0.00371,
0.00285,
0.00090,
-0.00069,   
-0.00132,
-0.00110,
-0.00042};   /* Unknown system coefficients*/

void main()
{
   int i,j,k,ii;
   double y_u,y_a,e,uen,twomu;
   e=0.0;
   uen=0.0;
   twomu=0.00000000005;

   /* initialize coefficient & x[n],d[n] signal buffer*/
   for(i = 0; i <N0; i++)
   {
      d_sys[i] = 0.0;
   }
   for(i = 0; i <N1; i++)
      d_fir[i] = 0.0;
   for(i=0;i<N0;i++)
       {w[i]=0.0;}
   for(ii=0;ii<10000;ii++)          /* Generate sample and identify unknow system */
   {
	  for(i = 0; i < Ns; i++)
       {
       in[i] = rand()-0x4000;/* Zero-mean */
       }
	   
       
	  for(j = 0; j < Ns; j++)
       {
               d_fir[0] = in[j];
               y_u = 0.0;
	       for(i = 0; i < N0; i++)
		    {
                      if(j<N1 && j>=i)
	              y_u += LP_coef[i]*d_fir[j-i]; /* Convolution of d_fir[n] with h[n] */
                    
		      if(j>=N1)
			y_u += LP_coef[i]*d_fir[N1-i-1];
                     }
			  d[j] = y_u;/* Store y[n] to output buffer */
		   }  
		   for(i = N1-1; i>0 ; i--)
		   d_fir[i] = d_fir[i-1];   /* Shift old data d_fir[n-i] */
       }
	   for(j = 0; j < Ns; j++)
	   {
	      d_sys[0] = in[j];
	      y_a = 0.0;
	      for(i = 0; i < N0; i++)
		    {if(j<N0 && j>=i)
		     y_a += w[i]*d_sys[j-i];} /* Convolution of d_sys[n] with w[n] */
                     if(j<N0 && j>=i)
                     e = d[j] - y_a;
		     uen = twomu * e;
			 for(k=0;k<j;k++)
			 {
                          w[k]=w[k] + (uen * d_sys[j-k]);
                         } /* update coefficient */
		     {if(j>=N0)
		       y_a += w[i]*d_sys[N0-i-1]; /* Convolution of d_sys[n] with w[n] */
		     }
                      if(j>=N0) 
                      e = d[j] - y_a;
		      uen = twomu * e;
			  for(k=0;k<N0;k++)
			    { w[k]+=(uen * d_sys[N0-k-1]);
                            } /* update coefficient */

		  }

              for(i = N0-1; i > 0; i--)
	          d_sys[i] = d_sys[i-1];   /* Shift old data d_sys[n-i] */
	 }
    }
		  for (i =0;i<N0;i++)
		      printf("%20.17f\n",w[i]);
}