14e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi#include <math.h>
24e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi#include <stdlib.h>
34e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi#include "fec.h"
44e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi
54e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi#define	MAX_RANDOM	0x7fffffff
64e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi
74e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi/* Generate gaussian random double with specified mean and std_dev */
84e213d510f437769f8a28578dd4f786fb7d16c4Bill Yidouble normal_rand(double mean, double std_dev)
94e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi{
104e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi  double fac,rsq,v1,v2;
114e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi  static double gset;
124e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi  static int iset;
134e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi
144e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi  if(iset){
154e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi    /* Already got one */
164e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi    iset = 0;
174e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi    return mean + std_dev*gset;
184e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi  }
194e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi  /* Generate two evenly distributed numbers between -1 and +1
204e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi   * that are inside the unit circle
214e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi   */
224e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi  do {
234e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi    v1 = 2.0 * (double)random() / MAX_RANDOM - 1;
244e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi    v2 = 2.0 * (double)random() / MAX_RANDOM - 1;
254e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi    rsq = v1*v1 + v2*v2;
264e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi  } while(rsq >= 1.0 || rsq == 0.0);
274e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi  fac = sqrt(-2.0*log(rsq)/rsq);
284e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi  gset = v1*fac;
294e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi  iset++;
304e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi  return mean + std_dev*v2*fac;
314e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi}
324e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi
334e213d510f437769f8a28578dd4f786fb7d16c4Bill Yiunsigned char addnoise(int sym,double amp,double gain,double offset,int clip){
344e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi  int sample;
354e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi
364e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi  sample = offset + gain*normal_rand(sym?amp:-amp,1.0);
374e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi  /* Clip to 8-bit offset range */
384e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi  if(sample < 0)
394e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi    sample = 0;
404e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi  else if(sample > clip)
414e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi    sample = clip;
424e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi  return sample;
434e213d510f437769f8a28578dd4f786fb7d16c4Bill Yi}
44