// system include files
//
#include <stdio.h>
#include <string.h>
#include <math.h>

// isip include files
//
#include "integral.h"

// constants
//
#define NORM  3.0518e-05
#define FRAME_DUR 0.020
#define SAMPLE_FREQ 8000
#define PI  3.14
#define MAX_DUR 60
#define NOMINAL_SIGNAL_NOISE 0.85
#define NOMINAL_NOISE 0.15
#define NUM_BIN 5000

main(){
 
  // declare local variables
  //
  int_4 window_dur_d = (int_4)(0.03*SAMPLE_FREQ);
  int_4 frame_dur_d  = (int_4)(0.02*SAMPLE_FREQ);
  int_4 in_bytes = int_4(FRAME_DUR * SAMPLE_FREQ * sizeof(int_2));
  int_2* signal_buf = new int_2[in_bytes];
  float_4* egy_variance = new float_4[ (int_4)(MAX_DUR/FRAME_DUR) ];
  float_4* sig_buf_d = (float_4*)NULL;
  sig_buf_d = new float_4[window_dur_d];
  memset(signal_buf, (int_4)0, in_bytes);
  int_4 *count;     
  count = new int_4[0.99/0.00001];
  int_4* cdf = new int_4[NUM_BIN];
  int_4 frame_num = 0;
  
  while (fread(signal_buf, 1, in_bytes,stdin) > 0) 
    { 
      // declare local variables
      //
      int_2* buf = (int_2*)signal_buf;
      float_4* window_ptr;
      float_4* signal_ptr;
      float_4 xn;
      float_4 xn_minus_1 = 0;
      
      // find the offset in sig_buf_d
      //
      if (window_dur_d > frame_dur_d)
	{window_ptr = sig_buf_d + (window_dur_d - frame_dur_d);}
      else
	{window_ptr = sig_buf_d;}
      
      // load the signal buffer
      //
      buf++;
      
      for (int i=0; i < frame_dur_d; i++) {
    
	// get a sample
	//
	xn = (float_4)(*buf) * NORM;
	
	// preemphasize
	//
	*window_ptr = xn - 0.95 * xn_minus_1;
	
	// store the previous sample
	//
	xn_minus_1 = xn;
	
	// increment counters
	//
	buf++ ;
	window_ptr++;
      }
      
      // compute energy after hamming windowing the data
      //
      float_4 energy_var = (float_4)0.0;
      float_4 energy = 0;
      float_4 mean = 0;
      signal_ptr = sig_buf_d + window_dur_d - 1;
      window_ptr = sig_buf_d + window_dur_d - 1;
      for (int j=0; j<window_dur_d; j++) {
        
	*signal_ptr=(*window_ptr)*(0.54 - 0.46*(float_4)
				   cos(2.0*PI*(int_4)j/window_dur_d));
	mean += *signal_ptr;	
	energy += (*signal_ptr) * (*signal_ptr);
        signal_ptr--;
	window_ptr--;
      }

      //pass energy value to the energy buffer
      //
      energy_var = energy ;//- mean*mean/window_dur_d;
      
      
      // shift the signal
      //
      if (window_dur_d > frame_dur_d) {
	int_4 shift = window_dur_d - frame_dur_d;
	window_ptr = sig_buf_d + frame_dur_d;
	memcpy(sig_buf_d, window_ptr, shift*sizeof(float_4));
      }
      
      *egy_variance = energy_var;
      int_4 kk = (int_4)((*egy_variance)/0.00001);
      // count[kk]++;
      for(int k = kk; k < NUM_BIN + 1;k++)
	{
	  cdf[k]++;
	}
      egy_variance++;      
      frame_num++;
    }
  
  //define thresholds in cdf to get nominal signal+noise level
  //and nominal noise level
  int_4 kk =(int_4)(cdf[NUM_BIN - 1] * NOMINAL_SIGNAL_NOISE);
  int_4 jj =(int_4)(cdf[NUM_BIN - 1] * NOMINAL_NOISE);
  float_4 signal_noise_level = 0;
  float_4 noise_level = 0;
  int flag = 0;
  for(int k =1;k < NUM_BIN; k++)
    {
      if(cdf[k] > jj && flag == 0 )
	{
	  noise_level = (float_4)(k*0.00001);
	  flag = 1;
	}
      
      if(cdf[k] > kk && flag == 1)
	{
	  signal_noise_level = (float_4)(k*0.00001);
	  break;
	}
    }
  
  // output the calculated snr value
  //   
  float_4 snr = 10*log10((signal_noise_level-noise_level)/noise_level);
  printf("\n %s %f %s\n\n","the snr of the input signal is ", snr," dB");     
  
}