/**********************************************************************************************************
/*  Created by Abdel Rabassa (AHR)                                                                        *
/*  This subroutine calculates the frame size for Speex (Quality 8 and Quality 2) encoding algorithm.     *
/*  The computation is based in the TES model, a model that generates data fitting both the marginal      * 
/*  distribution and the autocorrelation function of empirical data simultaneously. The first and second  *
/*  order statistical properties of the data generated follows that of original Speex data analyzed.      *
/*  Note: When implemented, variable Speex_Un must be have global scope.                                  *
/*                                                                                                        *
/*********************************************************************************************************/



#include <iostream>  
#include <cstdlib>
#include <ctime> 
using namespace std;

//************************************************************/
/*   Generate a random number between 0 and 1                */
/*   return a uniform number in [0,1].                       */
/*************************************************************/
 double unifRand()
{
    return rand() / double(RAND_MAX);
}

//************************************************************/
/** Generates a random number bounded by L and R         	**/
/*************************************************************/
double AHR_random_bounded(double L,double R)
{
	double r1;
	if (L>R)   //L<->R
	 {
		r1 = L;
		L = R;
		R = r1;
	 }//L always contains the smallest boundary 
	r1 = (R-L)* unifRand() + L; //generates random number between [L, R)
	return r1;
}

/************************************************************/
/** Receives a random number U[0,1] with autocorrelation = 0**/
/** and returns a Vn random number that follows the         **/   
/** innovation PDF to be used to generate a random number   **/
/** U[0,1] with certain autocorrelation function size that  **/ 
/** follows the PDF                                          */
/*************************************************************/
double AHR_Speex8_generates_Vn (double U)
		{

		double Vn;	
		/* INNOVATION DENSITY STEP FUNCTION FOR SPEEX Q8   (18 steps)
           	k	Lk      	Rk	  Pk
            1	-0.385,	-0.375	0.006427
	        2	-0.375,	-0.365	0.006427
	        3	-0.245,	-0.235	0.00153809
            4	-0.235,	-0.225	0.00175781
            5	-0.225,	-0.215	0.00175781
            6	-0.215,	-0.205	0.00175781
            7	-0.205,	-0.195	0.00175781
            8	-0.195,	-0.185	0.00175781
            9	-0.185,	-0.175	0.00175781
            10	-0.175,	-0.165	0.00230713
            11	-0.165,	-0.155	0.00324097
            12	-0.155,	-0.145	0.00324097
            13	-0.145,	-0.135	0.00324097
            14	-0.135,	-0.125	0.00324097
            15	-0.125,	-0.115	0.00274658
            16	-0.0355556,	-0.0177778	0.878826
            17	 0.145,	0.155	0.000439453
            18	0.164444,	0.171111	0.0777778  */  

            //Generates random numbers in the range [Ln,Rn) with probabilty Pk
			if (U<0.006427)                                         //k = 1
		     {
		      Vn = AHR_random_bounded(-0.385,	-0.375);  
		     }
	        else if (U>= 0.006427       && U <   0.012854     )     //k = 2
		     {
		      Vn = AHR_random_bounded(-0.375,	-0.365);
		     }
	        else if (U>= 0.012854      && U <  0.014392      )     //k = 3
		     {
		      Vn = AHR_random_bounded(-0.245,	-0.235);
		     }
	        else if (U>= 0.014392      && U <   0.01615     )     //k = 4
		     {
		      Vn = AHR_random_bounded(-0.235,	-0.225);
		     }
	        else if (U>=   0.01615    && U <   0.017908     )     //k = 5
		     {
		      Vn = AHR_random_bounded(-0.225,	-0.215);
		     }
	        else if (U>=   0.017908    && U <   0.019666     )     //k = 6
		     {
		      Vn = AHR_random_bounded(-0.215,	-0.205);
		     }
	        else if (U>=   0.019666    && U <    0.021423    )     //k = 7
		     {
		      Vn = AHR_random_bounded(-0.205,	-0.195);
		     }
	        else if (U>=   0.021423    && U <   0.023181     )     //k = 8
		     {
		      Vn = AHR_random_bounded(-0.195,	-0.185);
		     }
	        else if (U>=   0.023181    && U <   0.024939     )     //k = 9
		     {
		      Vn = AHR_random_bounded(-0.185,	-0.175);
		     }
	        else if (U>=   0.024939    && U <   0.027246     )     //k = 10
		     {
		      Vn = AHR_random_bounded(-0.175,	-0.165);
		     }
	        else if (U>=   0.027246    && U <   0.030487     )     //k = 11
		     {
		      Vn = AHR_random_bounded(-0.165,	-0.155);
		     }
	        else if (U>=   0.030487    && U <   0.033728     )     //k = 12
		     {
		      Vn = AHR_random_bounded(-0.155,	-0.145);
		     }
	        else if (U>=  0.033728     && U <     0.036969   )     //k = 13
		     {
		      Vn = AHR_random_bounded(-0.145,	-0.135);
		     }
	        else if (U>=   0.036969    && U <   0.04021     )     //k = 14
		     {
		      Vn = AHR_random_bounded(-0.135,	-0.125);
		     }
	        else if (U>=   0.04021    && U <    0.042957    )     //k = 15
		     {
		      Vn = AHR_random_bounded(-0.125,	-0.115);
		     }
	        else if (U>=    0.042957   && U <   0.921783     )     //k = 16
		     {
		      Vn = AHR_random_bounded(-0.0355556,	-0.0177778);
		     }
	        else if (U>=   0.921783    && U <   0.922222     )     //k = 17
		     {
		      Vn = AHR_random_bounded( 0.145,	0.155);
		     }
	        else                                                     //k = 18
		     {
		      Vn = AHR_random_bounded(0.164444,	0.171111);
		     }

		return(Vn);
		
		
		}

//************************************************************/
//**                AHR CREATED                             **/
/** Receives a random number U[0,1] with autocorrelation    **/
/** determined by the Innovation function and returns a frame*/
/** size that follows the PDF of Speex Q8 original  data      */
/*************************************************************/
double AHR_Speex8_histogram_inversion(double U)
{
  double f_s;         
   /*    SpeexQ8 Histogram		
   Cell Left boundary	Right boundary	Probability
   	1	4.2	    6.848	0.000335627
   	2	17.44	20.088	0.000167813
   	3	20.088	22.736	0.000335627
   	4	28.032	30.68	0.00050344
   	5	30.68	33.328	0.0062091
   	6	33.328	35.976	0.0241651
   	7	35.976	38.624	0.0345696
   	8	38.624	41.272	0.0689713
   	9	41.272	43.92	0.0823964
   	10	43.92	46.568	0.145662
   	11	46.568	49.216	0.0954858
   	12	49.216	51.864	0.140292
   	13	51.864	54.512	0.194328
   	14	54.512	57.16	0.00922974
   	15	57.16	59.808	0.0768585
   	16	59.808	62.456	0.0583991
   	17	62.456	65.104	0.0355764
   	18	65.104	67.752	0.0248364
   	19	67.752	70.4	0.00167813   */
   if (U<0.000335627)                                            //Cell = 1
     {
      f_s = AHR_random_bounded(4.2,  6.848);  
     }
    else if (U>= 0.000335627       && U <   0.00050344     )     //Cell = 2
     {
      f_s = AHR_random_bounded(17.44,  20.088); 
     }
    else if (U>= 0.00050344      && U <  0.000839067      )     //Cell = 3
     {
      f_s = AHR_random_bounded(20.088,  22.736); 
     }
    else if (U>= 0.000839067      && U <   0.001342507    )     //Cell = 4
     {
      f_s = AHR_random_bounded(28.032,  30.68); 
     }
    else if (U>=   0.001342507    && U <   0.007551607    )     //Cell = 5
     {
      f_s = AHR_random_bounded(30.68,  33.328); 
     }
    else if (U>=   0.007551607   && U <   0.031716707     )     //Cell = 6
     {
      f_s = AHR_random_bounded(33.328,  35.976); 
     }
    else if (U>=   0.031716707    && U <   0.066286307   )     //Cell = 7
     {
      f_s = AHR_random_bounded(35.976,  38.624); 
     }
    else if (U>=   0.066286307   && U <   0.135257607    )     //Cell = 8
     {
      f_s = AHR_random_bounded(38.624,  41.272); 
     }
    else if (U>=   0.135257607    && U <   0.217654007    )     //Cell = 9
     {
      f_s = AHR_random_bounded(41.272,  43.92); 
     }
    else if (U>=   0.217654007    && U <  0.363316007   )     //Cell = 10
     {
      f_s = AHR_random_bounded(43.92,  46.568); 
     }
    else if (U>=  0.363316007   && U <  0.458801807    )     //Cell = 11
     {
     f_s = AHR_random_bounded(46.568,  49.216); 
     }
    else if (U>=   0.458801807   && U <   0.599093807    )     //Cell = 12
     {
      f_s = AHR_random_bounded(49.216,  51.864); 
     }
    else if (U>=  0.599093807    && U <     0.793421807  )     //Cell = 13
     {
      f_s = AHR_random_bounded(51.864,  54.512); 
     }
    else if (U>=  0.793421807   && U <   0.802651547   )     //Cell = 14
     {
      f_s = AHR_random_bounded(54.512,  57.16); 
     }
    else if (U>=   0.802651547   && U <    0.879510047    )     //Cell = 15
     {
      f_s = AHR_random_bounded(57.16,  59.808); 
     }
    else if (U>=    0.879510047   && U <   0.937909147    )     //Cell = 16
     {
      f_s = AHR_random_bounded(59.808,  62.456); 
     }
    else if (U>=   0.937909147   && U <  0.973485547  )     //Cell = 17
     {
      f_s = AHR_random_bounded(62.456,  65.104); 
     }
    else if (U>=   0.973485547    && U <  0.998321947    )     //Cell = 18
     {
      f_s = AHR_random_bounded(65.104,  67.752); 
     }					
    else                                                     //Cell = 19
     {
      f_s = AHR_random_bounded(67.752,  70.4); 
     }
		
return (f_s);

}






/************************************************************/
//**                AHR CREATED                             **/
/** Receives a random number U[0,1] with autocorrelation = 0**/
/** and returns a Vn random number that follows the         **/   
/** innovation PDF to be used to generate a random number   **/
/** U[0,1] with certain autocorrelation function size that  **/ 
/** follows the PDF                                          */
/*************************************************************/
double AHR_Speex2_generates_Vn (double U)  
{

	double Vn;	
	/*INNOVATION DENSITY STEP FUNCTION FOR SPEEX Q2   (2 steps)
	k   Lk       Rk         Pk
	1 -0.035,  -0.0275,  0.504235
	2  0.075,    0.085,   0.495765  */
	//Generates random numbers in the range [Ln,Rn) with probabilty Pk
	if (U<0.504235)                                         //k = 1
	  {
		Vn = AHR_random_bounded(-0.035,  -0.0275); 
	  }
	else
	  {
		Vn = AHR_random_bounded(0.075,  0.085);
	  }
	return(Vn);
}

//************************************************************/
//**                AHR CREATED                             **/
/** Receives a random number U[0,1] with autocorrelation    **/
/** determined by the Innovation function and returns a frame*/
/** size that follows the PDF of Speex Q2 original  data      */
/*************************************************************/
double AHR_Speex2_histogram_inversion(double U)   
{
  double f_s;         
   /*    Speex Q2 Histogram		
   Cell Left boundary	Right boundary	Probability
    1     0,  0.224,  0.000320461
   	2	1.448,  2.672,  0.000160231		   
   	3	5.12,  6.344,  0.000160231		   
   	4	6.344,  7.568,  0.000480692		   
   	5	8.792,  10.016,  0.004166		   
   	6	10.016,  11.24,  0.00288415		   
   	7	11.24,  12.464,  0.0365326		   
   	8	12.464,  13.688,  0.0169845		   
   	9	13.688,  14.912,  0.0620093		   
   	10	14.912,  16.136,  0.195481		   
   	11	16.136,  17.36,  0.0990226		   
   	12	17.36,  18.584,  0.0722641		   
   	13	18.584,  19.808,  0.118411		   
   	14	19.808,  21.032,  0.0922929		   
   	15	21.032,  22.256,  0.167281		   
   	16	22.256,  23.48,  0.064573		   
   	17	23.48,  24.704,  0.00176254		   
   	18	24.704,  25.928,  0.0368531		   
   	19	25.928,  27.152,  0.0190675		   
   	20	27.152,  28.376,  0.00881269		   
   	21	28.376,  29.6,  0.000480692
   */

        if (U<0.000320461)                                            //Cell = 1
	     {
	      f_s = AHR_random_bounded(0,  0.224);  
	     }
        else if (U>= 0.000320461       && U <   0.000480692     )     //Cell = 2
	     {
	      f_s = AHR_random_bounded(1.448,  2.672); 
	     }
        else if (U>= 0.000480692      && U <  0.000640923      )     //Cell = 3
	     {
	      f_s = AHR_random_bounded(5.12,  6.344); 
	     }
        else if (U>= 0.000640923      && U <   0.001121615    )     //Cell = 4
	     {
	      f_s = AHR_random_bounded(6.344,  7.568); 
	     }
        else if (U>=   0.001121615    && U <   0.005287615    )     //Cell = 5
	     {
	      f_s = AHR_random_bounded(8.792,  10.016); 
	     }
        else if (U>=   0.005287615   && U <   0.008171765     )     //Cell = 6
	     {
	      f_s = AHR_random_bounded(10.016,  11.24); 
	     }
        else if (U>=   0.008171765    && U <   0.044704365   )     //Cell = 7
	     {
	      f_s = AHR_random_bounded(11.24,  12.464); 
	     }
        else if (U>=   0.044704365   && U <   0.061688865    )     //Cell = 8
	     {
	      f_s = AHR_random_bounded(12.464,  13.688); 
	     }
        else if (U>=   0.061688865    && U <   0.123698165    )     //Cell = 9
	     {
	      f_s = AHR_random_bounded(14.912,  16.136); //i swapped the frma cell values for cell 9 and 10
	     }
        else if (U>=   0.123698165    && U <  0.319179165   )     //Cell = 10
	     {
	      f_s = AHR_random_bounded(13.688,  14.912); 
	     }
        else if (U>=  0.319179165   && U <  0.418201765    )     //Cell = 11
	     {
	     f_s = AHR_random_bounded(16.136,  17.36); 
	     }
        else if (U>=   0.418201765   && U <   0.490465865    )     //Cell = 12
	     {
	      f_s = AHR_random_bounded(17.36,  18.584); 
	     }
        else if (U>=  0.490465865    && U <     0.608876865  )     //Cell = 13
	     {
	      f_s = AHR_random_bounded(18.584,  19.808); 
	     }
        else if (U>=  0.608876865   && U <   0.701169765   )     //Cell = 14
	     {
	      f_s = AHR_random_bounded(19.808,  21.032); 
	     }
        else if (U>=   0.701169765   && U <    0.868450765    )     //Cell = 15
	     {
	      f_s = AHR_random_bounded(21.032,  22.256); 
	     }
        else if (U>=    0.868450765   && U <   0.933023765    )     //Cell = 16
	     {
	      f_s = AHR_random_bounded(22.256,  23.48); 
	     }
        else if (U>=   0.933023765   && U <  0.934786305  )     //Cell = 17
	     {
	      f_s = AHR_random_bounded(23.48,  24.704); 
	     }
        else if (U>=   0.934786305    && U <  0.971639405    )     //Cell = 18
	     {
	      f_s = AHR_random_bounded(24.704,  25.928); 
	     }		
		 else if (U>=   0.971639405    && U <  0.990706905    )     //Cell = 19
	     {
	      f_s = AHR_random_bounded(25.928,  27.152); 
	     }					
		  else if (U>=   0.990706905    && U <  0.999519595    )   //Cell = 20
	     {
	      f_s = AHR_random_bounded(27.152,  28.376); 
	     }					
        else                                                     //Cell = 21
	     {
	      f_s = AHR_random_bounded(28.376,  29.6); 
	     }
return (f_s);

}

void main()
{
	
	const double fr_duration = 0.02;          // Frame duration = 20 ms in Speex
	int          Speex_quality;               // Speex Quality mode. Speex Quality =8 and Quality =2 allowed. 
	int          fr_packet;                   // number of voice frmaes to be encapsulated per packet
	int          fr_size;
    
	int          i;
	double       r1;
	double       Speex_Un;
	
	cout<<"Enter Speex quality mode [8:SpeexQ8, 2:SpeexQ2]: ";
	cin>>Speex_quality;
	cout<<"Enter the number of frames per packet: ";
	cin>>fr_packet;
    
	srand((unsigned)time(0)); 
    Speex_Un = (double) rand() / (RAND_MAX+1);
	
	for (i=0; i<fr_packet; i++) 
	{
	   //Calculates 1 frame size
	   
	   r1 = (double) rand() / (RAND_MAX+1);
	  
	   //TES phase 1. Obtain a U[0,1] with autocorrelation that follows innovation function
			
	  //Un = U(n-1) +   Vn       
	   if (Speex_quality ==8) 
	   {
	    Speex_Un = Speex_Un+ AHR_Speex8_generates_Vn (r1);   
	   }
	   else if (Speex_quality ==2)
	   {
        Speex_Un = Speex_Un+ AHR_Speex2_generates_Vn (r1);
	   }

	   //Modulo-1 operation
          Speex_Un = Speex_Un - (int)Speex_Un; 
	      if (Speex_Un<0)
		    Speex_Un = 1+ Speex_Un;
		   
		//TES phase 2. Sends the autocorrelated Un and receive a variable that follows Speex frame size PDF
	   if (Speex_quality ==8) 
	   {
	    fr_size = AHR_Speex8_histogram_inversion( Speex_Un );   
	   }
	   else if (Speex_quality ==2)
	   {
        fr_size = AHR_Speex2_histogram_inversion( Speex_Un );
	   }	
		
	   cout<< "Frame " <<i+1 << " size is: ";
	   cout<< fr_size << " bytes\nFrame duration is: ";
	   cout<<fr_duration*1000<< " ms\n";";
    }
}