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

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

/*
   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;
}