/**************************************************************************/
/*                       File: spectrum_main.cc                           */
/**************************************************************************/
/* Project: Computing spectrum of a given signal Class: EE 8993 Spring 96 */
/* Author: Neeraj Deshmukh                           Date: March 2, 1996  */
/**************************************************************************/

/*========================================================================
  This file defines the main function for spectral computation of any 
  given signal using FFT and LP analysis techniques
  ========================================================================*/
/*------------------------------------------------------------------------
  system and ISIP include files
  ------------------------------------------------------------------------*/
#include "spectrum_prms.h"

/*------------------------------------------------------------------------
  main function definition
  ------------------------------------------------------------------------*/
main (int_2 argc, String *argv)
{
  /*-----------------------------
    check for correct commandline
    -----------------------------*/
  if (argc != 5)
    {
      fprintf (stderr, "Usage: spectrum.exe <preemphasis (constant)> <window_duration (seconds)> <center_time> <lp_order> < <speech_file> | xmgr -source stdin \n");
      exit (-1);
    }
  
  /*--------------
    set parameters
    --------------*/
  float_4 preemphasis_d      = atof (argv[1]);
  float_4 window_length_d    = atof (argv[2]);
  float_4 win_center_time_d  = atof (argv[3]);
  int_2   lp_order_d         = atoi (argv[4]);

  // free memory
  //
  delete [] argv;

  // parameters related to signal in terms of sampling frequency
  //
  int_2   sample_frequency_d = (int_2) DEFAULT_SAMPLE_FREQUENCY;
  int_2   window_duration_d  = (int_2) (sample_frequency_d * window_length_d);

  /*-----------------------------------------------------------
    read in the signal frame centered around window center time
    -----------------------------------------------------------*/
  // compute the offset in input data and move the input pointer to that location
  //
  int_4 input_offset_d = (int_4) ((win_center_time_d - 0.5 * window_length_d) 
				  * sample_frequency_d);

  if (fseek (stdin, input_offset_d * sizeof (int_2), SEEK_SET) != 0)
    {
      fprintf (stderr, "Error! Unacceptable center time --- beyond duration of data.\n");
      fprintf (stderr, "Please enter a more suitable number for window center time and try again.\n");
      exit (-1);
    }
  
  // read the signal data
  //
  int_2 *input_frame_d = new int_2[window_duration_d];
  
  if (fread (input_frame_d, sizeof (int_2), window_duration_d, stdin) == 0)
    {
      fprintf (stderr, "Error! Unacceptable center time --- beyond duration of data.\n");
      fprintf (stderr, "Please enter a more suitable number for window center time and try again.\n");
      exit (-1);
    }
  
  /*-----------------
    debias the signal
    -----------------*/
  debias_signal_cc (input_frame_d, window_duration_d);

  /*-----------------------
    preemphasize the signal
    -----------------------*/
  float_4 this_sample_temp  = 0;
  float_4 prev_sample_temp  = 0;
  
  // run through the samples
  //
  float_4 *signal_buffer_d  = new float_4[window_duration_d];
  
  for (int_2 i = 0; i < window_duration_d; i++)
    {
      this_sample_temp    = input_frame_d[i] * DEFAULT_SAMPLE_SCALE_FACTOR;
      signal_buffer_d[i]  = this_sample_temp - preemphasis_d * prev_sample_temp;
      prev_sample_temp    = this_sample_temp;
    }
  
  // free memory
  //
  delete [] input_frame_d;
  
  /*---------------------------------------------------------------
    Hamming window the signal
    ---------------------------------------------------------------*/
  int_2   index_temp        = window_duration_d - 1;
  float_4 hamming_factor_d  = 2.0 * PI / index_temp;
  float_4 ham_temp          = 0;
  
  for (int_2 i = 0; i < window_duration_d / 2; i++)
    {
      ham_temp = 0.54 - 0.46 * cos (hamming_factor_d * i);
      signal_buffer_d[i]              *= ham_temp;
      signal_buffer_d[index_temp - i] *= ham_temp;
    }
  
  /*----------------------------------------------------
    compute the N-point FFT of this frame of signal data
    we use the radix-2 decimation-in-frequency algorithm
    by Burrus & Parks
    ----------------------------------------------------*/
  int_2   fft_order_d     = DEFAULT_FFT_ORDER;
  float_4 freq_scaling_d  = (float_4) sample_frequency_d / fft_order_d;
  
  if (fft_order_d < window_duration_d)
    {
      fprintf (stderr, "The order of FFT is not suitable for this window length.\n");
      fprintf (stderr, "Please try again with a suitable window length.\n");
      exit (-1);
    }
  
  // generate the twiddle factors for the FFT
  //
  float_4 *twiddle_real_d = new float_4[fft_order_d];
  float_4 *twiddle_imag_d = new float_4[fft_order_d];
  compute_twiddle_factors_cc (twiddle_real_d, twiddle_imag_d, fft_order_d);
  
  // initialize the FFT points with signal values and pad with zeros
  //
  float_4 *spectrum_d   = new float_4[fft_order_d];
  for (int_2 i = 0; i < window_duration_d; i++)
    spectrum_d[i] = signal_buffer_d[i];
  for (int_2 i = window_duration_d; i < fft_order_d; i++)
    spectrum_d[i] = 0;
  
  // compute the FFT magnitude spectrum in dB and print it
  //
  compute_fft_cc (spectrum_d, twiddle_real_d, twiddle_imag_d, fft_order_d);
  
  for (int_2 i = 0; i < fft_order_d / 2 + 1; i++)
    fprintf (stdout, "%.4f  %.4f\n", i * freq_scaling_d, spectrum_d[i]);
  
  /*--------------------------------------------------------
    compute the LP coefficients of this frame of signal data
    we use the Levinson-Durbin recursion algorithm
    --------------------------------------------------------*/
  // generate the autocorrelation terms for LP
  //
  float_4 *auto_corr_d = new float_4[lp_order_d + 1];
  for (int_2 i = 0; i < lp_order_d + 1; i++)
    {
      auto_corr_d[i] = 0;
      for (int_2 j = 0; j < window_duration_d - i; j++)
	auto_corr_d[i] += signal_buffer_d[j] * signal_buffer_d[j + i];
    }

  // free memory
  //
  delete [] signal_buffer_d;

  // compute the predictor coefficients for LP
  //
  float_4 *lp_predictor_d = new float_4[lp_order_d + 1];
  for (int_2 i = 0; i < lp_order_d + 1; i++)
    lp_predictor_d[i] = 0;

  float_4 lp_gain_d = compute_lp_cc (lp_predictor_d, lp_order_d, auto_corr_d);

  // free memory
  //
  delete [] auto_corr_d;
  
  // now compute values to plot the impulse response of G/A(z)
  // we use same order as the FFT for sampling the frequency domain for 
  // convenient point-to-point match
  //
  for (int_2 i = 0; i < fft_order_d; i++)
    spectrum_d[i] = 0;

  lp_response_cc (lp_predictor_d, twiddle_real_d, twiddle_imag_d, spectrum_d, 
		  lp_order_d, fft_order_d);
  
  /*----------------
  // alternatively, try computing the FFT of the LP-derived signal
  //
  spectrum_d[0] = 1.0;
  for (int_2 i = 1; i < lp_order_d + 1; i++)
    spectrum_d[i] = -1.0 * lp_predictor_d[i];
  for (int_2 i = lp_order_d + 1; i < fft_order_d; i++)
    spectrum_d[i] = 0;
  
  compute_fft_cc (spectrum_d, twiddle_real_d, twiddle_imag_d, fft_order_d);
  -------------------*/
  
  // now output this to plot
  //
  fprintf (stdout, "\n");
  for (int_2 i = 0; i < fft_order_d / 2 + 1; i++)
    fprintf (stdout, "%.4f  %.4f\n", i * freq_scaling_d, lp_gain_d - spectrum_d[i]);

  // free memory and exit
  //
  delete [] spectrum_d;
  delete [] twiddle_real_d;
  delete [] twiddle_imag_d;
  delete [] lp_predictor_d;
}