/* This program calculates the SNR of a speech file */
/* execution:  cat filename | snr.exe               */
#include     <math.h>
#include     <stdio.h>
#include     <math/integral/integral.h>

#define      sampling_freq    8000
#define      window_time      .03
#define      frame_time       .02
#define      PI               3.14159265359
#define      BINS             32767
#define      BYTES            2
#define      ALPHA            0.54
#define      NOM_SIGNAL_LEVEL 0.85
#define      NOM_NOISE_LEVEL  0.15

main()

{
  int_4      window_duration = (int)rint(sampling_freq * window_time);
  //  int_4      window_duration = 5;
  int_4      frame_duration  = (int)rint(sampling_freq * frame_time);
  //  int_4      frame_duration  = 3;
  int_4      overlap = (window_duration - frame_duration)/2;
  float_4    window[window_duration]; 
  float_4    window_buffer[window_duration];
  int_2      sample_buffer[frame_duration]; 

  float_4    alpha = 0.95;
  float_4    xm1 = 0.0;
                       
  float_4    energy = 0.0;
  float_4    rms; 
  int_4      i, Es = 0, En = 0;
  int_4      frames = 0;  
  int_4      pdf[BINS]; 
  float_4    cdf[BINS]; 
  float_4    SNR;

  float_4     xn;

  for (i=0;i<window_duration;i++){
    window[i] = 0;
  }
  for (i=0;i<BINS;i++){
    pdf[i] = 0;
    cdf[i] = 0;
  }


  // pre-fetch the data
  //
  fread(sample_buffer,BYTES,overlap,stdin);

  for (i=0;i<overlap;i++){
    xn = (float_4)sample_buffer[i];
    window_buffer[overlap+i] = xn - alpha * xm1;
    xm1 = xn;
  }
  
  while (fread(sample_buffer,BYTES,frame_duration,stdin) > 0){

    int offset = 2*overlap;
    for (i=0;i<frame_duration;i++){
      xn = (float_4)sample_buffer[i];
      window_buffer[offset+i] = xn - alpha * xm1;
      xm1 = xn;
    }

    //for (i=0; i<window_duration; i++) {
    //  printf("%d: %f\n",i,window_buffer[i]);
    //}
    //printf("... end\n");

    energy = 0.0;
    for (i=0;i<window_duration;i++){
      float val = window_buffer[i]*(ALPHA - .46*cos(2*PI*i/window_duration));
      window[i] = val*val;
      energy += window[i];
    }

    // compute rms
    //
    rms = sqrt(energy/window_duration);
    //    printf("%ld %f\n",frames,rms);
    
    if ((int_4)rms < BINS){
      pdf[(int_4)rint(rms)] += 1;
    }
    frames++;    
    
    // shift window
    //
    int num_to_cpy = window_duration - frame_duration;
    for (i=0;i<num_to_cpy;i++){
      window_buffer[i] = window_buffer[frame_duration + i];
    }
  }
  
  cdf[0] = 0;
  for (i=1;i<BINS;i++){
    cdf[i] = (float)pdf[i]/(float)frames + cdf[i-1];
  }

  //  int nominal_signal_level = cdf[BINS-1]*NOM_SIGNAL_LEVEL;
  //  int nominal_noise_level  = cdf[BINS-1]*NOM_NOISE_LEVEL;

  for (i=1;i<BINS;i++){
    if (cdf[i] > NOM_NOISE_LEVEL) {
      En = i;
      break;
    }
  }
  for (i=1;i<BINS;i++){
    if (cdf[i] > NOM_SIGNAL_LEVEL) {
      Es = i;
      break;
    }
  }
  //  printf("En, Es = %d, %d\n",En,Es);
  SNR = 20.0 * log10((float)(Es - En)/(float)En);
  printf("\nThe SNR is %2.2f dB\n\n", SNR);
}