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@iem.at [mailto:pd-list-bounces@iem.at] On Behalf Of Andy Farnell Sent: 16 March 2008 15:31 To: pd-list@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