/********************************************************************
*                                                                   *
* function    : LatticeWavefilters()                                *
* description : Lattice Wave Digital Filters                        *
*                                                                   *
*                                                                   *
*                                                                   *
* dependencies:                                                     *
* author      : Kenneth Blake -                                     *
*               I will maintain all rights applicable for this code,*
*               license and sofort according to swedish law         *
*               "upphovsratt" and what is called Copyright          *
*                                                                   *
*                                                                   *
*  input      :                                                     *
*                                                                   *
*                                                                   *
* history     :                                                     *
*                                                                   *
*               see                                                 *
*                   IEEE transactions on circuits and systems:      *
*                   vol CAS 32, no. 1, januari 1985                 *
*                   "Explicit formulas for lattice wave digital     *
*                    filters"                                       *
*                                                                   *
********************************************************************/

#include "math.h"


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

FILE *in1;        /* Pointer to file being used */
FILE *in2;        /* Pointer to file being used */
FILE *in3;        /* Pointer to file being used */
FILE *in4;        /* Pointer to file being used */
FILE *in5;        /* Pointer to file being used */
FILE *in6;        /* Pointer to file being used */
FILE *in7;        /* Pointer to file being used */

FILE *out1;        /* Pointer to file being used */

FILE *OUT1;        /* Pointer to file being used */
FILE *OUT2;        /* Pointer to file being used */
FILE *OUT3;        /* Pointer to file being used */
FILE *OUT4;        /* Pointer to file being used */

FILE *OUT5;        /* Pointer to file being used */
FILE *OUT6;        /* Pointer to file being used */
FILE *out8;        /* Pointer to file being used */



FILE *out   ;        /* Pointer to file being used */
char line[81];         /* Char array to hold lines read from file */



#define M_PI        3.14159265358979323846


#define N_MAX                21


#define KAPPA_NULL           0 
#define KAPPA_BIG_POSITIVE   1
#define KAPPA_SMALL_POSITIVE 2
#define KAPPA_BIG_NEGATIVE   3
#define KAPPA_SMALL_NEGATIVE 4

#define LW_LOWPASS     0
#define LW_HIGHPASS    1




typedef struct PORT 
{
  double alfa, kappa, delay;
  int property;

} port ;



typedef struct port2 
{
  double alfa, kappa, delay;
  int property;

} XFALL_latticewaveport ;






#define LW_LOWPASS     0
#define LW_HIGHPASS    1





double basic_first_degree_port(double in, port *p )
{
   double out, temp ;

   temp  =        -in            + p->delay ;
   out   =        p->kappa*temp  + p->delay ;
   p->delay =     out            - temp ;

   return out ;
}



double basic_second_degree_port(double in ,port *p1,port *p2)
{

   double out, temp ;

   temp  =        -p1->delay      + p2->delay ;
   out   =        p2->kappa*temp  + p2->delay ;
   p2->delay =     out            - temp ;

/* p1->delay = out ;                         */
   temp  =        -in             + out  ;
   out   =        p1->kappa*temp  + out  ;
   p1->delay =     out            - temp ;

   return out ;

}




void lattice_wave_5th_order(double in,  port pp[], double *lp,double *hp)
{

  double upper, lower ;

   upper = basic_first_degree_port( in,    &pp[0] ) ;
   upper = basic_second_degree_port(upper ,&pp[3],&pp[4]) ;

   lower = basic_second_degree_port(in ,&pp[1],&pp[2]) ;

   *lp = (lower+upper) * 0.5 ;
   *hp = (lower-upper) * 0.5 ;

}



void lattice_wave_5th_allpass(double in,  port pp[], double *up,double *low)
{

  double upper, lower ;

   upper = basic_first_degree_port( in,    &pp[0] ) ;
   upper = basic_second_degree_port(upper ,&pp[3],&pp[4]) ;

   lower = basic_second_degree_port(in ,&pp[1],&pp[2]) ;

   *up  = upper ;
   *low = lower ;

}


/*

This code you have to figure out yourself
it is from another implementation
but it shows the way of doing a highpass filter of order 
"order" = { 1,3,5,7,9,...  }


double filter(double in,int order)
{
  double upper, lower ;

  int i,j ;

  upper = first_degree_port(in, &pp[0]) ;
  
  lower = second_degree_port(in,    &pp[1],&pp[2]);

  for (i=2; i<=(order-1)/2;i++ )
  {
    j=2*i ;
    if (i%2==1)
     {
     lower = second_degree_port(lower, &pp[j-1],&pp[j]);
     }
    else
     upper = second_degree_port(upper, &pp[j-1],&pp[j]);
  }

 return ( (upper-lower) *0.5 ) ; 

}

*/

void clear_port(port p[], int order)
{
  short int i ;

  for (i=0 ; i< order  ; i++)
  {
     p[order].delay   =0.0 ;
  }
}



void  port_properties(port p[],int order)
{
    int i ;


     for (i=0 ; i<order ; i++)
     {
        if (p[i].kappa>0.0)
        {

           if (p[i].kappa <= 0.5)   p[i].property = KAPPA_SMALL_POSITIVE ;     /* 0.0  < kappa <= 0.5 */
           else                        p[i].property = KAPPA_BIG_POSITIVE   ;     /* 0.5  < kappa <  1.0 */
        }
        if(p[i].kappa<0.0)
        {
           if (p[i].kappa >= -0.5)  p[i].property = KAPPA_SMALL_NEGATIVE  ;    /* 0.0  < kappa <= -0.5 */
           else                        p[i].property = KAPPA_BIG_NEGATIVE    ;    /* -0.5 < kappa <  -1.0 */
        }

        if(p[i].kappa==0)
        {
                                       p[i].property = KAPPA_NULL     ;    /*  kappa = 0.0 */
        }

       switch (p[i].property)
       {
           case KAPPA_SMALL_POSITIVE : p[i].alfa =       p[i].kappa ; break ;  /* 0 < alfa<= 0.5 */
           case KAPPA_BIG_POSITIVE   : p[i].alfa = 1.0 - p[i].kappa ; break ;  /* 0 < alfa<= 0.5 */
           case KAPPA_SMALL_NEGATIVE : p[i].alfa =  fabs(p[i].kappa); break ;  /* 0 < alfa<= 0.5 */
           case KAPPA_BIG_NEGATIVE   : p[i].alfa = 1.0 + p[i].kappa ; break ;  /* 0 < alfa<= 0.5 */

       }

     }

}


void lw_butterworth_coefficients(port p[],double sampling_freq, double cutoff_freq, int order)
{
/*  alfa, property, kappa has the length order+1 */
     int i ;
     double fi_0, temp ;

     fi_0  = sin(M_PI*cutoff_freq/sampling_freq) / cos(M_PI*cutoff_freq/sampling_freq) ;

     p[0].kappa = (1.0 - fi_0 ) / (1.0 + fi_0 );

     for (i = 1; i <= (order-1)/2; i++) {

      temp = sin((2.0 *M_PI* cutoff_freq) / sampling_freq) * cos( (M_PI*i)/order );

          p[(2*i)-1].kappa = (temp - 1.0) / (temp + 1.0);

          p[(2*i)].kappa   = cos( (2.0 * M_PI* cutoff_freq)/sampling_freq );
     }

    port_properties(p,order);

}


void lw_chebyshev_coefficients(port p[],double F, double fp, int order)
{
     double Es, Ep ,fs, fi_p, fi_s,w,r,N, Ai,Bi;
     short int i ;
     Es = 0.01 ;
     Ep = 0.01 ;

    N = order ;
    fi_p = /* transformed frequency, passband */  tan(M_PI*fp/F) ;
    fi_s = /* transformed frequency, passband */  tan(M_PI*fs/F) ;



     /* help variables */

    w  = pow(     sqrt( sqrt(1/(Ep*Ep) + 1) + 1/Ep ), 1/N );
    r  = ( w-1/w) * fi_p;

    order=N ;
    p[0].kappa = (2-r) / (2+r) ;

     for (i = 1; i <= (order-1)/2; i++)
     {
          Ai = r*cos( M_PI*i/N) ;
          Bi = ( w*w + 1/(w*w) - 2.0*cos(2.0*M_PI*i/N) ) * fi_p*fi_p/4;
          p[(2*i)-1].kappa = ( Ai -Bi -1) /(Ai+Bi+1) ;
          p[(2*i)].kappa   = ( 1-Bi) / (1+Bi) ;
     }
    port_properties(p,order);

}


void lw_cauer_coefficients(port p[],double F, double fp, int order)


{
     double Es, Ep ,fs, fi_p, fi_s,w,r,N, Ai,Bi,Yi,fsmin,as;
     double r0,r1,r2;
     double x0,x1,x2,x3,x4;
     double q0,q1,q2,q3,q4;
     double m0,m1,m2,m3;
     double g1,g2,g3;
     double w0,w1,w2,w3,w4,w5;
     double c0,c1,c2,c3,c4;

     short int i ;
     Ep = 0.01 ;


     as = 6.0*order;

     Es = pow(  pow(10.0,(as/10.0))-1, 0.5 ) ;



    N = order ;
    fi_p = /* transformed frequency, passband */  tan(M_PI*fp/F) ;


    N = order ;
     /* help variables */
    r0 = pow( (Es/Ep), 0.5 )           ;
    r1 = r0 + pow ( pow(r0,4.0) -1, 0.5 );
    r2 = r1 + pow ( pow(r1,4.0) -1, 0.5 );

    x4 = 0.5* pow(  pow(2.0*r2,1.0/N),4.0)  ;
    x3 =    pow((x4+1/x4)*0.5,0.5) ;
    x2 =    pow((x3+1/x3)*0.5,0.5) ;
    x1 =    pow((x2+1/x2)*0.5,0.5) ;
    x0 =    pow((x1+1/x1)*0.5,0.5) ;

    fsmin= (F/M_PI) * atan( fi_p *x0*x0) ;


    fi_s = /* transformed frequency, passband */ tan(M_PI*fsmin/F) ;

    q0  = pow( ( fi_s/fi_p), 0.5 ) ;         /* q0 */
    q1  = q0*q0 + pow(  (q0*q0*q0*q0-1), 0.5 ) ;   /* q1 */
    q2  = q1*q1 + pow(  (q1*q1*q1*q1-1), 0.5 ) ;   /* q2 */
    q3  = q2*q2 + pow(  (q2*q2*q2*q2-1), 0.5 ) ;   /* q3 */
    q4  = q3*q3 + pow(  (q3*q3*q3*q3-1), 0.5 ) ;   /* q4 */

    m3  = 0.5 * ( pow(  pow(2.0*q4,0.5), N ) ) ;
    m2  = pow((0.5*(m3 +1/m3)),0.5);
    m1  = pow((0.5*(m2 +1/m2)),0.5);
    m0  = pow((0.5*(m1 +1/m1)),0.5);


    g1 = 1/Ep + pow( (1+1/(Ep*Ep)), 0.5);
    g2 = m1*g1 + pow( m1*g1*m1*g1+1,0.5 ) ;
    g3 = m2*g2 + pow( m2*g2*m2*g2+1,0.5 ) ;


    w5 = m3*m3/(g3*g3) + 1 ;
    w5 = pow( w5, 0.5 ) ;
    w5 = w5 + m3/g3;
    w5 = pow( w5, 1.0/N);

    w4 = ( w5 - 1/w5) /(2.0*q4) ;
    w3 = ( w4 - 1/w4) /(2.0*q3) ;
    w2 = ( w3 - 1/w3) /(2.0*q2) ;
    w1 = ( w2 - 1/w2) /(2.0*q1) ;
    w0 = ( w1 - 1/w1) /(2.0*q0) ;


    p[0].kappa = ( 1 + w0*q0*fi_p) / ( 1 - w0*q0*fi_p) ;





    order=N ;

     for (i = 1; i <= (order-1)/2; i++)
     {

         c4 = q4/(sin( i*M_PI/N) ) ;
         c3 = ( c4 +1/c4) / ( 2.0*q3) ;
         c2 = ( c3 +1/c3) / ( 2.0*q2) ;
         c1 = ( c2 +1/c2) / ( 2.0*q1) ;
         c0 = ( c1 +1/c1) / ( 2.0*q0) ;

         Yi =1/c0;

         Bi = ( w0*w0 + Yi*Yi)*( q0*q0*fi_p*fi_p)/( 1+ w0*w0*Yi*Yi) ;

         Ai = 1 - ( q0*q0+ 1/(q0*q0) -Yi*Yi)* Yi*Yi;
         Ai = pow(Ai,0.5);

         Ai = -Ai*2.0*w0*q0*fi_p /( 1+ w0*w0*Yi*Yi) ;

          p[(2*i)-1].kappa = ( Ai -Bi -1) /(Ai+Bi+1) ;
          p[(2*i)].kappa   = ( 1-Bi) / (1+Bi) ;
     }
    port_properties(p,order);

}






#define FFT_LEN       512 /* OBSERVE the method of using the Brigham        */
                                  /*         unpacking means that a fft length of   */
                                  /*         512 actually is a length of 1024       */
                                  /* Brigham: The fast fourier transform and its    */
                                  /* applications ISBN 0 13 307547-8 Prentice  Hall */
                          /* Page 188 onwards sect 9.3 fft alg for real data*/


port bu[7],ch[7], ca[7] ,high[7], allpass[7];

#define LENGTH 1030 
double temp[2000], t ,

Low_Pass_butterworth[LENGTH],
High_Pass_butterworth[LENGTH],
Low_Pass_chebyshev[LENGTH],
High_Pass_chebyshev[LENGTH],
Low_Pass_cauer[LENGTH],
High_Pass_cauer[LENGTH];


double temp[2000], t, 
       LP_magnitude[1024] ,  
       HP_magnitude[1024] ,  
       sinustable[FFT_LEN]; /* address of real value array*/
int   bitreverse[FFT_LEN]; /* address of real value array*/


/* ----------------------------------------------------------------- */
/* Calculate the table consisting of sinus values                    */
/* A table is always faster than a calculation                       */
/* ----------------------------------------------------------------- */
static 
void 
calc_sinus_table(int n)
{
  int i ;
   for (i=0; i<n ; i++) {
      sinustable[i] = sin(2.0*3.141592654 * i /n) ;
   }
}

/* ----------------------------------------------------------------- */
/* Calculate the bitreversed number                                  */
/* ----------------------------------------------------------------- */
static 
int 
bitrev(int rev, int fft_len)
{
 int upwards,downwards;

 int sum;
 sum=0;

      sum=0 ;

      upwards = 1 ;
      downwards = fft_len>>1 ;
  
      while( downwards > 0 )
      {
         if (upwards & rev) sum = sum + downwards ;
         upwards   = upwards   << 1;
            downwards = downwards >> 1;
        }
 return sum;
}

/* ----------------------------------------------------------------- */
/* Calculate the table consisting of bitreverse and zero values      */
/* The zero elements does not need swapping                          */
/* They have already been swapped                                    */
/* A table is always faster than a calculation                       */
/* ----------------------------------------------------------------- */
static 
void 
create_bitreverse_table(int len)
{
   int i,z ;

   for (i=0;i<len;i++) {
      z = bitrev(i,len);
      if (i<z) bitreverse[i] =z;
      else bitreverse[i] =0;
   }
}



/* ----------------------------------------------------------------- */
/* Cooley Tukey algorithm with output bitreversal                    */
/* The zero elements does not need swapping                          */
/* They have already been swapped                                    */
/* A table is always faster than a calculation                       */
/* ----------------------------------------------------------------- */
static 
void 
fft2DIF(double *re, double *im, int n)
{
    int p, q, n1, n2 , j ,w0,i,z,t,w;
    double wr, wi ;
    double tr,ti ;

    w0=1;
    n2 = n ;
    t=n/4;

    while (n2>0)
    {
            n1 = n2 ;
            n2 = n2 >> 1 ;
            j = 0 ;
            while (j != n2)
            {
                    w  = j*w0 ;
                    wr =  sinustable[t + w] ;
                    wi = -sinustable[w];
                    p=j;
                    while ( p<n )
                    {
                        q= p + n2 ;
                        tr       = re[p] - re[q] ;
                        ti       = im[p] - im[q] ;
                        re[p]   = re[p] + re[q] ;
                        im[p]   = im[p] + im[q] ;
                        re[q]   = tr*wr  - ti*wi  ;               /* x(p) =   x(p) + x(q)          */
                        im[q]   = tr*wi  + ti*wr  ;               /* x(q) = ( x(p) - x(q) )*WrN    */
                        p = p + n1 ;
                    }
                    j=j+1;
            }
          w0 = w0*2 ;
    }

/* ----------------------------------------------------------------- */
/* Swapp the elements according to bitreverse calculations           */
/* The zero elements does not need swapping                          */
/* They have already been swapped                                    */
/* A table is always faster than a calculation                       */
/* ----------------------------------------------------------------- */
    for(i=0 ; i < (n) ; i++ )
    {
        z  = bitreverse[i] ;
        if(z>0) {
            tr = re[i] ;
            ti = im[i] ;
            re[i] = re[z] ;
            im[i] = im[z] ;
            re[z] = tr ;
            im[z] = ti ;
      }
   }
}

void Brigham_split(double r[], double i[], short int N, double power[])
{
   double re,im, w ,scale ;
   int n;
/* ----------------------------------------------------------------- */
/* As the input is real, half the energy of the FFT is in the points */
/* 512 to 1023. We have thus to multiply the bin inbetween 0 to 511  */
/* with two to get the correct TOTAL POWER value                     */
/*                                                                   */
/* ----------------------------------------------------------------- */

   scale= 1.0/(2.0*N);
/*  scale= scale*scale*2.0 ; /* one side 0=> *2.0  */


   r[N] = r[0] ;
   i[N] = i[0] ;

   for ( n=0 ; n < N ; n++ ) {
       w  = M_PI*n/N ;
       re = (r[n]  + r[N-n] + cos(w) * ( i[n]  + i[N-n]) - sin(w) * (r[n] - r[N-n]) ) * 0.5   ;
       im = (i[n]  - i[N-n] - sin(w) * ( i[n]  + i[N-n]) - cos(w) * (r[n] - r[N-n]) ) * 0.5   ;
       power[n] = sqrt( (re*re + im*im) *scale );
   }
}



void  PSD(double *input, double *output)
{
         
  double   *profile_array ;
  double   rez[FFT_LEN+1] ;
  double   imz[FFT_LEN+1] ;
  int i;

/* ----------------------------------------------------------------- */
/* The input length thus contain               */
/* 1024 elements                                                     */
/*                                                                   */
/* The FFT is over 1024 points (1024 meter)                         */
/*                                                                   */
/* Copy the profile into the FFT array                               */
/*                                                                   */
/* ----------------------------------------------------------------- */


   profile_array = input;
   for (i=0 ; i< 512; i++)
   {
         rez[i]   = *profile_array++; 
         imz[i]   = *profile_array++;
   }

  fft2DIF( rez,imz, FFT_LEN) ;
  Brigham_split(rez,imz, FFT_LEN, output) ;
}



void ChangeDot_to_comma_and_delete_temp_txt_file(char *filename)
{

  FILE *in1;        /* Pointer to file being used */
  FILE *OUT1;        /* Pointer to file being used */

   if ((OUT1= fopen( filename, "w")) == NULL) {
      printf("Error opening text file for writing\n");
      exit(0);
   }

   if ((in1 = fopen("temp.txt", "r")) == NULL) {
      printf("Error opening text file for writing\n");
      exit(0);
   }

/* change '.' to ',' in text file so that excel can use it*/
   {
   int byte ;
     do 
     {
        byte = fgetc(in1) ;
        if ( 0 != byte ) 
        {
          if ( '.' == byte )
          byte = ',' ;
         }
        if (EOF != byte) fputc( byte, OUT1); 
          
     }
     while  ( EOF != byte) ;  
   }  
   
   
  fclose(OUT1);
  fclose(in1);
  remove("temp.txt")   ;

}

port butterworth[7],chebyshev[7], cauer[7];




void main()
{
  int i,even,odd,j , order, t;

  double in, lp,hp ,upper, lower;

calc_sinus_table(FFT_LEN)        ;
create_bitreverse_table(FFT_LEN) ;

 order = 5 ;
 lw_butterworth_coefficients(&butterworth[0],    20000.0, 5000.0,  order) ;
 lw_chebyshev_coefficients(  &chebyshev[0] ,     20000.0, 5000.0, order) ;
 lw_cauer_coefficients    (  &cauer[0],          20000.0, 5000.0, order) ;


  clear_port(&butterworth[0],    order);
  clear_port(&chebyshev[0], order);
  clear_port(&cauer[0],     order);



/*   forward filtering */

     for (i=0; i<500;i++)
     {
       lattice_wave_5th_order(0.0, &butterworth[0], &Low_Pass_butterworth[i], &High_Pass_butterworth[i]) ;
     }

/* Dirac pulse */
     lattice_wave_5th_order(32.0, &butterworth[0], &Low_Pass_butterworth[i], &High_Pass_butterworth[i]) ;

     for (i=501; i<1024;i++)
     {
       lattice_wave_5th_order(0.0,&butterworth[0], &Low_Pass_butterworth[i], &High_Pass_butterworth[i]) ;
     }

   if ((OUT1= fopen("temp.txt", "w")) == NULL) {
      printf("Error opening text file for writing\n");
      exit(0);
   }

     fprintf(OUT1,"Butterworth_NON_LINEAR_Lowpass NON_LINEAR_highpass\n");
     for (i=0; i<1024; i++)
     {
       fprintf(OUT1,"%f %f\n", Low_Pass_butterworth[i], High_Pass_butterworth[i]);
     }

     fclose(OUT1);
   
    ChangeDot_to_comma_and_delete_temp_txt_file("Butterworth Impulse response NONLinear_phase.txt") ;

/* ************************************************ */

    PSD(&Low_Pass_butterworth[0], &LP_magnitude[0]) ;
    PSD(&High_Pass_butterworth[0], &HP_magnitude[0]) ;

    if ((OUT1= fopen("temp.txt", "w")) == NULL) {
      printf("Error opening text file for writing\n");
      exit(0);
   }

     fprintf(OUT1,"Butterworth_NON_LINEAR_Lowpass Butterworth_NON_LINEAR_highpass\n");
     for (i=0; i<FFT_LEN; i++)
     {
       fprintf(OUT1,"%f %f\n", LP_magnitude[i],HP_magnitude[i]);
     }

     fclose(OUT1);
   
    ChangeDot_to_comma_and_delete_temp_txt_file("Butterworth NONlinear magnitude.txt") ;
/* ************************************************ */





/* rewerse filtering Highpass*/
     clear_port(&butterworth[0], 5);
     for (i=1024; i>0; i--)
     {
       lattice_wave_5th_order(High_Pass_butterworth[i],&butterworth[0] , &in, &High_Pass_butterworth[i]) ;
     }


/* rewerse filtering low pass*/
     clear_port(&butterworth[0], 5);

     for (i=1024; i>0; i--)
     {
       lattice_wave_5th_order(Low_Pass_butterworth[i],&butterworth[0] , &Low_Pass_butterworth[i], &in) ;
     }


   if ((OUT1= fopen("temp.txt", "w")) == NULL) {
      printf("Error opening text file for writing\n");
      exit(0);
   }

     fprintf(OUT1,"Butterworth_LINEAR_Lowpass Butterworth_LINEAR_highpass\n");
     for (i=0; i<1024; i++)
     {
       fprintf(OUT1,"%f %f\n", Low_Pass_butterworth[i], High_Pass_butterworth[i]);
     }

     fclose(OUT1);
   
    ChangeDot_to_comma_and_delete_temp_txt_file("Butterworth Impulse response Linear_phase.txt") ;


/* ************************************************ */

    PSD(&Low_Pass_butterworth[0], &LP_magnitude[0]) ;
    PSD(&High_Pass_butterworth[0], &HP_magnitude[0]) ;

    if ((OUT1= fopen("temp.txt", "w")) == NULL) {
      printf("Error opening text file for writing\n");
      exit(0);
   }

     fprintf(OUT1,"Butterworth_LINEAR_Lowpass Butterworth_LINEAR_highpass\n");
     for (i=0; i<FFT_LEN; i++)
     {
       fprintf(OUT1,"%f %f\n", LP_magnitude[i],HP_magnitude[i]);
     }

     fclose(OUT1);
   
    ChangeDot_to_comma_and_delete_temp_txt_file("Butterworth linear magnitude.txt") ;
/* ************************************************ */



/*





  CHEBYSHEV






*/



/*   forward filtering */

     for (i=0; i<500;i++)
     {
       lattice_wave_5th_order(0.0, &chebyshev[0], &Low_Pass_chebyshev[i], &High_Pass_chebyshev[i]) ;
     }

/* Dirac pulse */
     lattice_wave_5th_order(32.0, &chebyshev[0], &Low_Pass_chebyshev[i], &High_Pass_chebyshev[i]) ;

     for (i=501; i<1024;i++)
     {
       lattice_wave_5th_order(0.0,&chebyshev[0], &Low_Pass_chebyshev[i], &High_Pass_chebyshev[i]) ;
     }

   if ((OUT1= fopen("temp.txt", "w")) == NULL) {
      printf("Error opening text file for writing\n");
      exit(0);
   }

     fprintf(OUT1,"Chebyshev_NON_LINEAR_Lowpass Chebyshev_NON_LINEAR_highpass\n");
     for (i=0; i<1024; i++)
     {
       fprintf(OUT1,"%f %f\n", Low_Pass_chebyshev[i], High_Pass_chebyshev[i]);
     }

     fclose(OUT1);
   
    ChangeDot_to_comma_and_delete_temp_txt_file("chebyshev Impulse response NONLinear_phase.txt") ;

/* ************************************************ */

    PSD(&Low_Pass_chebyshev[0], &LP_magnitude[0]) ;
    PSD(&High_Pass_chebyshev[0], &HP_magnitude[0]) ;

    if ((OUT1= fopen("temp.txt", "w")) == NULL) {
      printf("Error opening text file for writing\n");
      exit(0);
   }

     fprintf(OUT1,"Chebyshev_NON_LINEAR_Lowpass Chebyshev_NON_LINEAR_highpass\n");
     for (i=0; i<FFT_LEN; i++)
     {
       fprintf(OUT1,"%f %f\n", LP_magnitude[i],HP_magnitude[i]);
     }

     fclose(OUT1);
   
    ChangeDot_to_comma_and_delete_temp_txt_file("chebyshev NONlinear magnitude.txt") ;
/* ************************************************ */





/* rewerse filtering high pass*/
     clear_port(&chebyshev[0], 5);
     for (i=1024; i>0; i--)
     {
       lattice_wave_5th_order(High_Pass_chebyshev[i],&chebyshev[0] , &in, &High_Pass_chebyshev[i]) ;
     }


/* rewerse filtering low pass*/
     clear_port(&chebyshev[0], 5);

     for (i=1024; i>0; i--)
     {
       lattice_wave_5th_order(Low_Pass_chebyshev[i],&chebyshev[0] , &Low_Pass_chebyshev[i], &in) ;
     }


   if ((OUT1= fopen("temp.txt", "w")) == NULL) {
      printf("Error opening text file for writing\n");
      exit(0);
   }

     fprintf(OUT1,"Chebyshev_LINEAR_Lowpass Chebyshev_LINEAR_highpass\n");
     for (i=0; i<1024; i++)
     {
       fprintf(OUT1,"%f %f\n", Low_Pass_chebyshev[i], High_Pass_chebyshev[i]);
     }

     fclose(OUT1);
   
    ChangeDot_to_comma_and_delete_temp_txt_file("chebyshev Impulse response Linear_phase.txt") ;


/* ************************************************ */

    PSD(&Low_Pass_chebyshev[0], &LP_magnitude[0]) ;
    PSD(&High_Pass_chebyshev[0], &HP_magnitude[0]) ;

    if ((OUT1= fopen("temp.txt", "w")) == NULL) {
      printf("Error opening text file for writing\n");
      exit(0);
   }

     fprintf(OUT1,"Chebyshev_LINEAR_Lowpass Chebyshev_LINEAR_highpass\n");
     for (i=0; i<FFT_LEN; i++)
     {
       fprintf(OUT1,"%f %f\n", LP_magnitude[i],HP_magnitude[i]);
     }

     fclose(OUT1);
   
    ChangeDot_to_comma_and_delete_temp_txt_file("chebyshev linear magnitude.txt") ;
/* ************************************************ */



/*


  CAUER


*/




/*   forward filtering */

     for (i=0; i<500;i++)
     {
       lattice_wave_5th_order(0.0, &cauer[0], &Low_Pass_cauer[i], &High_Pass_cauer[i]) ;
     }

/* Dirac pulse */
     lattice_wave_5th_order(32.0, &cauer[0], &Low_Pass_cauer[i], &High_Pass_cauer[i]) ;

     for (i=501; i<1024;i++)
     {
       lattice_wave_5th_order(0.0,&cauer[0], &Low_Pass_cauer[i], &High_Pass_cauer[i]) ;
     }

   if ((OUT1= fopen("temp.txt", "w")) == NULL) {
      printf("Error opening text file for writing\n");
      exit(0);
   }

     fprintf(OUT1,"CAUER_NON_LINEAR_Lowpass CAUER_NON_LINEAR_highpass\n");
     for (i=0; i<1024; i++)
     {
       fprintf(OUT1,"%f %f\n", Low_Pass_cauer[i], High_Pass_cauer[i]);
     }

     fclose(OUT1);
   
    ChangeDot_to_comma_and_delete_temp_txt_file("cauer Impulse response NONLinear_phase.txt") ;

/* ************************************************ */

    PSD(&Low_Pass_cauer[0], &LP_magnitude[0]) ;
    PSD(&High_Pass_cauer[0], &HP_magnitude[0]) ;

    if ((OUT1= fopen("temp.txt", "w")) == NULL) {
      printf("Error opening text file for writing\n");
      exit(0);
   }

     fprintf(OUT1,"CAUER_NON_LINEAR_Lowpass CAUER_NON_LINEAR_highpass\n");
     for (i=0; i<FFT_LEN; i++)
     {
       fprintf(OUT1,"%f %f\n", LP_magnitude[i],HP_magnitude[i]);
     }

     fclose(OUT1);
   
    ChangeDot_to_comma_and_delete_temp_txt_file("cauer NONlinear magnitude.txt") ;
/* ************************************************ */





/* rewerse filtering high pass*/
     clear_port(&cauer[0], 5);
     for (i=1024; i>0; i--)
     {
       lattice_wave_5th_order(High_Pass_cauer[i],&cauer[0] , &in, &High_Pass_cauer[i]) ;
     }


/* rewerse filtering low pass*/
     clear_port(&cauer[0], 5);

     for (i=1024; i>0; i--)
     {
       lattice_wave_5th_order(Low_Pass_cauer[i],&cauer[0] , &Low_Pass_cauer[i], &in) ;
     }


   if ((OUT1= fopen("temp.txt", "w")) == NULL) {
      printf("Error opening text file for writing\n");
      exit(0);
   }

     fprintf(OUT1,"CAUER_LINEAR_Lowpass CAUER_LINEAR_highpass\n");
     for (i=0; i<1024; i++)
     {
       fprintf(OUT1,"%f %f\n", Low_Pass_cauer[i], High_Pass_cauer[i]);
     }

     fclose(OUT1);
   
    ChangeDot_to_comma_and_delete_temp_txt_file("cauer Impulse response Linear_phase.txt") ;


/* ************************************************ */

    PSD(&Low_Pass_cauer[0], &LP_magnitude[0]) ;
    PSD(&High_Pass_cauer[0], &HP_magnitude[0]) ;

    if ((OUT1= fopen("temp.txt", "w")) == NULL) {
      printf("Error opening text file for writing\n");
      exit(0);
   }

     fprintf(OUT1,"CAUER_LINEAR_Lowpass CAUER_LINEAR_highpass\n");
     for (i=0; i<FFT_LEN; i++)
     {
       fprintf(OUT1,"%f %f\n", LP_magnitude[i],HP_magnitude[i]);
     }

     fclose(OUT1);
   
    ChangeDot_to_comma_and_delete_temp_txt_file("cauer linear magnitude.txt") ;
/* ************************************************ */






















 exit(0) ;




   
}