Code for generating 1/f^alpha pink noise

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. �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 athttp://www.iki.fi/sampo.niskanen/PinkNoise/
>
> Feedback is welcome. �:)

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. �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 athttp://www.iki.fi/sampo.niskanen/PinkNoise/
>
> > Feedback is welcome. �:)
>
> 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.
�����������������������������������������������������������������������
** 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. �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.

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/     \
        \     ________________________________________\___
         \___/___________________________________________/