Code for generating 1/f^alpha pink noise

Started by Sampo Niskanen October 13, 2008
Hi,

In a recent project I required a source of 1/f^(5/3) pink noise, but was 
unable to find any code or explanation on how to perform that.  The DSP 
Generation of Pink Noise page[1] contains extensive information on 
generating 1/f pink noise, but has no mention of 1/f^alpha pink noise.

I therefore made my own code for the purpose.  It produces pink noise by 
filtering white noise with an IIR filter described by N. Jeremy Kasdin, 
Proceedings of the IEEE, Vol. 83, No. 5, May 1995, p. 822.  (The 
reference was found in an old thread on this group.)

Because there seemed to be no freely available code to produce 1/f^alpha 
pink noise on the net, I made the Java class available.  The procedure 
is very simple to re-implement in other languages as well.  The code and 
a more thorough analysis of the procedure is available at 
http://www.iki.fi/sampo.niskanen/PinkNoise/


Feedback is welcome.  :)


[1] http://www.firstpr.com.au/dsp/pink-noise/

-- 
  __________________________________________________
 /____\   Sampo Niskanen <=> sampo.niskanen@iki.fi  \
       \       http://www.iki.fi/sampo.niskanen/     \
        \     ________________________________________\___
         \___/___________________________________________/
Sampo Niskanen wrote:
> Hi, > > In a recent project I required a source of 1/f^(5/3) pink noise, but was > unable to find any code or explanation on how to perform that. &#2013266080;The
DSP
> Generation of Pink Noise page[1] contains extensive information on > generating 1/f pink noise, but has no mention of 1/f^alpha pink noise. > > I therefore made my own code for the purpose. &#2013266080;It produces pink noise
by
> filtering white noise with an IIR filter described by N. Jeremy Kasdin, > Proceedings of the IEEE, Vol. 83, No. 5, May 1995, p. 822. &#2013266080;(The > reference was found in an old thread on this group.) > > Because there seemed to be no freely available code to produce 1/f^alpha > pink noise on the net, I made the Java class available. &#2013266080;The
procedure
> is very simple to re-implement in other languages as well. &#2013266080;The code
and
> a more thorough analysis of the procedure is available
athttp://www.iki.fi/sampo.niskanen/PinkNoise/
> > Feedback is welcome. &#2013266080;:)
Hi Sampo I didn't look at your implementation, but it's great that you are making this publicly available. A slight technical nitpick: by truncating the fractional integrator IIR filter it becomes an FIR filter, so your method is actually FIR based. This distinction is important, because there are quite a few IIR approximations to pink noise (notably the ones where a parallel structure of first order lowpass filters is used), which are fundamentally different than your approach. By the way: is there an error in your formula under Methodology? You write that the FIR filter coefficients are generated by the recursion a_0 = 1 a_k = (k-1-alpha/2) a_{k-1} / k However, Kasdin writes (equation (104)) a_0 = 1 a_k = (alpha/2+k-1) a_{k-1} / k Regards, Andor
I wrote:
> Sampo Niskanen wrote: > > Hi, > > > In a recent project I required a source of 1/f^(5/3) pink noise, but was > > unable to find any code or explanation on how to perform that. &#2013266080;The
DSP
> > Generation of Pink Noise page[1] contains extensive information on > > generating 1/f pink noise, but has no mention of 1/f^alpha pink noise. > > > I therefore made my own code for the purpose. &#2013266080;It produces pink
noise by
> > filtering white noise with an IIR filter described by N. Jeremy Kasdin, > > Proceedings of the IEEE, Vol. 83, No. 5, May 1995, p. 822. &#2013266080;(The > > reference was found in an old thread on this group.) > > > Because there seemed to be no freely available code to produce 1/f^alpha > > pink noise on the net, I made the Java class available. &#2013266080;The
procedure
> > is very simple to re-implement in other languages as well. &#2013266080;The code
and
> > a more thorough analysis of the procedure is available
athttp://www.iki.fi/sampo.niskanen/PinkNoise/
> > > Feedback is welcome. &#2013266080;:) > > Hi Sampo > > I didn't look at your implementation, but it's great that you are > making this publicly available. A slight technical nitpick: by > truncating the fractional integrator IIR filter it becomes an FIR > filter, so your method is actually FIR based. This distinction is > important, because there are quite a few IIR approximations to pink > noise (notably the ones where a parallel structure of first order > lowpass filters is used), which are fundamentally different than your > approach. > > By the way: is there an error in your formula under Methodology? You > write that the FIR filter coefficients are generated by the recursion > > a_0 = 1 > a_k = (k-1-alpha/2) a_{k-1} / k > > However, Kasdin writes (equation (104)) > > a_0 = 1 > a_k = (alpha/2+k-1) a_{k-1} / k
Oops. You _are_ using an IIR filter and you are using equation (116) to compute the purely recursive filter coefficients. So what you wrote above all makes sense. Did you compare the AR with the MA approach? Regards, Andor
Hi Sampo,
Thanks for the presenting the material, very interesting and useful.
Steve
1) When alpha = 1, how do your results compare to the Voss-McCartney
algorithm?

2) Is this the formula for the rolloff?
     rolloff = (alpha * -10dB/decade)?

3) How do you convert from db/decade to db/octave?
   i.e. -10dB/decade is 3.0102999 dB/octave.


Thanks for posting your code.
Stacy wrote:
> 1) When alpha = 1, how do your results compare to the Voss-McCartney > algorithm? > > 2) Is this the formula for the rolloff? > rolloff = (alpha * -10dB/decade)?
That's how I interpret it.
> 3) How do you convert from db/decade to db/octave? > i.e. -10dB/decade is 3.0102999 dB/octave.
Multiply by .30102999 (.3 is close enough). Jerry -- Engineering is the art of making what you want from things you can get. &#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095; ** Posted from http://www.teranews.com **
Thanks for all the input.


> 1) When alpha = 1, how do your results compare to the Voss-McCartney > algorithm?
I haven't compared this method to the Voss-McCartney method. In my application I required 1/alpha^(5/3) noise, which cannot be generated by the Voss-McCartney algorithm (or at least I couldn't figure out how to modify it in such a way). I haven't compared the AR and MA methods either. I'm no DSP expert and once I had a working filter that produced good enough results I went back to my original project. :) In fact, I decided to use only a two-pole filter in my application since I didn't want the values deviating from zero for long periods of time. Anyway, I hope the code is useful to others as well. -- __________________________________________________ /____\ Sampo Niskanen <=> sampo.niskanen@iki.fi \ \ http://www.iki.fi/sampo.niskanen/ \ \ ________________________________________\___ \___/___________________________________________/
Sampo Niskanen wrote:

> I haven't compared the AR and MA methods either. &#2013266080;I'm no DSP expert
and
> once I had a working filter that produced good enough results I went > back to my original project. &#2013266080;:) &#2013266080;In fact, I decided to
use only a
> two-pole filter in my application since I didn't want the values > deviating from zero for long periods of time.
Hi Sampo This sounds a bit contradicting. Either you are intersted in 1/f^{5/3} noise, or you are worried about values in the series deviating from zero for long periods. You can't have both. This deviaton for longer periods from zero is one of the very characteristics that makes 1/f^alpha noise "statistically difficult" to handle. This phenomenon was called the "Joseph Effect" by Mandelbrot, because of Joseph's dream that Egypt was to experience seven years of plentiful harvests followed by seven years of famine (long periods of deviation from the zero mean). Particularly climatological and hydrological time series can be modeled very well by 1/f^alpha noise. The Nile river data has an estimated Hurst exponent of H = 0.91, which translates to an alpha of 0.82. Note that your alpha is > 1, indicating non-stationary long range correlation. For that case long deviations from zero is the rule, not the exception. Just a thought.
> Anyway, I hope the code is useful to others as well.
I'm sure it is :-). Regards, Andor
> Note that your alpha is > 1, indicating non-stationary long > range correlation. For that case long deviations from zero > is the rule, not the exception.
You can consider alpha = 2 (Wiener process). In that case, the expected time between zero crossings is infinite. 5/3 isn't much smaller than 6/3.
Andor <andor.bariska@gmail.com> wrote:
> This sounds a bit contradicting. Either you are intersted in 1/f^{5/3} > noise, or you are worried about values in the series deviating from > zero for long periods.
> You can't have both.
The application in my case was simulating wind turbulence. Wind turbulence freqency has a power spectrum density proportional to 1/(1+K*f)^(5/3), which for large f is equal to 1/f^(5/3). The simulation time scale, however, is quite short, so I don't want a gust of wind to last the whole simulation duration (that would effectively change the average wind speed). By choosing the number of poles suitably one can choose how much low frequency components to include. The original (empirical) formula doesn't go to infinity at f=0 either, so one could estimate it with a suitable number of poles. In my case I wanted even less low-frequency components, and using 2 poles with a 20Hz sampling rate yields maximum wind gust lengths of approximately 3-5 seconds (the spectrum turns flat below 0.3Hz). So in my application, it is desireable that the noise is pink at the high end of the spectrum, but flat at the low end. -- __________________________________________________ /____\ Sampo Niskanen <=> sampo.niskanen@iki.fi \ \ http://www.iki.fi/sampo.niskanen/ \ \ ________________________________________\___ \___/___________________________________________/