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

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?

• How about memory? If that's not a concern but computation is, I'd say do a random-phase iDFT of the desired frequency curve and save it as a static const wavetable in your device. – leftaroundabout Sep 27 '11 at 12:47
• @leftaroundabout - Or would multiplying the DFT of random white noise by the 1/f frequency curve, then doing an IDFT have better randomness characteristics? – hotpaw2 Sep 28 '11 at 6:24
• White noise is essentially a random-phase iFT of the constant function, so it shouldn't make much difference. – leftaroundabout Sep 28 '11 at 9:28

• If that site goes down your answer disappears, if you were to bring over the basics of each solution the answer would be much improved, using the site as reference. – Kortuk Oct 9 '11 at 17:25
• @Kortuk: The answer is community wiki, so feel free to do so yourself! The information there should be enough to point to other web references (like datageist's answer for the first option). I agree, though, that more details would be good. – Peter K. Oct 9 '11 at 18:23

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.

• Any idea how he came up with these filter coefficients? I guess it's just a nonlinear fit to the desired slope, but I'd be quite interested to know if there is a more specific algorithm. – nibot Apr 26 '12 at 20:07
• My best guess would be one of the approximation techniques mentioned in his thesis. It's a great read either way. – datageist Apr 26 '12 at 20:23
• Wow, that is quite a document! Thanks for the link. – nibot Apr 26 '12 at 22:17
• The problem with the filter white noise method is you don't get the same magnitude phase relationships as you do with an autocorrelated time series. Thus if you are trying to emulate natural processes you should not generate white noise and filter it. You should actually create autocorrelated noise as a time series, i.e. current value depends on previous value + noise. See in statistics "AR" processes. You can test this by generating noise using both methods, then FFT, and plot real versus imaginary (complex plane of the frequency domain). You'll notice a big difference in the pattern – Paul S Feb 17 '16 at 1:10
• Hi Paul, welcome to DSP.SE. If you just care about how the noise sounds (in audio work, for example) then the magnitude spectrum is the primary concern. It would be great if you could detail your thoughts in a new answer, though. I don't think we have anything on the site describing that technique yet. – datageist Feb 17 '16 at 2:52

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$$.

• Corsini and Saletti state "This filter consists of N cascaded first-order sections, each with a real pole-zero pair" and the N poles are "uniformly distributed with respect to the frequency logarithm with a density of h poles per frequency decade (p/d), and the N zeros follow accordingly." The Discussion section of the paper was exceptionally nicely done, so there was no problem just programming what they said to do. All I have is my old hardcopy and a scanned copy of it. For the PSD above, I used 3 poles/decade and the PSD is mean-subtracted and unilateral. Ed V – Ed V Apr 23 at 22:11