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