/* Multichannel FIR filter*/

#include "stdio.h"
#include "math.h"
main () {
  int data[1025][30], dout[1025][30], bg[1025];
  float b[15], b1[15], xm[32], initial, sum, pn, ps, snr;
  int i, j, k;
  FILE *fp1, *fp2, *fp3; 
  static char rec_file[]="record01.dat";


  fp1=fopen ("9294all.asc","r");
  for (i=0; i<1025; i++) {
     for (j=0; j<29; j++) {
	fscanf (fp1, "%d ", &data[i][j]);  
     }
     fscanf (fp1, "%d", &data[i][29]);
  }
  fclose (fp1);

  fp2=fopen ("9294all.check", "w");
  for (i=0; i<1025; i++) {
     for (j=0; j<29; j++) {
	fprintf (fp2, "%d ", data[i][j]);
     }
     fprintf (fp2, "%d\n", data[i][29]);
  }
  fclose (fp2);

  fp3=fopen ("b", "r");
  for (i=0; i<15; i++) {
     fscanf (fp3, "%f ", &b[i]);
  }
  fclose (fp3);

  for (i=0; i<1025; i++) {
     initial+=data[i][0];
  }
  initial=initial/1025;

  for (i=0; i<1025; i++) {
     for (k=0; k<15; k++) {
	xm[k]=initial/3;
     }
     for (j=0; j<30; j++) {
	xm[0]=data[i][j];
	sum=0;
	for (k=14; k>0; k--) {
	  sum=sum+b[k]*xm[k];
	  xm[k]=xm[k-1];
	}
	sum=sum+b[0]*xm[0]; 
	dout[i][j]=sum+0.5;
     }
  }
     
  for (i=0; i<1025; i++) {
     initial+=dout[i][0];
  }
  initial=initial/1025;

  for (j=0; j<30; j++) {
     for (k=0; k<32; k++) {
	xm[k]=initial/2;
     }
     for (i=0; i<1025; i++) {
        if (dout[i][j]<8*initial) {
	   xm[0]=dout[i][j];
	}
	else {
	   xm[0]=initial;
	}
	sum=0;
	for (k=59; k>0; k--) {
	  sum=sum+xm[k];
	  xm[k]=xm[k-1];
	}
	sum=sum+xm[0]; 
	sum=sum/60;
	bg[i]=sum+0.5;
     }

     sum=0;
     for (i=0; i<1025; i++) {
	dout[i][j]=dout[i][j]-bg[i];
	if (dout[i][j]<0) {
	   dout[i][j]=0;
	}
	sum+=dout[i][j];
     }
     sum=sum/1025;

     for (i=0; i<1025; i++) {
        dout[i][j]=dout[i][j]-sum;
        if (dout[i][j]<0) {
           dout[i][j]=0;
        }
     }

     sum=0;
     for (i=0; i<1025; i++) { 
        dout[i][j]=dout[i][j]*dout[i][j];
        dout[i][j]=dout[i][j]*dout[i][j];
	sum+=dout[i][j];
     }
     sum=sum/1025; 
 
     for (i=0; i<1025; i++) { 
        dout[i][j]=dout[i][j]-sum;     
        if (dout[i][j]<0) {    
           dout[i][j]=0;       
        }      
     } 

     sum=0; 
     for (i=0; i<1025; i++) { 
	sum+=dout[i][j];
     } 
     sum=sum/1025;

     for (i=0; i<1025; i++) { 
        dout[i][j]=dout[i][j]-sum;     
        if (dout[i][j]<0) {    
           dout[i][j]=0;       
        }      
	dout[i][j]=sqrt ((double)dout[i][j]);
	dout[i][j]=sqrt ((double)dout[i][j]);
     } 
	
  }

  fp2=fopen ("bg.out", "w");
  for (i=0; i<1025; i++) {
	fprintf (fp2, "%d\n", bg[i]);
  }
  fclose (fp2);

  fp2=fopen ("9294all.out", "w");
  for (i=0; i<1025; i++) {
     for (j=0; j<29; j++) {
	fprintf (fp2, "%d ", dout[i][j]);
     }
     fprintf (fp2, "%d\n", dout[i][29]);
  }
  fclose (fp2);

  for (i=0; i<30; i++) {
     if (i<10) {
       j=i+48; 
       rec_file[7]=j;
     }
     if (i>9) {
       rec_file[6]='1';
        j=i-10+48;
       rec_file[7]=j;
     }
     if (i>19) {
       rec_file[6]='2';
        j=i-20+48;
       rec_file[7]=j;
     }
     fp2=fopen (rec_file,"w");
     for (j=0; j<1025; j++)
	fprintf (fp2, "%d\n",dout[j][i]);
     fclose (fp2);
  }

/* Calculate SNR */
/* Original spectra SNR */
  printf ("\n*** Signal and noise ratio (SNR): ***\n\n");

  pn=0;
  ps=0;
  for (j=0; j<30; j++) {
     for (i=0; i<428; i++) {
        pn+=data[i][j]*data[i][j];
     }
     for (i=433; i<508; i++) {
        pn+=data[i][j]*data[i][j];
     }
     for (i=511; i<571; i++) {
        pn+=data[i][j]*data[i][j];
     }
     for (i=574; i<1025; i++) {
        pn+=data[i][j]*data[i][j];
     }
  }
  pn=pn/(1014*30);

  for (i=428; i<433; i++) {
     for (j=0; j<30; j++) {
        ps+=data[i][j]*data[i][j];
     }
  }

  for (i=508; i<511; i++) {
     for (j=0; j<30; j++) {
        ps+=data[i][j]*data[i][j];
     }
  }

  for (i=571; i<574; i++) {
     for (j=0; j<30; j++) {
        ps+=data[i][j]*data[i][j];
     }
  }    
  ps=ps/(11*30);
 
  snr=10*(log10(ps/pn));
  printf ("      Original spetra:   SNR = %f  dB\n", snr);


/* Processed spectra SNR */

  pn=0;
  ps=0;
  for (j=0; j<30; j++) {
     for (i=0; i<428; i++) {
        pn+=dout[i][j]*dout[i][j];
     }
     for (i=433; i<508; i++) {
        pn+=dout[i][j]*dout[i][j];
     }
     for (i=511; i<571; i++) {
        pn+=dout[i][j]*dout[i][j];
     }
     for (i=574; i<1025; i++) {
        pn+=dout[i][j]*dout[i][j];
     }
  }
  pn=pn/(1014*30);

  for (i=428; i<433; i++) {
     for (j=0; j<30; j++) {
        ps+=dout[i][j]*dout[i][j];
     }
  }

  for (i=508; i<511; i++) {
     for (j=0; j<30; j++) {
        ps+=dout[i][j]*dout[i][j];
     }
  }

  for (i=571; i<574; i++) {
     for (j=0; j<30; j++) {
        ps+=dout[i][j]*dout[i][j];
     }
  }    
  ps=ps/(11*30);
  snr=10*(log10(ps/pn));
 
  printf ("      Processed spetra:  SNR = %f  dB\n\n", snr);
}