/**************************************************************************
*  echogen.c - C program that simulates echo path for generating echo and
*              mixing with near-end speech
*
***************************************************************************
*  System simulation configuration:
*
*                  ______________
*      farEndIn    |            |      speakerOut
*    ------------->|            |----------------
*                  |            |           ____|____
*                  |            |           |       |
*                  | Speakphone |           |  A(z) |
*                  |   System   |           |_______|
*                  |            |               | acousticEcho 
*                  |            |               |
*                  |            |              +|  +   nearEndIn
*                  |            |<------------(sum)<------------
*                  |____________| microphoneIn
*
*     farEndIn is a far-end speech from a data file
*     nearEndIn is a near-end speech from a data file
*     acousticEcho is acoustic echo generated A(z)
*     microphoneIn is the near-end speech plus acoustic echo
*     A(z) is an acoustic echo path
*
*--------------------------------------------------------------------------
*  Synopsis
*
*    This program is executed by keying in the following command:
*
*        echogen NEARFILE FARFILE MICFILE <enter>
*
*    where
*      NEARFILE  - file name of near-end speech
*      FARFILE  - file name of far-end speech
*      MICFILE - file name of microphoneIn signal
*
*--------------------------------------------------------------------------
*  Functions called by this program
*
*    fir, shift, uran, echopath
*
***************************************************************************
*    Specify auxiliary routines and parameter values
**************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define MAX  1024             // maximum order of arrays

void main(int argc, char *argv[])
{
  void shift();               // function to update signal vector
  void echopath();            // function to generate echo path
  float fir();                // function to perform FIR filtering
  float uran();               // function to generate white noise

/**************************************************************************
*    Define variable arrays, define and clear variables
**************************************************************************/

  int i;                      // integers for indexes
  int AzOrder = 512;          // order of A(z)

  float farEndIn = 0.;        // new data from the far-end
  float nearEndIn = 0.;       // new data from the near-end talker
  float acousticEcho = 0.;    // new data of acoustic echo
  float speakerOut = 0.;      // signal send to loudspeaker
  float microphoneIn = 0.;    // signal pick up by microphone

  float AzCoef[MAX];          // coefficient vector of A(z)
  float AzBuf[MAX];           // signal buffer for A(z) 
 
/**************************************************************************
*   Define variables for generating line and acoustic impulse responses
**************************************************************************/

  int acousticDelay = 10;            // acoustic echo path delay
  float acousticVar = (float)0.054;  // variance
  float acousticDecay = (float)0.012;// decay of acoustic echo

/**************************************************************************
*    Declare file pointers
**************************************************************************/

  FILE *farin;                // file pointer of far_end signal
  FILE *nearin;               // file pointer of near_end signal
  FILE *micout;               // file pointer of transmit signal

/**************************************************************************
*    Define usage of executable file
**************************************************************************/

  if ( argc < 3 )
  {
    fprintf(stderr,"%s\n*** Usage ***\n",argv[0]);
    fprintf(stderr,"echopath NEARFILE FARFILE MICFILE\n");
    fprintf(stderr,"    NEARFILE is file name of near-end speech\n");
    fprintf(stderr,"    FARFILE is file name of far-end speech\n");
    fprintf(stderr,"    MICFILE is file name of microphoneIn signal\n");
    exit(0);
  }

/**************************************************************************
*    Input data and parameters at execution time
**************************************************************************/

  if (( nearin = fopen(argv[1],"r")) == NULL ) // open near-end speech file
  {
    fprintf(stderr,"Can't open the near-end data file %s!\n",argv[1]);
    exit(0);
  }
  if (( farin = fopen(argv[2],"r")) == NULL )  // open far-end speech file
  {
    fprintf(stderr,"Can't open the far-end data file %s!\n",argv[2]);
    exit(0);
  }
  micout = fopen(argv[3],"w");    // open MIC signal data file for storage


/**************************************************************************
*    Clear arrays
**************************************************************************/

  for (i=0; i<AzOrder; ++i)
  {
    AzCoef[i] = 0.;           // clear A(z) coefficient vector
    AzBuf[i] = 0.;            // clear A(z) signal buffer
  }

/**************************************************************************
*    Generate acoustic echo path model for simulation
**************************************************************************/

  echopath(AzCoef, acousticDelay, AzOrder, acousticDecay, acousticVar);

/**************************************************************************
*    Start of main program
***************************************************************************
*    Generate microphoneIn for speakerphone simulation
**************************************************************************/

  while( (fscanf(farin,"%f",&farEndIn)) != EOF)  // read farEndIn data file
  {
    fscanf(nearin,"%f",&nearEndIn);    // read nearEndIn from data file
    speakerOut = farEndIn;             // insert RX_processing here later
 	
    shift(AzBuf, AzOrder, speakerOut); // update signal buffer of A(z)
    acousticEcho = fir(AzBuf,AzCoef,AzOrder);// compute acoustic echo
    microphoneIn = nearEndIn + acousticEcho;
	                      // MIC signal is near-end speech + acoustic echo

    
    fprintf(micout,"%d\n",(int)(microphoneIn+0.5)); // write farEndOut to file

  }
  fcloseall();
}

/**************************************************************************
*  ECHOPATH - function for generating echo path impulse response
***************************************************************************
*  Functions called by this program - uran
*  Algorithm:
*    h(n) = 0,            for 0 <= n < D
*         = v(n)*exp(-c*n), for D <= n < echo_length
*  where
*    v(n) is a white noise with variance = var_v
*    D is echo path delay
*    echo_length is length of echo tail
*    c << 1 is a small constant that control decay rate of echo
*************************************************************************/
void echopath(h, D, echoLength, c, varV)

float *h;
int D, echoLength;
float c, varV;
{
  float uran();                 // random noise generator
  int i;                        // integer for index
  float var;                    // variance of noise
  float uu = 0.;
  var = (float)sqrt(varV);
  for (i=0; i<D; i++)
  {
     h[i] = 0.;                 // h(n) = 0, for 0 <= n < D
  }     
  for (i=D; i<echoLength; i++)  // for D <= n < echo_length
  {
    uu = (float)(var*(uran() - 0.5));
    h[i] = (float)(uu * exp(-c*i)); // h(n) = v(n)*exp(-c*n) 
  }
}

/**************************************************************************
*  URAN - This function generates pseudo-random numbers
**************************************************************************/
static  long  n = (long)12357;// seed I(0) = 12357
float uran()

{
  float ran;                  // random noise r(n)
  n = (long)2045 * n + 1L;    // I(n) = 2045 * I(n-1) + 1
  n -= (n / 1048576L) * 1048576L;
                              // I(n) = I(n) - INT[I(n)/1048576] * 1048576
  ran = (float)(n + 1L) / (float)1048577;
                              // r(n) = FLOAT[I(n) + 1] / 1048577
  return(ran);                // return r(n) to main function
}

/**************************************************************************
*  SHIFT - This function updates signal vector
*          data stored as [x(n)  x(n-1) ... x(n-N+1)]
**************************************************************************/
void  shift(x, N, new)

float *x;
int N;
float new;
{
  int i;                      // index
  for (i=N-1; i>0; --i)
  {
    x[i] = x[i-1];            // shift old data x(n-i), i = 1, 2, ... N-1
  }
  x[0] = (float)new;          // insert new data x(n)
}

/**************************************************************************
*  FIR - This function performs FIR filtering
*          ntap-1
*    y(n) = sum  wi * x(n-i)
*           i=0
**************************************************************************/
float fir(x, w, ntap)

float *x, *w;
int ntap;
{
  float yn;                    // output of FIR filter
  int i;                       // index
  yn = 0.0;                    // y(n) = 0.
  for (i=0; i<ntap; ++i)
  {
    yn += w[i]*x[i];           // FIR filtering of x(n) to get y(n)
  }
  return(yn);                  // return y(n) to main function
}