The companion page discusses some presumably new ideas about using correlated rather than white noise generator stages for producing digitally-synthesized approximate pink noise, with a randomization strategy that eliminates the consistent spectral inaccuracies that the Voss-McCartney (V-M) algorithm produces. The V-M algorithm still has some advantages: it is efficient, its has very uniform processor loading for predictable timing in real-time applications, it is extensible, and the spectral anomalies, though large, are localized and probably not audible. For a description of that algorithm see the definitive site for pink noise generation, http://www.firstpr.com.au/dsp/pink-noise/.

With a randomized updating policy rather than a fixed one, it is not possible to know in advance which generator stages will update at any given time, so the computational loading at any given instant varies from nothing to the worst case that every generator stage is updated. The goal here is to apply the correlated generator approach but reduce the uneven loading and its drawbacks. If possible: improve accuracy and increase the spectrum coverage.

I will have to report only modest success. There is some
residual approximation error, and very little remaining freedom
to try to improve upon it. I consider the design as shown in
this page the * end of the line * — but of course,
there could be some new theoretical twists that bring some new life
to the idea.

Like the Voss-McCartney algorithm, this approach maintains a set of generator stages, updated at different rates. The difference is that the updating rate is only a statistical average, and is never performed on a predictable schedule (which is the source of the spectrum anomalies that plague the V-M results).

The new twist is this: the updating of the generator stages is
made conditional, rather than independent. Stage 1 is updated
*only* on the condition that Stage 2 through Stage N are
not updated! Similar idea for the rest of the stages. This
produces a different random process than when the stages are
updated completely independently, but the generators retain the
same probabilities of updating. By updating at most one stage at
any given update opportunity, almost uniform processor loading is
achieved.

There is a price for this. Whether one generator can update at any given time depends on what the other generators are doing, so there is no longer complete statistical independence. This lack of independence shows up as signals are combined. There is just enough difference to make the helpful theory about signal energy distributions unusable to guide the design. So I gave my best shot for normalizing the residual errors manually. The results shown here are not optimal, but close enough that you gain much by reaching a true optimum.

Okay, here is the algorithm that includes a coefficient set to
cover approximately 9 octaves conforming to the "pink" 1/f power
distribution. First, we need to define the parameters. The
`pP`

and the `pA`

parameters are the
*amplitude scaling* and *probability of update*
parameters described in the previous note.

pA = [ 3.8024 2.9694 2.5970 3.0870 3.4006 ] pP = [ 0.00198 0.01280 0.04900 0.17000 0.68200 ] pSUM = [ 0.00198 0.01478 0.06378 0.23378 0.91578 ]

The `pSUM` values are obtained by forming a running
sum of the terms in the `pP` array. You should note that
the values must sum to less than 1.0 for this algorithm to work.
You can think of this as a *cumulative probability distribution.*

You also need two other things. You need a pseudo-random
number generator `rand`

that produces reasonably random values
in the range 0.0 to 1.0, suitable for signal generation. You need an
accumulator array `contrib`

that maintains the values currently
contributed by each generator stage. You can initialize it randomly
or fill with zeroes.

contrib = zeros(Nterms)

Now for each new time sample, perform the following update to determine the pink noise generator value.

ur1 = rand; for stage=1:Nstages if (ur1 <= pSUM(stage)) ur2 = rand; contrib(istage)=2*(ur2-0.5)*pA(istage); break for; endif endfor genout = sum(contrib);

The update always uses one random value to determine which stages (if any) to update. There is a small probability that none of the stages are updated. An update operation uses a second random value. Auto-correlation results when values are not updated, with statistical properties discussed in the previous note. The number of random values required is 1 or 2, averaging a little less than 2, so the cost of evaluation is almost the same as the V-M algorithm.

If you are using exact (fixed point) arithmetic, you can avoid summing the values of the generator stages to obtain the final output. Maintain the combined sum in another accumulator variable. To update, after a generator stage is selected, subtract the previous value of its contribution from the output accumulator, save the new contribution level, and add the new value into the accumulator. The generator output then equals the value of the accumulator variable.

The scaling shown is very arbitrary. The worst possible case for a fixed point implementation is when all of the random values reach 1.0 at the same time. With the scaling parameters shown, the generator would output its maximum level of 15.8564. Observing that 32767/15.8564 = 2066.48, you can scale all of the amplitude coefficients by 2066.48 to obtain a generator that does not exceed a signed 16-bit fixed point range.

You can balance the loading to a uniform 2 random evaluations per update by generating and discarding a value in the special case where none of the stages are updated.

If you have a 32-bit random number generator, you can apply some additional arithmetic to avoid one of the random number evaluations. A 32-bit generator might cost more than a 16-bit generator twice as shown. Otherwise, observe:

- If stage 1 is updated, we know the value of
`ur1`

is in the interval 0.0000 to 0.00198. Values of a uniform distribution are uniformly distributed in a subset of the interval, so we can derive`ur2 = ur1/pP(1)`

- If stage 2 is updated, we know the value of
`ur1`

is in the interval 0.00198 to 0.1280. Values of a uniform distribution are uniformly distributed in a subset of the interval, so we can derive`ur2 = (ur1-pSUM(1))/pP(2)`

- Similar for stage 3,
`urv = (uru-pSUM(2))/pP(3)`

- ... and so forth

You can get slightly faster execution by reversing the order of
the stages. I present them from lowest frequency to highest, but the
number of branches taken during evaluation is minimized using the
reverse order. If you do this, be sure that the `pSUM` array
is recalculated to match (don't just reverse the numbers)!

The previous note provided 3-stage, 4-stage and 5-stage generators. There is really no point in the 3-stage or 4-stage generators any more. If it turns out that there is too much low-frequency power, you can high-pass filter the noise sequence.

Andreas Robinson has observed that the pinkness property of the result depends on the updating schedule, not on the statistical distribution from which the original zero-mean uncorrelated random numbers are taken. The uniform distribution is relatively arbitrary choice. Instead of uniform random generators, the generators could be Gaussian and everything still works. Or going to the other extreme, the distribution can be a balanced two-valued discrete distribution. This particular choice requires only one bit of precision for its representation, and some interesting things happen. The updating formula can be greatly simplified for this case as follows.

ur2 = 0 with probability 0.5; 1 with probability 0.5 contrib(istage)= // general form = 2*(ur2-0.5)*pA(istage); contrib(istage) // simplified form = +pA(istage) if ur2 == 1 = -pA(istage) if ur2 == 0

The multiply operation is no longer needed for the updating. Furthermore, as there are only two possible levels contributed to the composite generator output by each stage, the contributions of all stages combine to only 32 possible distinct output levels. These can be precomputed to avoid the rest of the real-time arithmetic, making the algorithm very friendly indeed for implementations using a processor with limited memory, small register width, and limited arithmetic capability. Look here for more details about his approach (including a hardware implementation that avoids the precomputed table as well!)

Not all pink noise signals are created equal, and I suspect that this is not they way I would build a pink signal to drive a voltage-controlled lowpass filter, for example. But for many testing and direct audio applications, the signal spectrum is all that really matters and this might work very well.

Efficiency is close to as good as you can get of course. While there is no place in the spectrum where the signal energy is off to the degree that it is with the V-M algorithm, the anomalies there are highly localized, whereas here the anomalies in this algorithm are smaller and spread.

A plot of approximation errors is shown below. Using the parameter values given above for the generator, the x-axis range is in natural logarithms and corresponds to roughly the range 30 Hz to 18 kHz at 44kHz sampling, a little over 9 octaves. You can see the response droop near zero frequency, as expected because it is not possible to fit the singularity at frequency zero. You can also see where the levels starts to deviate from the 1/f model at the higher frequency end, approaching the Nyquist frequency. There is little power left at the very high frequencies so this shouldn't make any difference. The signal power should equal the 1/f model, hence the ratio should be 1.0. Computing the decibels for the ratio between actual response and the ratio curve, the plot should show a level line at 0.0 dB. As you can see, the deviations are about 0.25 dB, which corresponds to maximum amplitude deviation of about 2.9 percent. This should be close enough for all but the most demanding laboratory applications.

Just to give a little getter perspective on this, the new algorithm (5 generator terms) was run back-to-back against the Voss algorithm (12 generator terms). The Voss curve is artificially lowered so the two data sets don't overlap. Otherwise conditions are identical. The fuzz is uncertainty in the estimates, due to the limited sample size (just a few million points!) but the systematic variations are clearly visible. Draw your own conclusions.

Whether the new algorithm is the right one for you probably depends on how you intend to use it. It is easier to implement than the patterned scheduling of the V-M algorithm, with audibility good enough that you probably won't be able to distinguish it from a "true" pink process. If you need a laboratory grade noise source so you can depend on the correct energy distribution everywhere in the spectrum very accurately over the long term, this is not the best way to go.

Site: RidgeRat's Barely Operating Site http://home.earthlink.net/~ltrammell Author: Larry Trammell Created: Jan 19, 2006 Revised: Jan 07, 2008 Status: Experimental Contact: NOSPAM ltramme1476 AT earthlink DOT net NOSPAM Related: See the 'definitive pink noise site' at http://www.firstpr.com.au/dsp/pink-noise/ Restrictions: This information is public domain. Thanks to Andreas Robinson for contributing his variant to this site.