Pink ($1/f$) pseudo-random noise generation

Question

What are some algorithms for generating a good pseudo-random approximation to $1/f$ (pink) noise, yet suitable for implementation with low computational cost on an integer DSP?

Answer 1

There are several. This site has a reasonable (but possibly old) list:

• Filtering white noise to make it pink.
• The Voss algorithm.
• The Voss-McCartney algorithm.

Answer 2

Linear Filtering

The first approach in Peter's answer (i.e. filtering white noise) is a very straightforward approach. In Spectral Audio Signal Processing, JOS gives a low-order filter that can be used to produce a decent approximation, along with an analysis of how well the resulting power spectral density matches the ideal. Linear filtering will always yield an approximation, but that may not matter in practice. To paraphrase JOS:

There is no exact (rational, finite-order) filter which can produce pink noise from white noise. This is because the ideal amplitude response of the filter must be proportional to the irrational function $1/\sqrt{f}$, where $f$ denotes frequency in Hz. However, it is easy enough to generate pink noise to any desired degree of approximation, including perceptually exact.

The coefficients of the filter he gives are as follows:

B = [0.049922035, -0.095993537, 0.050612699, -0.004408786];
A = [1, -2.494956002, 2.017265875, -0.522189400];

They're formatted as parameters to the MATLAB filter function, so for the sake of clarity, they correspond to the following transfer function:

$$ H(z) = { .041 - .096z^{-1} + .051z^{-2} - .004z^{-3} \over{} 1 -2.495z^{-1} + 2.017z^{-2} - .522z^{-3}} $$

Obviously, it's better to use the full precision of the coefficients in practice. Here's a link to what the pink noise generated using that filter sounds like:

For fixed-point implementation, since it's usually more convenient to work with coefficients in the range [-1,1), some reworking of the transfer function will be in order. Generally, the recommendation is to break things up into second-order sections, but part of the reason for that (as opposed to using first-order sections) is for the convenience of working with real coefficients when the roots are complex. For this particular filter, all the roots are real, and combining then into second-order sections would probably still yield some denominator coefficients >1, so three first-order sections is a reasonable choice, as follows:

$$ H(z) = { 1 - b_1z^{-1}\over{1 - a_1z^{-1}} } \space { 1-b_2z^{-1} \over{1-a_2z^{-1} } } \space {1 - b_3z^{-1} \over{ 1 - a_3z^{-1} } } $$

where

$$ b_1 = 0.98223157, \space b_2 = 0.83265661, \space b_3 = 0.10798089 $$ $$ a_1 = 0.99516897, \space a_2 = 0.94384177, \space a_3 = 0.55594526 $$

Some judicious choice of sequencing for those sections, combined with some choice of gain factors for each section will be required to prevent overflow. I haven't tried any of the other filters given in the link in Peter's answer, but similar considerations would probably apply.

White Noise

Obviously, the filtering approach requires a source of uniform random numbers in the first place. If a library routine isn't available for a given platform, one of the simplest approaches is to use a linear congruential generator. One example of an efficient fixed-point implementation is given by TI in Random Number Generation on a TMS320C5x (pdf). A detailed theoretical discussion of various other methods can be found in Random Number Generation and Monte Carlo Methods by James Gentle.

Resources

Several sources based on following links in Peter's answer are worth highlighting.

• The first filter-based chunk of code references Introduction to Signal Processing by Orfanidis. The full text is available at that link, and [in Appendix B] it has coverage of both pink and white noise generation. As the comment mentions, Orfanidis mostly covers the Voss algorithm.

• The Spectrum Produced by the Voss-McCartney Pink Noise Generator. Way down near the bottom of the page, after extensive discussion of variants of the Voss algorithm, this link is referenced in giant pink letters. It's much easier reading than some of the preceding ASCII diagrams.

• A Bibliograph on 1/f Noise by Wentian Li. This is referenced both in Peter's source and by JOS. It has a dizzying number of references on 1/f noise in general, dating all the way back to 1918.

Answer 3

I have been using Corsini and Saletti's algorithm since 1990: G. Corsini, R. Saletti, "A 1/f^gamma Power Spectrum Noise Sequence Generator", IEEE Transactions on Instrumentation and Measurement, 37(4), December, 1988, 615-619. The gamma exponent is between -2 and +2. It works well for my purposes. Ed

If this attempt to add a screenshot works, the figure below shows an example of how well the Corsini and Saletti algorithm performs (at least as I programmed it back in 1990). The sampling frequency was 1 kHz, gamma = 1, and 1000 32k FFT PSDs were averaged.

---

This follows up to my previous post on the Corsini and Saletti (C&S) noise generator. The next two figures show how well the C&S generator performs in regard to generation of low frequency (gamma > 0) and high frequency (gamma < 0) noises. The third figure compares the 1/f noise PSDs of the C&S generator (same as my first post) and the Example B.9 1/f generator given in Prof. Orfanidis’s excellent book (eqn B.29, p. 736). All of these PSDs are averages of 1000 32k FFT PSDs. They are all unilateral and mean-subtracted. For the C&S PSDs, I used 3 poles/decade and specified 4 decades (0.05 to 500 Hz) as the desired usable range. So the C&S generator had n = 12 pole and zero pairs. The sampling frequency was 1 kHz, Nyquist was 500 Hz, and the resolution element was just over 0.0305 Hz. Ed V

As Corsini & Saletti state in their paper, $f_c ≥ 10f_M$, where $f_c$ is the sampling frequency and $f_M$ is the “upper limit of the frequency band on which we will generate noise samples”. The digital filter coefficients are given by their equations (5.1): $$a_i = exp[-2\pi10^{(i-N)/h - \gamma/2h - c}]$$ $$b_i = exp[-2\pi10^{(i-N)/h - c}]$$ where c = 1. To obtain C&S PSDs like those shown above, let c = 0 and $f_M = 0.5f_c$.