[PD] FWIW: Box Muller Gaussian noise

[Atwood, Robert C]

FWIW here's a C subroutine with a little 'main' routine for testing, I
did a few years ago for the Box Muller method, pretty much directly
based upon what I read in Numerical Recipes. The text of the section in
Num.Recip. Is a very clear explanation of it , I find , your mileage may
vary! 



      1 #include <math.h>
      2 #include <stdio.h>
      3 #include "machine.h" <---- just neds to choose datatype CA_FLOAT
as double or float or other real-type value
                             <----- and define LOG accordingly... Can be
gotten rid of with some small edits.
      4
      5 /*
      6    BOX-MULLER routine for returning a random number normally
distributed
      7    about a mean of zero, with a standard deviation of 1. Adapted
from Numerical Recipes in C by Press and Flannery.
      8 */
      9
     10 CA_FLOAT box_muller ()
     11 {
     12   static int flag = 0;
     13   static CA_FLOAT gset = 0;
     14   CA_FLOAT v1, v2, fac, gdev, r;
     15
     16   if (!flag) {
     17     do {
     18       v1 = 2 * drand48 () - 1;
     19       v2 = 2 * drand48 () - 1;
     20       r = v1 * v1 + v2 * v2;
     21     } while (r >= 1 || r == 0);
     22     fac = sqrt ((-2 * LOG (r) / r));
     23     gset = v1 * fac;
     24     gdev = v2 * fac;
     25     flag = 1;
     26   } else {
     27     gdev = gset;
     28     flag = 0;
     29   }
     30   return (gdev);
     31 }                               /* end of box_muller */
     32
     33 /* return gaussian devates for a given mean and standard
deviation */
     34 /* by calling the box_muller routine. */
     35
     36 CA_FLOAT gaussdev (CA_FLOAT * params)
     37 {
     38   CA_FLOAT stdev;
     39   CA_FLOAT mean;
     40
     41   if (params == NULL || (params + 1) == NULL) {
     42     fprintf (stderr, "ERROR: gaussdev: You messed up, NULL
pointer here.\n");
     43     exit (1);
     44   }
     45
     46   mean = params[0];
     47   stdev = params[1];
     48   return ((stdev * box_muller ()) + mean);
     49 }
     50
     51 #ifdef TEST_BOX_MULLER
     52 void main ()
     53 {
     54   int i;
     55   CA_FLOAT stdev = 2;
     56   CA_FLOAT mean = 10;
     57   CA_FLOAT params[2];
     58
     59   params[0] = mean; 
     60   params[1] = stdev;
     61   for (i = 0; i < 10000; i++) {
     62     printf ("%.5g\n", gaussdev (params));
     63   }
     64
     65 }
     66 #endif
     67
    
> -----Original Message-----
> From: pd-list-bounces at iem.at [mailto:pd-list-bounces at iem.at] 
> On Behalf Of Andy Farnell
> Sent: 16 March 2008 15:31
> To: pd-list at iem.at
> Subject: [PD] Box Muller Gaussian noise
> 
> Could a stats mathematician please help me check this. (attached) 
> 
> I'm following the Box Muller formula for getting  a cheap Gaussian 
> distribution (instead of adding up 12 sources a la central 
> limit method).
> 
> http://www.dspguru.com/howto/tech/wgn.htm
> 
> Does this look right?
> 
> Also, I have no idea how to check the distribution. It sounds 
> the same as
> uniform noise and looks the same in the spectrograph? What gives?
> Do I need to average over a very long time or something to see any
> difference?
> 
> Cheers all,
> 
> Andy
> 
> -- 
> Use the source
>