diff tools/noise.c @ 0:5242703e91d3 tip

Initial checkin for AIM92 aimR8.2 (last updated May 1997).

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