tomwalters@0: /*
tomwalters@0:    noise.c   generate samples of a normally distributed deviate (white noise)
tomwalters@0:    -------   in binary shorts or floats.
tomwalters@0: 
tomwalters@0: 
tomwalters@0:    The routine gasdev returns normally distributed numbers with zero mean
tomwalters@0:    and unit variance, using the Box-Muller method described in Numerical
tomwalters@0:    recipes [p203]. It is based upon a random number generator ran1 which
tomwalters@0:    generates uniformly distributed numbers in the range [0,1].
tomwalters@0: 
tomwalters@0:    The period of the random sequences is claimed to be infinite,
tomwalters@0:    [Numerical recipes, p196].
tomwalters@0:    Random sequences which repeat can be heard as periodic if the sequence
tomwalters@0:    repeats every 4secs or less. Therefore the period of a random sequence
tomwalters@0:    should be greater than 5secs, or 100000 samples (at 20kHz samplerate).
tomwalters@0: 
tomwalters@0:    The seed of each random sequence is the variable "seed", which should be
tomwalters@0:    set to any negative value to initialize or re-initialize the sequence,
tomwalters@0:    [Numerical recipes, p196]. The same seed gives the same random sequence.
tomwalters@0:    Otherwise, (and by default when seed=off or seed=0), the system call
tomwalters@0:    getpid() supplies a new seed each run.
tomwalters@0: 
tomwalters@0:    Each sample is scaled using the given mean and variance parameters so that
tomwalters@0:    samples are drawn from a normal distribution of given mean and variance.
tomwalters@0: */
tomwalters@0: 
tomwalters@0: 
tomwalters@0: #include <stdio.h>
tomwalters@0: #include <math.h>
tomwalters@0: #include "options.h"
tomwalters@0: #include "units.h"
tomwalters@0: #include "strmatch.h"
tomwalters@0: 
tomwalters@0: char applic[]  = "generate white noise. " ;
tomwalters@0: char usage[]   = "noise [options]" ;
tomwalters@0: 
tomwalters@0: static char *helpstr, *debugstr, *sampstr, *durstr, *meanstr, *varstr, *seedstr, *datastr ;
tomwalters@0: 
tomwalters@0: static Options option[] = {
tomwalters@0:     {   "help"      ,   "off"       ,  &helpstr     ,   "help"                      , DEBUG   },
tomwalters@0:     {   "debug"     ,   "off"       ,  &debugstr    ,   "debugging switch"          , DEBUG   },
tomwalters@0:     {   "samplerate",   "20kHz"     ,  &sampstr     ,   "samplerate "               , VAL     },
tomwalters@0:     {   "duration"  ,   "500ms"     ,  &durstr      ,   "duration of waveform"      , VAL     },
tomwalters@0:     {   "mean"      ,   "0"         ,  &meanstr     ,   "mean value of waveform"    , VAL     },
tomwalters@0:     {   "variance"  ,   "500"       ,  &varstr      ,   "variance of values"        , VAL     },
tomwalters@0:     {   "seed"      ,   "off"       ,  &seedstr     ,   "seed for random sequence"  , VAL     },
tomwalters@0:     {   "type"      ,   "short"     ,  &datastr     ,   "o/p datatype (short/float)", VAL     },
tomwalters@0:    ( char * ) 0 } ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: int     samplerate ;
tomwalters@0: int     duration   ;
tomwalters@0: float   mean       ;
tomwalters@0: float   variance   ;
tomwalters@0: int     seed       ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: main (argc,argv)
tomwalters@0: int    argc;
tomwalters@0: char **argv;
tomwalters@0: {
tomwalters@0:     float gasdev();
tomwalters@0:     float n2 ;
tomwalters@0:     short n1 ;
tomwalters@0: 
tomwalters@0:     getopts( option,argc,argv ) ;
tomwalters@0:     if ( !isoff( helpstr ) )
tomwalters@0: 	helpopts3( helpstr, argv[0], applic, usage, option ) ;
tomwalters@0: 
tomwalters@0:     samplerate = to_Hz( sampstr, 0 ) ;
tomwalters@0:     duration   = to_p( durstr, samplerate ) ;
tomwalters@0:     mean       = atof( meanstr ) ;
tomwalters@0:     variance   = atof( varstr ) ;
tomwalters@0: 
tomwalters@0:     if ( isoff( seedstr ) )
tomwalters@0: 	seed = ( -getpid() ) ;
tomwalters@0:     else
tomwalters@0: 	seed = atoi( seedstr ) ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     if ( iststr( datastr, "short" ) )
tomwalters@0: 	while ( duration-- ) {
tomwalters@0: 	    n1 = (short)( gasdev( &seed ) * variance + mean ) ;
tomwalters@0: 	    fwrite( &n1, sizeof(short), 1, stdout ) ;
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0:     else if ( iststr( datastr, "float" ) )
tomwalters@0: 	while ( duration-- ) {
tomwalters@0: 	    n2 = gasdev( &seed ) * variance + mean ;
tomwalters@0: 	    fwrite( &n2, sizeof(float), 1, stdout ) ;
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0:     else
tomwalters@0: 	fprintf(stderr,"noise: unknown datatype [%s]\n", datastr) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: float gasdev( idum )
tomwalters@0: int *idum;
tomwalters@0: {
tomwalters@0:     static int   iset=0;
tomwalters@0:     static float gset;
tomwalters@0:     float  fac,r,v1,v2;
tomwalters@0:     float  ran1();
tomwalters@0: 
tomwalters@0:     if ( iset == 0 ) {
tomwalters@0: 	do {
tomwalters@0: 	    v1 = 2.0 * ran1( idum ) - 1.0 ;
tomwalters@0: 	    v2 = 2.0 * ran1( idum ) - 1.0 ;
tomwalters@0: 	    r  = v1*v1+v2*v2;
tomwalters@0: 	} while (r >= 1.0);
tomwalters@0: 	fac  = sqrt(-2.0*log(r)/r);
tomwalters@0: 	gset = v1*fac;
tomwalters@0: 	iset = 1;
tomwalters@0: 	return v2*fac;
tomwalters@0:     }
tomwalters@0:     else {
tomwalters@0: 	iset = 0;
tomwalters@0: 	return gset;
tomwalters@0:     }
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: #define M1 259200
tomwalters@0: #define IA1 7141
tomwalters@0: #define IC1 54773
tomwalters@0: #define RM1 (1.0/M1)
tomwalters@0: #define M2 134456
tomwalters@0: #define IA2 8121
tomwalters@0: #define IC2 28411
tomwalters@0: #define RM2 (1.0/M2)
tomwalters@0: #define M3 243000
tomwalters@0: #define IA3 4561
tomwalters@0: #define IC3 51349
tomwalters@0: 
tomwalters@0: float ran1( idum )
tomwalters@0: int *idum;
tomwalters@0: {
tomwalters@0:     static long ix1,ix2,ix3;
tomwalters@0:     static float r[98];
tomwalters@0:     float temp;
tomwalters@0:     static int iff=0;
tomwalters@0:     int j;
tomwalters@0: 
tomwalters@0:     if ( *idum < 0 || iff == 0 ) {
tomwalters@0: 	iff=1;
tomwalters@0: 	ix1=(IC1-(*idum)) % M1;
tomwalters@0: 	ix1=(IA1*ix1+IC1) % M1;
tomwalters@0: 	ix2=ix1 % M2;
tomwalters@0: 	ix1=(IA1*ix1+IC1) % M1;
tomwalters@0: 	ix3=ix1 % M3;
tomwalters@0: 	for (j=1;j<=97;j++) {
tomwalters@0: 	    ix1=(IA1*ix1+IC1) % M1;
tomwalters@0: 	    ix2=(IA2*ix2+IC2) % M2;
tomwalters@0: 	    r[j]=(ix1+ix2*RM2)*RM1;
tomwalters@0: 	}
tomwalters@0: 	*idum=1;
tomwalters@0:     }
tomwalters@0:     ix1=(IA1*ix1+IC1) % M1;
tomwalters@0:     ix2=(IA2*ix2+IC2) % M2;
tomwalters@0:     ix3=(IA3*ix3+IC3) % M3;
tomwalters@0:     j=1 + ((97*ix3)/M3);
tomwalters@0:     if (j > 97 || j < 1) {
tomwalters@0: 	fprintf(stderr, "noise: This cannot happen\n");
tomwalters@0: 	exit(1);
tomwalters@0:     }
tomwalters@0:     temp=r[j];
tomwalters@0:     r[j]=(ix1+ix2*RM2)*RM1;
tomwalters@0:     return temp;
tomwalters@0: }
tomwalters@0: