#include <stdio.h>
#include <math.h>
#include <fcntl.h>
#include <string.h>


int fd;
int steps;


#define MAXLENGTH 10*44100         /* 3 sec. samples enough ? */
#define MSEC(t) (int)(t*44.1)

short os[5][MAXLENGTH];

int consume_until(fd,str,len)
int fd, len;
char *str;
{
   int i;
   char b[4];
   read(fd,b,len);
/*   printf("[c,%s,%d]:%c%c%c%c,%d\n",str,len,b[0],b[1],b[2],b[3],
      strncmp(b,str,len)); */
   while (strncmp(b,str,len) != 0) {
      for (i=3; i>0; i--) b[i-1]=b[i];
      /* printf("[c ]:%d %d %d %d\n",(int) b[0],(int) b[1], (int) b[2], (int) b[3] ); */

      if (read(fd,b+3,1) != 1) {
         printf("Error: file format error\n");
         return(-1);
      }
   }
   return(0);                  /* match */
}

int readwav(s)                 /* Proto wav header, n=nnodes, filename */
char *s;
{
   char bb[4];
   short is, channels;
   int ii, length, samplerate;
   struct sample *p;
   fd = open(s,O_RDONLY|O_BINARY);
   if (fd < 0) return(-1);

   consume_until(fd,"RIFF",4);
   read(fd,&ii,4);
   consume_until(fd,"WAVE",4);

   consume_until(fd,"fmt ",4);
   read(fd,bb,4);   /* 0x10 */
   read(fd,bb,2);   /* 0x01 */
   read(fd,&channels,2);  /* (short) 1*/              /* Mono */
   if (channels != 1) {
      printf("Error: number of channels must be one.\n");
      return(-1);
   }
   read(fd,&samplerate,4);
   if (samplerate != 44100) {
      printf("Error: Sample rate should be 44100 Hz\n");
   }
   read(fd,&ii,4);
   read(fd,&is,2);
   read(fd,&is,2);

   consume_until(fd,"data",4);
   read(fd,&length,4);  /* 2* steps (length of data) */
   if (length/2 >= MAXLENGTH) {
      printf("Error: Memory limit %d samples exceeded\n",MAXLENGTH);
      exit(-1);
   }

   steps = length/2;
   if (read(fd,&(os[0][0]),length) != length) {
      printf("Error: not all data could be read of %s\n",s);
   }
   close(fd);
   return(0);
}



iwrite(fd,n,l)
int fd,l;
unsigned int n;
{
   write(fd,&n,l);
}

int initwav(s)                 /* Proto wav header, n=nnodes, filename */
char *s;
{
   fd = open(s,O_WRONLY|O_CREAT|O_TRUNC|O_BINARY);
   if (fd == NULL) return(-1);

   write(fd,"RIFF",4);
   iwrite(fd,(2*steps+36),4);
   write(fd,"WAVE",4);

   write(fd,"fmt ",4);
   iwrite(fd,(0x10),4);
   iwrite(fd,((short) 0x01),2);
   iwrite(fd,((short) 1),2);              /* Mono */
   iwrite(fd,(44100/1),4);                  /* Sample rate */
   iwrite(fd,(2*44100/1),4);
   iwrite(fd,((short) 2),2);
   iwrite(fd,((short) 16),2);

   write(fd,"data",4);
   iwrite(fd,(2*steps),4);
   return(0);
}


writewav(fd)
{
   int i;

   write(fd,&os[4][i],2*steps);
   close(fd);
}

clear_os()
{
   int i;
   for (i=0; i<steps; i++) {
       os[4][i] = os[3][i] = os[2][i] = os[1][i] = os[0][i] = 0;
   }
}

oldf()
{
   int i;
   int a,b,c,d, m, f,r, t, u;
   double tf, co1, co2, cu ;
      cu = ((float)(steps-i) * co1 + i*co2 ) /steps;
      t = (int) ((0x1000) * os[0][i]) - a ;
      /* t = 32000 - a; */
      tf = (142.48e-6 * (double) cu / 2.0);
      /* printf("t %d tf %f",t,tf); */
      tf *= (double) t;
      /* printf(" tf %f a %d",tf,(int) a); */
      a += (int) tf;
      /* os[4][i] += (short) ((float) t  * 142.48e-6 * cu / 2.0); */
      os[4][i] = (short) (a/0x1000);
      /* printf(" os %d a %d\n",(int) os[4][i],a); */
      /* a+=(f*(int)os[0][i]-(a)/0x10 s);
      os[4][i] = (short) (a/m); */
}

filter(co1, co2, re1, re2)
float co1, co2, re1, re2;
{
   int i;
   int a,b,c,d, m, f,r, t, u;
   double tf, cu, re;
   a=b=c=d=0;
   /* f = (int) (cu*((float) 0x1));
   r = (int) (re*((float) 0x100)); */
   m = 0x1000;
   /* co1 = 20000;
   co2 = 20; */
   for (i=0; i<steps; i++) {
      cu = ((float)(steps-i) * co1 + i*co2 ) /steps;
      re = ((float)(steps-i) * re1 + i*re2 ) /steps;
      t = (int) ((0x1000) * os[0][i]) - a ;
      tf = (142.48e-6 * (double) cu / 2.0);
      tf *= (double) t;
      a += (int) tf;
      os[4][i] = (short) (a/0x1000);
   }
    /*  os[i][4] = (os[i][0] + os[i+1][0] + os[i+2][0] + os[i+3][0])/4; */
}

filtern(co1, co2, re1, re2,n)
float co1, co2, re1, re2;
{
   int i, j;
   int b,c,d, m, f,r, t, u;
   int a[4];
   double tf, cu, re, mc;
   b=c=d=0;
   a[0] = a[1] = a[2] = a[3] = 0;
   /* f = (int) (cu*((float) 0x1));
   r = (int) (re*((float) 0x100)); */
   m = 0x1000;
   /* co1 = 20000;
   co2 = 20; */
   for (i=1; i<steps; i++) {
      cu = ((float)(steps-i) * co1 + i*co2 ) /steps;
      re = ((float)(steps-i) * re1 + i*re2 ) /steps;
      mc = ((4.0-fabs(re))/4.0);
      /* mc = 0.0; */
    for (j=0; j<n; j++) {
      t = (int) ((0x1000) * os[j][i]*mc) + re*((0x1000) * (os[(j-1+4)%4][i] + os[(j-1+4)%4][i-1])/2) - a[j] ;
      tf = (142.48e-6 * (double) cu / 2.0);
      tf *= (double) t;
      a[j] += (int) tf;
      os[4][i] = os[j+1][i] = (short) (a[j]/0x1000);
    }
   }
    /*  os[i][4] = (os[i][0] + os[i+1][0] + os[i+2][0] + os[i+3][0])/4; */
}


int main(argc,argv)
int argc;
char *argv[];
{
   float co1, co2, re1,re2;
   if (argc != 7) {
      printf("Usage: %s <infile.wav> <outfile.wav> <cutoff_1> <cutoff_2> <reso_1> <reso_2>\n",argv[0]);
      exit(-1);
   }
   clear_os();
   readwav(argv[1]);

   sscanf(argv[3],"%f",&co1);
   sscanf(argv[4],"%f",&co2);
   sscanf(argv[5],"%f",&re1);
   sscanf(argv[6],"%f",&re2);
   filtern(co1,co2,re1,re2,4);

   if (initwav(argv[2]) != 0) {
      printf("Cannot open output file %s.\n", argv[2]);
      exit(-3);
   }
   writewav(fd);
   return(0);
}