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
Hi,
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?
I agree that this method is better than adding up 12 uniform distributions, which is not only inefficient, but not that accurate.
Your implementation looks about right, it would get the X but not the Y part of the random normal generation. I am not too sure about the multiplication by 0.5 at the very end, but maybe there is some pd related reason for that. It's not part of the usual formula, but would just have an effect on the scale of the distribution, not the shape.
Also, I have no idea how to check the distribution.
My pd knowledge is a bit weak (still a relative beginner), but if you could sample from that stream of random numbers, save the samples into a file and plot a histogram, you could check for approximate normality that way. There are free programs like R (http://www.r-project.org/) available that draw histograms and are available for all platforms, many others as well.
It sounds the same as uniform noise and looks the same in the spectrograph? What gives?
In my experience, it is extremely hard to distinguish different distributions by ear if they are all used for random noise streams. I teach statistics, and one year, trying to combine my music hobby with my work, tried to see if I could get students to recognize different distributions by their sounds...its very difficult if you just let them generate random noise. Some composers like Xenakis used different distributions to control other parameters such as pitch choice, and there, over time, you might recognize say a uniform distribution of 12 tones from one using an approximate normal distribution discretized to these 12 tones, as some tones might be recognized as appearing more often than others. That's just one simple example of random number to musical parameter function, of course.
Do I need to average over a very long time or something to see any difference?
If you want to SEE the difference, you just need a few values from each distribution, plot histograms, and it will be obvious which is a normal, which is a uniform, or Poisson, etc. However, to hear the difference is a trickier proposition.
I hope this helps, Lawrence
Andy Farnell wrote:
Could a stats mathematician please help me check this. (attached)
Well I'm not one but...
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?
The way I understand it, white noise has the characteristic that if you count the number of occurrences of each possible value they will be normally distributed (the bell curve, with most values being near the peak, since all values are equally probable) with any particular value being located at random on the curve -- there is no correlation between a given value and its probability of occurrence.
Gaussian noise is distributed so that the values themselves follow the bell curve with the extreme values at the tails of the curve and the median values in the middle (on average of course, any given trial can give weird results).
So, to check your distribution you could write some values to a table and then bin them according to value to see if the points around zero occur a lot more often than the points near one and minus one.
I tried your patch with [noise~] and [gaussianoise] writing to tables and visually compared the results and it's true you can't hear the difference but you can definitely see it.
(gaussianoise has occasional values that exceed [-1 ... 1], which I suppose is normal...white noise is always on [-1...1])
With white noise there should be no preferred values but gaussian noise will have a lot more hovering around zero.
Martin
Cheers all,
Andy
PD-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list
On Sun, Mar 16, 2008 at 11:16 AM, Martin Peach martin.peach@sympatico.ca wrote:
(gaussianoise has occasional values that exceed [-1 ... 1], which I suppose is normal...white noise is always on [-1...1])
That's true. With the Box-Muller method, there is the log(~U1) term, but you can always just add a small value to U1, which will truncate your distribution. The size of the small value can be calculated to fit with any given threshold.
With white noise there should be no preferred values but gaussian noise will have a lot more hovering around zero.
Martin
Cheers all,
Andy
PD-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list
PD-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list
Charles Henry wrote:
On Sun, Mar 16, 2008 at 11:16 AM, Martin Peach martin.peach@sympatico.ca wrote:
(gaussianoise has occasional values that exceed [-1 ... 1], which I suppose is normal...white noise is always on [-1...1])
That's true. With the Box-Muller method, there is the log(~U1) term, but you can always just add a small value to U1, which will truncate your distribution. The size of the small value can be calculated to fit with any given threshold.
I think it's really because the Box-Muller method selects random numbers in pairs which map to points in a unit square on the plane, but then selects only those points which are inside the unit circle, something that the pd patch doesn't do (how to resample points in a dsp vector until they are in range?). The attached patch shows the straightforward way of doing it by simply selecting a random radius and angle and returning the resulting y coordinate as the random number. The results are always on [-1,1]. I don't think sin~ will be any slower than log~.
Martin
#N canvas 767 150 450 351 10; #X obj 86 2 noise~; #X obj 135 192 tabwrite~ bbb; #X obj 145 159 table bbb; #X obj 135 137 bng 15 250 50 0 empty empty empty 17 7 0 10 -262144 -1 -1; #X obj 75 258 dac~; #X obj 274 43 noise~; #X obj 274 64 *~ 3.14159; #X obj 275 98 sin~; #X obj 85 115 *~; #X text 151 135 write to table; #X obj 84 193 *~; #X obj 32 30 vsl 15 128 0 1 0 0 empty empty empty 0 -9 0 10 -4034 -13381 -1 2000 1; #X text 89 80 random radius on [-1 , 1]; #X text 278 79 random angle on [-pi , +pi]; #X connect 0 0 8 0; #X connect 3 0 1 0; #X connect 5 0 6 0; #X connect 6 0 7 0; #X connect 7 0 8 1; #X connect 8 0 1 0; #X connect 8 0 10 0; #X connect 10 0 4 0; #X connect 10 0 4 1; #X connect 11 0 10 1;
Oh no that's wrong isn't it :( The log is necessary to keep the distribution normal, and the range is going to get wider the closer to zero the radius is allowed to get. The attached patch has a scale adjustment... Still I wonder what kind of distribution gaussianoise2 gives, it's not just white.
Martin
Martin Peach wrote:
Charles Henry wrote:
On Sun, Mar 16, 2008 at 11:16 AM, Martin Peach martin.peach@sympatico.ca wrote:
(gaussianoise has occasional values that exceed [-1 ... 1], which I suppose is normal...white noise is always on [-1...1])
That's true. With the Box-Muller method, there is the log(~U1) term, but you can always just add a small value to U1, which will truncate your distribution. The size of the small value can be calculated to fit with any given threshold.
I think it's really because the Box-Muller method selects random numbers in pairs which map to points in a unit square on the plane, but then selects only those points which are inside the unit circle, something that the pd patch doesn't do (how to resample points in a dsp vector until they are in range?). The attached patch shows the straightforward way of doing it by simply selecting a random radius and angle and returning the resulting y coordinate as the random number. The results are always on [-1,1]. I don't think sin~ will be any slower than log~.
Martin
PD-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list
#N canvas 767 150 450 461 10; #X obj 65 -59 noise~; #X obj 135 192 tabwrite~ bbb; #X obj 145 159 table bbb; #X obj 135 137 bng 15 250 50 0 empty empty write 17 7 0 10 -4034 -13381 -13381; #X obj 75 258 dac~; #X obj 274 43 noise~; #X obj 274 64 *~ 3.14159; #X obj 275 98 sin~; #X obj 85 115 *~; #X obj 84 193 *~; #X obj 21 34 vsl 15 128 0 1 0 1 empty empty volume 0 -9 0 10 -4034 -13381 -13381 3800 1; #X text 278 79 random angle on [-pi , +pi]; #X obj 64 -38 abs~; #X obj 64 23 ln~; #X obj 64 77 sqrt~; #X obj 64 46 *~ -2; #X obj 64 -18 *~ 0.999; #X obj 65 3 +~ 0.001; #X text 103 77 random radius on [0 , big]; #X obj 241 -14 -; #X obj 230 -62 t b f; #X msg 231 -41 1; #X floatatom 241 8 5 0 0 0 - - -; #X floatatom 230 -83 5 0 0 0 - - -; #X obj 233 -105 hsl 128 15 0 1 0 1 empty empty scale -2 -8 0 10 -4034 -13381 -13381 2000 1; #X obj 135 -101 metro 100; #X obj 135 -121 tgl 15 0 empty empty plot 17 7 0 10 -4034 -13381 -13381 0 1; #X connect 0 0 12 0; #X connect 3 0 1 0; #X connect 5 0 6 0; #X connect 6 0 7 0; #X connect 7 0 8 1; #X connect 8 0 1 0; #X connect 8 0 9 0; #X connect 9 0 4 0; #X connect 9 0 4 1; #X connect 10 0 9 1; #X connect 12 0 16 0; #X connect 13 0 15 0; #X connect 14 0 8 0; #X connect 15 0 14 0; #X connect 16 0 17 0; #X connect 17 0 13 0; #X connect 19 0 22 0; #X connect 19 0 17 1; #X connect 20 0 21 0; #X connect 20 1 19 1; #X connect 21 0 19 0; #X connect 23 0 20 0; #X connect 23 0 16 1; #X connect 24 0 23 0; #X connect 25 0 3 0; #X connect 26 0 25 0;
from the first equation that andy posted, I produced a gem representation. the box muller noise seems wrong, because it does not use the whole range but is shifted to the negative side. note, this is not a distribution of frequencies, but of noise values.. marius.
Martin Peach wrote:
Oh no that's wrong isn't it :( The log is necessary to keep the distribution normal, and the range is going to get wider the closer to zero the radius is allowed to get. The attached patch has a scale adjustment... Still I wonder what kind of distribution gaussianoise2 gives, it's not just white.
Martin
Martin Peach wrote:
Charles Henry wrote:
On Sun, Mar 16, 2008 at 11:16 AM, Martin Peach martin.peach@sympatico.ca wrote:
(gaussianoise has occasional values that exceed [-1 ... 1], which I suppose is normal...white noise is always on [-1...1])
That's true. With the Box-Muller method, there is the log(~U1) term, but you can always just add a small value to U1, which will truncate your distribution. The size of the small value can be calculated to fit with any given threshold.
I think it's really because the Box-Muller method selects random numbers in pairs which map to points in a unit square on the plane, but then selects only those points which are inside the unit circle, something that the pd patch doesn't do (how to resample points in a dsp vector until they are in range?). The attached patch shows the straightforward way of doing it by simply selecting a random radius and angle and returning the resulting y coordinate as the random number. The results are always on [-1,1]. I don't think sin~ will be any slower than log~.
Martin
PD-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list
PD-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list
GEM is broken here, but thanks for the info Marius. I'm reading through the docs for R at the moment. It makes lovely plots, but haven't figured how to get my data in to it yet...
JFYI the application is rainfall. Many papers I read describe rainfall as Gaussian.
I know from physical analysis that raindrops are uniform in size and velocity for any local sample, so I've realised this distribution is about how they fall within an area and pondering how a distribution can be Gaussian in 2D.
Thing is, I can't figure out any good reason why rain should by anything other than uniformly distributed ! :(
When I use Martins second patch with a thresholding function to trigger droplet sounds, it does sound a lot more like real rainfall than a uniformly triggered model.
I'm in one of those grey areas where I half understand what I'm doing, which is a dangerous place to be.
Anybody know of cool papers I might have missed on the distribution of rain drops and the effect on their sound?
Thanks,
Andy
On Sun, 16 Mar 2008 15:43:34 -0400 marius schebella marius.schebella@gmail.com wrote:
from the first equation that andy posted, I produced a gem representation. the box muller noise seems wrong, because it does not use the whole range but is shifted to the negative side. note, this is not a distribution of frequencies, but of noise values.. marius.
Martin Peach wrote:
Oh no that's wrong isn't it :( The log is necessary to keep the distribution normal, and the range is going to get wider the closer to zero the radius is allowed to get. The attached patch has a scale adjustment... Still I wonder what kind of distribution gaussianoise2 gives, it's not just white.
Martin
Martin Peach wrote:
Charles Henry wrote:
On Sun, Mar 16, 2008 at 11:16 AM, Martin Peach martin.peach@sympatico.ca wrote:
(gaussianoise has occasional values that exceed [-1 ... 1], which I suppose is normal...white noise is always on [-1...1])
That's true. With the Box-Muller method, there is the log(~U1) term, but you can always just add a small value to U1, which will truncate your distribution. The size of the small value can be calculated to fit with any given threshold.
I think it's really because the Box-Muller method selects random numbers in pairs which map to points in a unit square on the plane, but then selects only those points which are inside the unit circle, something that the pd patch doesn't do (how to resample points in a dsp vector until they are in range?). The attached patch shows the straightforward way of doing it by simply selecting a random radius and angle and returning the resulting y coordinate as the random number. The results are always on [-1,1]. I don't think sin~ will be any slower than log~.
Martin
PD-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list
PD-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list
Here's a histogram generator (binner) that shows the distribution of [gaussianoise]. Using it I can quickly see that [gaussianoise2] is too peaked around zero and that [gaussianoise3] chops the tails off when the scale is low. If you have uniformly distributed raindrops falling, any given area will receive a number of raindrops that clusters about the mean in a normal distribution, just as if you first bin the number of occurrences of each value of white noise, then bin the resulting counts, the histogram of the counts will look like a bell curve centered at the mean count.
Martin
Andy Farnell wrote:
GEM is broken here, but thanks for the info Marius. I'm reading through the docs for R at the moment. It makes lovely plots, but haven't figured how to get my data in to it yet...
JFYI the application is rainfall. Many papers I read describe rainfall as Gaussian.
I know from physical analysis that raindrops are uniform in size and velocity for any local sample, so I've realised this distribution is about how they fall within an area and pondering how a distribution can be Gaussian in 2D.
Thing is, I can't figure out any good reason why rain should by anything other than uniformly distributed ! :(
When I use Martins second patch with a thresholding function to trigger droplet sounds, it does sound a lot more like real rainfall than a uniformly triggered model.
I'm in one of those grey areas where I half understand what I'm doing, which is a dangerous place to be.
Anybody know of cool papers I might have missed on the distribution of rain drops and the effect on their sound?
Thanks,
Andy
On Sun, 16 Mar 2008 15:43:34 -0400 marius schebella marius.schebella@gmail.com wrote:
from the first equation that andy posted, I produced a gem representation. the box muller noise seems wrong, because it does not use the whole range but is shifted to the negative side. note, this is not a distribution of frequencies, but of noise values.. marius.
Martin Peach wrote:
Oh no that's wrong isn't it :( The log is necessary to keep the distribution normal, and the range is going to get wider the closer to zero the radius is allowed to get. The attached patch has a scale adjustment... Still I wonder what kind of distribution gaussianoise2 gives, it's not just white.
Martin
Martin Peach wrote:
Charles Henry wrote:
On Sun, Mar 16, 2008 at 11:16 AM, Martin Peach martin.peach@sympatico.ca wrote:
(gaussianoise has occasional values that exceed [-1 ... 1], which I suppose is normal...white noise is always on [-1...1])
That's true. With the Box-Muller method, there is the log(~U1) term, but you can always just add a small value to U1, which will truncate your distribution. The size of the small value can be calculated to fit with any given threshold.
I think it's really because the Box-Muller method selects random numbers in pairs which map to points in a unit square on the plane, but then selects only those points which are inside the unit circle, something that the pd patch doesn't do (how to resample points in a dsp vector until they are in range?). The attached patch shows the straightforward way of doing it by simply selecting a random radius and angle and returning the resulting y coordinate as the random number. The results are always on [-1,1]. I don't think sin~ will be any slower than log~.
Martin
PD-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list
PD-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list
#N canvas 318 170 625 753 10; #X obj -243 -312 tabread bbb; #X obj -243 -376 tgl 15 0 empty empty empty 17 7 0 10 -262144 -1 -1 0 1; #X obj -243 -290 moses 0; #X obj -243 -334 f; #X obj -219 -334 + 1; #X obj -195 -334 sel 100; #X msg -146 -334 0; #X obj -204 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -257985 -13381 -1 12700 1; #X obj -189 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4034 -13381 -1 12700 1; #X obj -174 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4034 -13381 -1 12700 1; #X obj -159 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4034 -13381 -1 12700 1; #X obj -144 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4034 -13381 -1 12700 1; #X obj -129 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4034 -13381 -1 12700 1; #X obj -114 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4034 -13381 -1 12700 1; #X obj -99 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4034 -13381 -1 12700 1; #X obj -84 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4034 -13381 -1 12700 1; #X obj -243 -269 * -1; #X obj -69 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4034 -13381 -1 0 1; #X obj -204 -261 bin 0.1; #X obj -101 -340 s reset; #X obj -101 -361 bng 15 250 50 0 empty empty empty 17 7 0 10 -4034 -13381 -1; #X obj -189 -241 bin 0.2; #X obj -174 -220 bin 0.3; #X obj -159 -199 bin 0.4; #X obj -144 -179 bin 0.5; #X obj -129 -159 bin 0.6; #X obj -114 -139 bin 0.7; #X obj -99 -118 bin 0.8; #X obj -54 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4032 -13381 -1 0 1; #X obj -39 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4032 -13381 -1 0 1; #X obj -24 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4032 -13381 -1 0 1; #X obj -9 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4032 -13381 -1 0 1; #X obj 6 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4032 -13381 -1 0 1; #X obj 21 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4032 -13381 -1 0 1; #X obj 36 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4032 -13381 -1 0 1; #X obj 51 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4032 -13381 -1 0 1; #X obj 66 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4032 -13381 -1 0 1; #X obj 81 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4032 -13381 -1 0 1; #X obj -219 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -257985 -13381 -1 12700 1; #X obj -234 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4034 -13381 -1 12700 1; #X obj -249 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4034 -13381 -1 12700 1; #X obj -264 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4034 -13381 -1 12700 1; #X obj -279 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4034 -13381 -1 12700 1; #X obj -294 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4034 -13381 -1 12700 1; #X obj -309 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4034 -13381 -1 12700 1; #X obj -324 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4034 -13381 -1 12700 1; #X obj -339 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4034 -13381 -1 12700 1; #X obj -354 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4034 -13381 -1 300 1; #X obj -219 -199 bin 0.1; #X obj -234 -179 bin 0.2; #X obj -249 -159 bin 0.3; #X obj -264 -139 bin 0.4; #X obj -279 -119 bin 0.5; #X obj -294 -99 bin 0.6; #X obj -309 -79 bin 0.7; #X obj -324 -59 bin 0.8; #X obj -369 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4032 -13381 -1 0 1; #X obj -384 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4032 -13381 -1 0 1; #X obj -399 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4032 -13381 -1 0 1; #X obj -414 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4032 -13381 -1 0 1; #X obj -429 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4032 -13381 -1 0 1; #X obj -444 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4032 -13381 -1 0 1; #X obj -459 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4032 -13381 -1 0 1; #X obj -474 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4032 -13381 -1 0 1; #X obj -489 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4032 -13381 -1 0 1; #X obj -504 210 vsl 15 128 0 127 0 0 empty empty empty 0 -9 0 10 -4032 -13381 -1 0 1; #X obj -84 -95 bin 0.9; #X obj -69 -77 bin 1; #X obj -54 -55 bin 1.1; #X obj -39 -35 bin 1.2; #X obj -24 -14 bin 1.3; #X obj -9 6 bin 1.4; #X obj 6 26 bin 1.5; #X obj 21 46 bin 1.6; #X obj 36 67 bin 1.7; #X obj 51 87 bin 1.8; #X obj 66 109 bin 1.9; #X obj 81 129 bin 2; #X obj -339 -39 bin 0.9; #X obj -354 -19 bin 1; #X obj -369 1 bin 1.1; #X obj -384 21 bin 1.2; #X obj -399 41 bin 1.3; #X obj -414 61 bin 1.4; #X obj -429 81 bin 1.5; #X obj -444 101 bin 1.6; #X obj -459 121 bin 1.7; #X obj -474 141 bin 1.8; #X obj -489 163 bin 1.9; #X obj -504 183 bin 2; #X obj -243 -356 metro 1; #X obj -448 -204 tabwrite~ bbb; #X obj -438 -237 table bbb; #X obj -448 -259 bng 15 250 50 0 empty empty empty 17 7 0 10 -262144 -1 -1; #X text -432 -261 write to table; #X obj -367 -360 metro 100; #X obj -367 -335 delay 10; #X obj -367 -375 tgl 15 0 empty empty make_histogram 17 7 0 10 -4034 -13381 -13381 0 1; #X obj -461 -302 gaussianoise; #X connect 0 0 2 0; #X connect 1 0 90 0; #X connect 2 0 16 0; #X connect 2 1 18 0; #X connect 3 0 4 0; #X connect 3 0 0 0; #X connect 4 0 5 0; #X connect 5 0 6 0; #X connect 5 1 3 1; #X connect 6 0 3 1; #X connect 6 0 1 0; #X connect 16 0 48 0; #X connect 18 0 7 0; #X connect 18 1 21 0; #X connect 20 0 19 0; #X connect 21 0 8 0; #X connect 21 1 22 0; #X connect 22 0 9 0; #X connect 22 1 23 0; #X connect 23 0 10 0; #X connect 23 1 24 0; #X connect 24 0 11 0; #X connect 24 1 25 0; #X connect 25 0 12 0; #X connect 25 1 26 0; #X connect 26 0 13 0; #X connect 26 1 27 0; #X connect 27 0 14 0; #X connect 27 1 66 0; #X connect 48 0 38 0; #X connect 48 1 49 0; #X connect 49 0 39 0; #X connect 49 1 50 0; #X connect 50 0 40 0; #X connect 50 1 51 0; #X connect 51 0 41 0; #X connect 51 1 52 0; #X connect 52 0 42 0; #X connect 52 1 53 0; #X connect 53 0 43 0; #X connect 53 1 54 0; #X connect 54 0 44 0; #X connect 54 1 55 0; #X connect 55 0 45 0; #X connect 55 1 78 0; #X connect 66 0 15 0; #X connect 66 1 67 0; #X connect 67 0 17 0; #X connect 67 1 68 0; #X connect 68 0 28 0; #X connect 68 1 69 0; #X connect 69 0 29 0; #X connect 69 1 70 0; #X connect 70 0 30 0; #X connect 70 1 71 0; #X connect 71 0 31 0; #X connect 71 1 72 0; #X connect 72 0 32 0; #X connect 72 1 73 0; #X connect 73 0 33 0; #X connect 73 1 74 0; #X connect 74 0 34 0; #X connect 74 1 75 0; #X connect 75 0 35 0; #X connect 75 1 76 0; #X connect 76 0 36 0; #X connect 76 1 77 0; #X connect 77 0 37 0; #X connect 78 0 46 0; #X connect 78 1 79 0; #X connect 79 0 47 0; #X connect 79 1 80 0; #X connect 80 0 56 0; #X connect 80 1 81 0; #X connect 81 0 57 0; #X connect 81 1 82 0; #X connect 82 0 58 0; #X connect 82 1 83 0; #X connect 83 0 59 0; #X connect 83 1 84 0; #X connect 84 0 60 0; #X connect 84 1 85 0; #X connect 85 0 61 0; #X connect 85 1 86 0; #X connect 86 0 62 0; #X connect 86 1 87 0; #X connect 87 0 63 0; #X connect 87 1 88 0; #X connect 88 0 64 0; #X connect 88 1 89 0; #X connect 89 0 65 0; #X connect 90 0 3 0; #X connect 93 0 91 0; #X connect 95 0 93 0; #X connect 95 0 96 0; #X connect 96 0 1 0; #X connect 97 0 95 0; #X connect 98 0 91 0;
#N canvas 773 84 262 273 10; #X obj 103 68 inlet; #X obj 103 184 outlet; #X obj 148 184 outlet; #X obj 103 91 moses $1; #X obj 103 138 f; #X obj 125 138 + 1; #X obj 46 91 r reset; #X msg 46 121 0; #X obj 103 114 b; #X connect 0 0 3 0; #X connect 3 0 8 0; #X connect 3 1 2 0; #X connect 4 0 5 0; #X connect 4 0 1 0; #X connect 5 0 4 1; #X connect 6 0 7 0; #X connect 7 0 4 0; #X connect 8 0 4 0;
Wow, that's a gorgeous demonstraton Martin!
Everything becomes clear as time -> infinity :)
And somehow our little Earthling brains are able to spot this signature distribution as we listen to rainfall.
Now I'm getting how uniform fall leads to a Gaussian bell around the mean for an area over time.
Thanks.
Chuck, I'm sorry I couldn't follow all of your derivation of Box Muller, but thanks for the analysis. I think we agree it's a neat trick for an efficient source of WGN.
thanks all,
Andy
On Sun, 16 Mar 2008 16:54:29 -0400 Martin Peach martin.peach@sympatico.ca wrote:
Here's a histogram generator (binner) that shows the distribution of [gaussianoise]. Using it I can quickly see that [gaussianoise2] is too peaked around zero and that [gaussianoise3] chops the tails off when the scale is low. If you have uniformly distributed raindrops falling, any given area will receive a number of raindrops that clusters about the mean in a normal distribution, just as if you first bin the number of occurrences of each value of white noise, then bin the resulting counts, the histogram of the counts will look like a bell curve centered at the mean count.
Martin
Andy Farnell wrote:
GEM is broken here, but thanks for the info Marius. I'm reading through the docs for R at the moment. It makes lovely plots, but haven't figured how to get my data in to it yet...
JFYI the application is rainfall. Many papers I read describe rainfall as Gaussian.
I know from physical analysis that raindrops are uniform in size and velocity for any local sample, so I've realised this distribution is about how they fall within an area and pondering how a distribution can be Gaussian in 2D.
Thing is, I can't figure out any good reason why rain should by anything other than uniformly distributed ! :(
When I use Martins second patch with a thresholding function to trigger droplet sounds, it does sound a lot more like real rainfall than a uniformly triggered model.
I'm in one of those grey areas where I half understand what I'm doing, which is a dangerous place to be.
Anybody know of cool papers I might have missed on the distribution of rain drops and the effect on their sound?
Thanks,
Andy
On Sun, 16 Mar 2008 15:43:34 -0400 marius schebella marius.schebella@gmail.com wrote:
from the first equation that andy posted, I produced a gem representation. the box muller noise seems wrong, because it does not use the whole range but is shifted to the negative side. note, this is not a distribution of frequencies, but of noise values.. marius.
Martin Peach wrote:
Oh no that's wrong isn't it :( The log is necessary to keep the distribution normal, and the range is going to get wider the closer to zero the radius is allowed to get. The attached patch has a scale adjustment... Still I wonder what kind of distribution gaussianoise2 gives, it's not just white.
Martin
Martin Peach wrote:
Charles Henry wrote:
On Sun, Mar 16, 2008 at 11:16 AM, Martin Peach martin.peach@sympatico.ca wrote: > (gaussianoise has occasional values that exceed [-1 ... 1], which I > suppose is normal...white noise is always on [-1...1]) That's true. With the Box-Muller method, there is the log(~U1) term, but you can always just add a small value to U1, which will truncate your distribution. The size of the small value can be calculated to fit with any given threshold.
I think it's really because the Box-Muller method selects random numbers in pairs which map to points in a unit square on the plane, but then selects only those points which are inside the unit circle, something that the pd patch doesn't do (how to resample points in a dsp vector until they are in range?). The attached patch shows the straightforward way of doing it by simply selecting a random radius and angle and returning the resulting y coordinate as the random number. The results are always on [-1,1]. I don't think sin~ will be any slower than log~.
Martin
PD-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list
PD-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list
I just neatened that up into an abstration + help
All vanilla Replaced [abs~] More efficient [q8-sqrt~] seens fine No need for pi multiplier as is implicit in [cos~] radians (?) Martins histogram in separate GOP abs
If I made a mistake please correct and repost.
a.
On Sun, 16 Mar 2008 21:13:04 +0000 Andy Farnell padawan12@obiwannabe.co.uk wrote:
Wow, that's a gorgeous demonstraton Martin!
Everything becomes clear as time -> infinity :)
And somehow our little Earthling brains are able to spot this signature distribution as we listen to rainfall.
Now I'm getting how uniform fall leads to a Gaussian bell around the mean for an area over time.
Thanks.
Chuck, I'm sorry I couldn't follow all of your derivation of Box Muller, but thanks for the analysis. I think we agree it's a neat trick for an efficient source of WGN.
thanks all,
Andy
On Sun, 16 Mar 2008 16:54:29 -0400 Martin Peach martin.peach@sympatico.ca wrote:
Here's a histogram generator (binner) that shows the distribution of [gaussianoise]. Using it I can quickly see that [gaussianoise2] is too peaked around zero and that [gaussianoise3] chops the tails off when the scale is low. If you have uniformly distributed raindrops falling, any given area will receive a number of raindrops that clusters about the mean in a normal distribution, just as if you first bin the number of occurrences of each value of white noise, then bin the resulting counts, the histogram of the counts will look like a bell curve centered at the mean count.
Martin
Andy Farnell wrote:
GEM is broken here, but thanks for the info Marius. I'm reading through the docs for R at the moment. It makes lovely plots, but haven't figured how to get my data in to it yet...
JFYI the application is rainfall. Many papers I read describe rainfall as Gaussian.
I know from physical analysis that raindrops are uniform in size and velocity for any local sample, so I've realised this distribution is about how they fall within an area and pondering how a distribution can be Gaussian in 2D.
Thing is, I can't figure out any good reason why rain should by anything other than uniformly distributed ! :(
When I use Martins second patch with a thresholding function to trigger droplet sounds, it does sound a lot more like real rainfall than a uniformly triggered model.
I'm in one of those grey areas where I half understand what I'm doing, which is a dangerous place to be.
Anybody know of cool papers I might have missed on the distribution of rain drops and the effect on their sound?
Thanks,
Andy
On Sun, 16 Mar 2008 15:43:34 -0400 marius schebella marius.schebella@gmail.com wrote:
from the first equation that andy posted, I produced a gem representation. the box muller noise seems wrong, because it does not use the whole range but is shifted to the negative side. note, this is not a distribution of frequencies, but of noise values.. marius.
Martin Peach wrote:
Oh no that's wrong isn't it :( The log is necessary to keep the distribution normal, and the range is going to get wider the closer to zero the radius is allowed to get. The attached patch has a scale adjustment... Still I wonder what kind of distribution gaussianoise2 gives, it's not just white.
Martin
Martin Peach wrote:
Charles Henry wrote: > On Sun, Mar 16, 2008 at 11:16 AM, Martin Peach > martin.peach@sympatico.ca wrote: >> (gaussianoise has occasional values that exceed [-1 ... 1], which I >> suppose is normal...white noise is always on [-1...1]) > That's true. With the Box-Muller method, there is the log(~U1) term, > but you can always just add a small value to U1, which will truncate > your distribution. The size of the small value can be calculated to > fit with any given threshold. > I think it's really because the Box-Muller method selects random numbers in pairs which map to points in a unit square on the plane, but then selects only those points which are inside the unit circle, something that the pd patch doesn't do (how to resample points in a dsp vector until they are in range?). The attached patch shows the straightforward way of doing it by simply selecting a random radius and angle and returning the resulting y coordinate as the random number. The results are always on [-1,1]. I don't think sin~ will be any slower than log~.
Martin
PD-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list
PD-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list
-- Use the source
PD-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list
Okay thanks for nfo.
@ Miller, please, could we get [abs~] and [ln~] into vanilla. I think we all agree they are bread and butter objects.
cheers,
Andy
On Mon, 17 Mar 2008 16:59:08 +0100 Steffen Juul stffn@dibidut.dk wrote:
Hey Andy,
On 16/03/2008, at 23.12, Andy Farnell wrote:
I just neatened that up into an abstration + help
Thanks for wrapping it up.
All vanilla
I don't think [ln~] is vanilla. But [expr~ ln($v1)] could maybe do,
as it's shipped with vanilla.Best, Steffen
On Mon, 2008-03-17 at 16:21 +0000, Andy Farnell wrote:
@ Miller, please, could we get [abs~] and [ln~] into vanilla. I think we all agree they are bread and butter objects.
i do agree.
btw: thank you all for this very interesting thread. it was very insightful. i really liked following it, although i didn't actively participate.
roman
___________________________________________________________ Der frühe Vogel fängt den Wurm. Hier gelangen Sie zum neuen Yahoo! Mail: http://mail.yahoo.de
marius schebella wrote:
Andy Farnell wrote:
Okay thanks for nfo.
@ Miller, please, could we get [abs~] and [ln~] into vanilla. I think we all agree they are bread and butter objects.
and >~, >=~, <~, <=~
I can't get those objects to instantiate with pd-extended on Ubuntu, although some other zexy objects load OK, maybe it relates to the funny characters in the name? They also don't work on WinXP but that might be for another reason... Also I just committed a rojo~ object in svn externals/mrpeach, it uses the Box Muller method to generate red noise, and has scale and alpha inlets to adjust the quality with. The output looks good using the binner abstraction.
Martin
# Hi, sorry for jumping into the conversation. I am doing my statistical physics homework and can not read whole mails just skim read them and saw terms like "uniform distribution" and "gaussian distribution". I thought it could be worth to mention the "central limit theorem" which says that the cumulative effect of every kind of distribution will become a gaussian distribution:
http://en.wikipedia.org/wiki/Illustration_of_the_central_limit_theorem http://en.wikipedia.org/wiki/Central_limit_theorem -ugur guney-
On Sun, Mar 16, 2008 at 11:13 PM, Andy Farnell padawan12@obiwannabe.co.uk wrote:
Wow, that's a gorgeous demonstraton Martin!
Everything becomes clear as time -> infinity :)
And somehow our little Earthling brains are able to spot this signature distribution as we listen to rainfall.
Now I'm getting how uniform fall leads to a Gaussian bell around the mean for an area over time.
Thanks.
Chuck, I'm sorry I couldn't follow all of your derivation of Box Muller, but thanks for the analysis. I think we agree it's a neat trick for an efficient source of WGN.
thanks all,
Andy
On Sun, 16 Mar 2008 16:54:29 -0400
Martin Peach martin.peach@sympatico.ca wrote:
Here's a histogram generator (binner) that shows the distribution of [gaussianoise]. Using it I can quickly see that [gaussianoise2] is too peaked around zero and that [gaussianoise3] chops the tails off when the scale is low. If you have uniformly distributed raindrops falling, any given area will receive a number of raindrops that clusters about the mean in a normal distribution, just as if you first bin the number of occurrences of each value of white noise, then bin the resulting counts, the histogram of the counts will look like a bell curve centered at the mean count.
Martin
Andy Farnell wrote:
GEM is broken here, but thanks for the info Marius. I'm reading through the docs for R at the moment. It makes lovely plots, but haven't figured how to get my data in to it yet...
JFYI the application is rainfall. Many papers I read describe rainfall as Gaussian.
I know from physical analysis that raindrops are uniform in size and velocity for any local sample, so I've realised this distribution is about how they fall within an area and pondering how a distribution can be Gaussian in 2D.
Thing is, I can't figure out any good reason why rain should by anything other than uniformly distributed ! :(
When I use Martins second patch with a thresholding function to trigger droplet sounds, it does sound a lot more like real rainfall than a uniformly triggered model.
I'm in one of those grey areas where I half understand what I'm doing, which is a dangerous place to be.
Anybody know of cool papers I might have missed on the distribution of rain drops and the effect on their sound?
Thanks,
Andy
On Sun, 16 Mar 2008 15:43:34 -0400 marius schebella marius.schebella@gmail.com wrote:
from the first equation that andy posted, I produced a gem representation. the box muller noise seems wrong, because it does not use the whole range but is shifted to the negative side. note, this is not a distribution of frequencies, but of noise values.. marius.
Martin Peach wrote:
Oh no that's wrong isn't it :( The log is necessary to keep the distribution normal, and the range is going to get wider the closer to zero the radius is allowed to get. The attached patch has a scale adjustment... Still I wonder what kind of distribution gaussianoise2 gives, it's not just white.
Martin
Martin Peach wrote:
Charles Henry wrote: > On Sun, Mar 16, 2008 at 11:16 AM, Martin Peach > martin.peach@sympatico.ca wrote: >> (gaussianoise has occasional values that exceed [-1 ... 1], which I >> suppose is normal...white noise is always on [-1...1]) > That's true. With the Box-Muller method, there is the log(~U1) term, > but you can always just add a small value to U1, which will truncate > your distribution. The size of the small value can be calculated to > fit with any given threshold. > I think it's really because the Box-Muller method selects random numbers in pairs which map to points in a unit square on the plane, but then selects only those points which are inside the unit circle, something that the pd patch doesn't do (how to resample points in a dsp vector until they are in range?). The attached patch shows the straightforward way of doing it by simply selecting a random radius and angle and returning the resulting y coordinate as the random number. The results are always on [-1,1]. I don't think sin~ will be any slower than log~.
Martin
PD-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list
PD-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list
-- Use the source
PD-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list
One thing I've been thinking of lately for pd is to synthesize a impulse train with the poisson process. Many natural phenomena, like the timing of rain drops, can be simulated with the poisson process. It's just one of those ubitquitous relations that pops up everywhere.
A poisson process models a chain of events that occur with a density with respect to time. Each of the inter-event times is an exponential variable, which is easily calculatable, just as in the previous example. Take U1 a uniform variable on (0,1] , and calculate ln(U1). You've simulated a draw (observed value) from the exponential distribution, just like that. So, it's easy to make a message chain which does this calculation and outputs a series of bangs. What confuses me is whether it will be accurate in timing, and how to turn it into a series of impulses. I'm not sure I know how to do the rest of this problem.
My interest in making an impulse train is to make the impulses, which will have subtle deviations about a central frequency, and process them through a filter (formant).
I wonder what the spectrum of the poisson process looks like...
Chuck
On Sun, 16 Mar 2008 18:11:50 -0500 "Charles Henry" czhenry@gmail.com wrote:
My interest in making an impulse train is to make the impulses, which will have subtle deviations about a central frequency, and process them through a filter (formant).
Would be useful. It's certainly common, saw this recently looking at engines. A petrol engine is timed by the spark, but the sputter you get with a desiel engine seems to come from the ignition being a thermodynamic property of the compression, sometimes its a bit early or late in the cycle, but the chance of being very early or late falls off rapidly around a point.
(sorry for the duplicate, Martin!)
The Box-Muller method (I've always thought it was just Ross who did this one) is a classic trick. It probably goes back even to Gauss--who knows and who cares
pdf of Gaussian:
f(x)= k * e^-(x^2/ (2*sigma^2)) k is a normalization constant--which we will determine
The trick is, you can't integrate f(x) directly. What we do is create a joint probability distribution of g(x,y)=f(x)*f(y). x and y are independently distributed.
Then, if we integrate g(x,y) on the whole plane, we get the square of the distribution, and we can find the square of the normalization constant.
g(x,y)=f(x)= k * e^-(x^2/ (2*sigma^2)) * k * e^-(y^2/ (2*sigma^2)) =k^2 * e^-(x^2 + y^2 / (2*sigma^2))
And it's simple to change to polar coordinates: x^2+y^2=r^2 and We won't even need to care about the rest.
integral( r=0 to inf, theta=0 to 2pi, g(r,theta)*r*dr ) <-the r*dr term makes this integrable =integral( r=0 to inf, theta=0 to 2pi, k^2 * e^-(r^2 / (2*sigma^2)) * r * dr dtheta )
Integrate by theta, and make the substitution u=r^2 / 2, du= r*dr
= 2*pi * k^2 * integral( u=0 to inf, e^-(u/sigma^2) du ) = - 2*pi * sigma^2 * k^2 * ( e^-(u/sigma^2) evaluated at u=inf and u=0) = 2*pi * sigma^2 * k^2
so, for normalization: 1 = 2*pi * sigma^2 * k^2 , k=1/sqrt(2*pi*sigma^2)
So, here's the method, and the reason it works.
In g(r, theta) , r and theta are independent. Theta doesn't even appear in the distribution. And this means that theta is uniformly distributed on [0,2*pi).
So, we can find g(r) independently by integrating with respect to theta.
we get g(r)= 1/sigma^2 * e^-(r^2/sigma^2)
Now, r^2 is exponentially distributed on 0 to inf. What we're going to do is create the cumultive distribution function of r^2 which takes values on [0,1]. Then, we can choose a value on [0,1] and find the corresponding value of r^2 (and inverse function!)
our cdf (cumulative dist function) G(Z)=P( Z<r^2 )= integral ( r^2 = 0 to Z, 1/sigma^2 * e^-(r^2/sigma^2) d(r^2)
G(Z)=1 - e^-(Z/sigma^2)
Now, we can just take sigma to be one, here and proceed with the following inversion. You'll come up with the same formula regardless.
Take U1=Z on [0,1] , take U2 on [0,1]
And we're going to take those values on the unit square, map them onto r^2 and theta, and then transform back to x and y, giving us two independently distributed normal variables.
U1=1 - e^-(Z/sigma^2)
e^-(Z/sigma^2) = 1 - U1
Z/sigma^2 = - ln (1-U1)
Z= r^2= -1*sigma^2 * ln (1-U1)
r = sqrt(-1*sigma^2 * ln (1-U1) ) r = sigma * sqrt(-1 * ln (1-U1) )
theta=2*pi*U2
x=r*cos(theta) = sigma * sqrt(-1 * ln (1-U1) ) * cos (2*pi*U2) y=r*sin(theta) = sigma * sqrt(-1 * ln (1-U1) ) * sin (2*pi*U2)
Notice, the sigma comes out in front, and the variable (1-U1) is distributed uniformly on [0,1] also. Hence, it can be simplified to another uniform variable U3, or whatever.
This is probably the *most* classic trick in statistics. Enjoy!
Chuck
On Sun, Mar 16, 2008 at 2:11 PM, Martin Peach martin.peach@sympatico.ca wrote:
Oh no that's wrong isn't it :( The log is necessary to keep the distribution normal, and the range is going to get wider the closer to zero the radius is allowed to get. The attached patch has a scale adjustment... Still I wonder what kind of distribution gaussianoise2 gives, it's not just white.
Martin
So, we can find g(r) independently by integrating with respect to theta.
we get g(r)= 1/sigma^2 * e^-(r^2/sigma^2)
Ah, here's a missing factor of 2 in the exponent. That should be
g(r)= 1/sigma^2 * e^-(r^2/ ( 2*sigma^2) )
And the correct formulae are
x=r*cos(theta) = sigma * sqrt(-2 * ln (1-U1) ) * cos (2*pi*U2) y=r*sin(theta) = sigma * sqrt(-2 * ln (1-U1) ) * sin (2*pi*U2)
I thought I missed something...
On Sun, Mar 16, 2008 at 10:30 AM, Andy Farnell padawan12@obiwannabe.co.uk wrote:
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?
It should sound the same and look the same spectrally. There was a previous discussion on the list about how noise distributions (uniform vs gaussian pdfs) should sound. Namely, the same :)
To see the difference, you could try to calculate kurtosis values! I've never had a good use for that one. http://en.wikipedia.org/wiki/Kurtosis That statistic will show you the difference between Gaussian and uniformly distributed noise, and confirm your distribution.
The Box Muller method should be good--that's a very reliable method.
Cheers, Chuck
Cheers all,
Andy
-- Use the source
PD-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list
-
Andy Farnell -
Charles Henry -
Lawrence Joseph -
marius schebella -
Martin Peach -
Roman Haefeli -
Steffen Juul -
Uğur Güney