1#include <math.h>
2#include <stdlib.h>
3#include "fec.h"
4
5#define	MAX_RANDOM	0x7fffffff
6
7/* Generate gaussian random double with specified mean and std_dev */
8double normal_rand(double mean, double std_dev)
9{
10  double fac,rsq,v1,v2;
11  static double gset;
12  static int iset;
13
14  if(iset){
15    /* Already got one */
16    iset = 0;
17    return mean + std_dev*gset;
18  }
19  /* Generate two evenly distributed numbers between -1 and +1
20   * that are inside the unit circle
21   */
22  do {
23    v1 = 2.0 * (double)random() / MAX_RANDOM - 1;
24    v2 = 2.0 * (double)random() / MAX_RANDOM - 1;
25    rsq = v1*v1 + v2*v2;
26  } while(rsq >= 1.0 || rsq == 0.0);
27  fac = sqrt(-2.0*log(rsq)/rsq);
28  gset = v1*fac;
29  iset++;
30  return mean + std_dev*v2*fac;
31}
32
33unsigned char addnoise(int sym,double amp,double gain,double offset,int clip){
34  int sample;
35
36  sample = offset + gain*normal_rand(sym?amp:-amp,1.0);
37  /* Clip to 8-bit offset range */
38  if(sample < 0)
39    sample = 0;
40  else if(sample > clip)
41    sample = clip;
42  return sample;
43}
44