Generating 1/f (flicker/pink) noise

Hi Allan,

On Wed, 26 Mar 2008 00:49:49 +1100, Allan Herriman wrote:

> You might like to point out on your webpage that this method of
> generating pink noise has has O(nlog(n)) time complexity and O(n)
> storage requirements for generating n samples.  This limits it to
> relatively small sample sizes (which, nevertheless, may be large enough
> for a particular need).
> 
> Other methods (e.g. additive or subtractive synthesis, mentioned
> elsewhere in this thread) have O(1) time storage requirements, and
> are suitable for streaming continuous data.

[Allan means O(1) time per sample, of course.  Anyone doing streaming 
signal processing thinks per-sample as a matter of course.  Although in 
the context of a batch algorithm, that probably needs a bit of 
clarification.]

I've not read Max's algorithm, but I'm pretty sure that I can rely on the 
independence of the pink stream from any output process to use arbitrary-
precision FIR filtering of white noise to make streaming pink noise as 
nice as you could want, with O(n) time and O(1) storage requirements, 
batch-fashion.  The constants of proportionality might be a bit larger 
than yours...  On the other hand, if I care to use overlap-add FFT 
convolution to run the FIR filter, then I suspect that I can make my 
amortized computation-per-sample cost smaller than any of the proposed 
pink-specific filtering algorithms.  The cost is enough memory to get the 
per-sample compute cost and spectrum accuracy right, and a multi-threaded 
execution environment that lets me keep computing new blocks of noise in 
the background.  I accept that there are plenty of environments where 
neither of those costs can be borne.

Cheers,

-- 
Andrew

On 25 Mrz., 19:04, Allan Herriman <allanherri...@hotmail.com> wrote:
> On Wed, 26 Mar 2008 04:29:21 +1100, Allan Herriman
>
>
>
>
>
> <allanherri...@hotmail.com> wrote:
> >On Tue, 25 Mar 2008 08:03:25 -0700 (PDT), Andor
> ><andor.bari...@gmail.com> wrote:
>
> >>On 25 Mrz., 15:43, Allan Herriman <allanherri...@hotmail.com> wrote:
> >>> On Tue, 25 Mar 2008 07:23:21 -0700 (PDT), Andor
>
> >>> <andor.bari...@gmail.com> wrote:
> >>> >On 25 Mrz., 14:49, Allan Herriman <allanherri...@hotmail.com> wrote:
> >>> >> On Fri, 21 Mar 2008 10:38:39 -0700 (PDT), maxlittle2...@googlemail.com
> >>> >> wrote:
>
> >>> >> >Dear All
>
> >>> >> >After having spent a frustrating hour searching around to find some
> >>> >> >simple code that would generate 1/f noise, or more generally 1/f^alpha
> >>> >> >noise, I decided to write it myself. So here goes:
>
> >>> >> >http://www.eng.ox.ac.uk/samp/powernoise_soft.html
>
> >>> >> >I hope you enjoy it. Comments and feedback welcome.
>
> >>> >> >All the best,
> >>> >> >Max
>
> >>> >> You might like to point out on your webpage that this method of
> >>> >> generating pink noise has has O(nlog(n)) time complexity and O(n)
> >>> >> storage requirements for generating n samples. &#4294967295;This limits it to
> >>> >> relatively small sample sizes (which, nevertheless, may be large
> >>> >> enough for a particular need).
>
> >>> >> Other methods (e.g. additive or subtractive synthesis, mentioned
> >>> >> elsewhere in this thread) have O(1) time and storage requirements, and
> >>> >> are suitable for streaming continuous data.
>
> >>> >I find that hard to believe. The storage and the processing time is
> >>> >independent of the number of samples to be computed?
>
> >>> Short term brain malfunction, sorry. &#4294967295;
>
> >>> I meant to say: O(n) time, and O(1) storage (if streaming) or O(n)
> >>> storage (if saving the data somewhere).
>
> >>Ok, that makes more sense. Inspired by this thread (and other
> >>circumstances), I have read the Kasdin reference that I posted above
> >>somewhere. There it is shown that the LTI filtering techniques to
> >>construct 1/f noise either
>
> >>1. have a finite number of coefficients (either transversal or
> >>recursive) and consequently, the 1/f-law breaks down for frequencies
> >>below some cutoff frequency
>
> >>or
>
> >>2. use filters with growing number of filter coefficients (either
> >>transversal or recursive) as the simulation proceeds, resulting in
> >>sequences with asymptotically unbiased PSDs (true 1/f spectra).
>
> >>This seems natural with my understanding of the autocorrelation
> >>function of 1/f-noises, which only decays hyperbolically. Using
> >>filtered sequences with a constant number of coefficients, the
> >>autocorrelation function is either truncated (FIR filters) or decays
> >>exponentially (stable rational transfer function filter) or does not
> >>decay at all (unstable rational transfer function filters, for example
> >>to construct 1/f^2 noise).
>
> >>Are any of the methods that only require O(n) computation time
> >>asymptotically unbiased?
>
> >Yes, the Voss pink noise generator (additive synthesis) requires O(n)
> >computation time (specifically: two calls to rand() and a small amount
> >of overhead per output sample), regardless of the output bandwidth.
>
> Minor clarification: for each output sample, there are two calls to
> rand(); a small amount of overhead and the summation of a list of
> saved random numbers.
> The length of that list is proportional to the number of octaves
> generated. &#4294967295;So I guess my claim of "regardless of the output
> bandwidth" is false if a naive implementation is used.

Some interesting algorithms are presented here:

http://planck.lal.in2p3.fr/IMG/pdf/astro-ph-full.pdf

The "random midpoint"-algorithm requires only one randn() call per
sample! However, it recursively generates a 1/f^alpha sequence by
starting with the start and the endpoint, so it's not usable for
streaming.

The main algorithm is to sum 6dB/oct lowpass noises with
logarithmically spaced cutoff frequencies (a rather old idea). The
fractional differencing method (Kasdin, Hosking) is also mentioned.

>
> Since only one element in the list changes each sample, one could use
> the old "boxcar filter recursive FIR" trick (as used in CIC filters)
> to give a constant amount of processing per output sample.
>
>
>
> >The storage requirement is, however, directly proportional to the
> >number of octaves generated. &#4294967295;
> >As an engineer (rather than a mathematician) I could claim that a
> >modest amount of storage would get you close enough to any desired
> >accuracy at any frequency close to (but not including) DC.

The one thing I remember from my statistics class was that the
defining characteristic was the asymptotic behaviour of the PSD
towards DC. Anything with bounded PSD misses the nature of long-range
correlated 1/f noises.

Regards,
Andor

On Tue, 25 Mar 2008 15:19:43 -0700, Andor wrote:

> The one thing I remember from my statistics class was that the defining
> characteristic was the asymptotic behaviour of the PSD towards DC.
> Anything with bounded PSD misses the nature of long-range correlated 1/f
> noises.

While that's true, and a useful property for lots of estimation and 
theoretical applications, one of the really big uses of "pink" noise is 
as a standard test signal in audio systems.  (And that's because it 
almost matches the PSD of lots of more interesting signals.)  The systems 
are generally known to be band-limited to something like human hearing, 
so it's OK if the signal is too.  (In fact, since they're also usually 
bounded in an absolute sense, it's not too useful to be generating values 
asymptotically approaching infinity.)

Cheers,

-- 
Andrew

On 25 Mar, 14:51, dbd <d...@ieee.org> wrote:
> On Mar 24, 4:51 pm, maxlittle2...@googlemail.com wrote:
>
>
>
> > On 23 Mar, 18:33, dbd <d...@ieee.org> wrote:
>
> > > On Mar 22, 4:03 pm, maxlittle2...@googlemail.com wrote:
>
> > > > ...
> > > > Ah, well I did find those resources, but no simple "one pager" of
> > > > Matlab code for generating perfect 1/f^alpha noise.
> > > > ...
> > > > Max
>
> > > There can be no page that is 'simple' and 'perfect' for all
> > > applications.
>
> > > One characteristic of your implementation choice is that the power
> > > spectrum does not just have the desired expectation value in any
> > > frequency region, but it equals the expectation value. While this is a
> > > possible power spectum value, it is not a likely one, and is unlikely
> > > to repeat. This can be useful in certain testing situations. It might
> > > be inappropriate in others. An alternative is to use independent
> > > Gaussian samples for real and imaginary frequency components and
> > > weight them with the desired frequency weighting. An example
> > > implementation of this, modified from your implementation, follows
> > > below. This could provide more realistic statistical behavior in some
> > > applications.
>
> > > Dale B. Dalrymplehttp://dbdimages.com
>
> > Dale
>
> > Thanks, these are excellent points. I've added an option to the
> > routine so that a choice can be made between deterministic and
> > stochastic power.  I've also added an option that normalises the
> > signal amplitude so that it lies in the range -1 to +1. This may be
> > useful in certain circumstances. Grab the code from:
>
> >http://www.eng.ox.ac.uk/samp/powernoise_soft.html
>
> > Max
>
> Max
>
> If you really think a normalization and a peak limiting capability
> should be included inside the function, I suggest you do so with a
> gain and a clipping capability instead of putting the mean and scale
> of the output under the control of the extremal behavior of the noise
> generated.
>
> There are test purposes where the original function or my addition
> could be called multiple times to generate longer sequences of test
> data. Even these special cases no longer work when gain and offset
> change from block to block.
>
> Dale B. Dalrymple

Thanks for the thoughtful suggestions Dale.

I think the normalization to the range [-1, 1] is required for some
applications (particularly nonlinear time series analysis): at least
one member of this group has found that some gain is desireable, and I
reason that it's better to have a signal that uses a standard number
range. Otherwise, as you imply, there probably isn't much point in
putting this scaling in the function. If by clipping I understand you
to mean some kind of nonlinear "clamping" of values outside a certain
range, I don't agree with this because it will distort the spectrum.

This also leads on to your second paragraph. Unfortunately this
method, as it stands, can't be called multiple times to generate
longer sequences: there will be several problems. One of them is that
the signal is not guaranteed to be continuous across multiple calls.
The resulting discontinuities will introduce spectral distortions.

Of course, one could imagine a solution to fix this problem that
matches the endpoints of each successive call, but then you'd also
need to match as many derivatives of the signal as are required to
ensure the appropriate differentiability for the chosen "roughness" of
the signal (which, informally, is inversely proportional to alpha). I
don't think such an approach would be elegant though, and I'm sure
there's a better one, although I can't think of it right now.

Max

On 25 Mar, 22:19, Andor <andor.bari...@gmail.com> wrote:

> Some interesting algorithms are presented here:
>
> http://planck.lal.in2p3.fr/IMG/pdf/astro-ph-full.pdf
>
> The "random midpoint"-algorithm requires only one randn() call per
> sample! However, it recursively generates a 1/f^alpha sequence by
> starting with the start and the endpoint, so it's not usable for
> streaming.

Nice reference, Andor, and a very nice contribution to the discussion
here.

I particularly like the logistic map algorithm suggested in the paper.

Max

On Mar 25, 5:58 pm, maxlittle2...@googlemail.com wrote:
> On 25 Mar, 14:51, dbd <d...@ieee.org> wrote:
>
>
>
> > On Mar 24, 4:51 pm, maxlittle2...@googlemail.com wrote:
>
> > > On 23 Mar, 18:33, dbd <d...@ieee.org> wrote:
>
> > > > On Mar 22, 4:03 pm, maxlittle2...@googlemail.com wrote:
>
> > > > > ...
> > > > > Ah, well I did find those resources, but no simple "one pager" of
> > > > > Matlab code for generating perfect 1/f^alpha noise.
> > > > > ...
> > > > > Max
>
> > > > There can be no page that is 'simple' and 'perfect' for all
> > > > applications.
>
> > > > One characteristic of your implementation choice is that the power
> > > > spectrum does not just have the desired expectation value in any
> > > > frequency region, but it equals the expectation value. While this is a
> > > > possible power spectum value, it is not a likely one, and is unlikely
> > > > to repeat. This can be useful in certain testing situations. It might
> > > > be inappropriate in others. An alternative is to use independent
> > > > Gaussian samples for real and imaginary frequency components and
> > > > weight them with the desired frequency weighting. An example
> > > > implementation of this, modified from your implementation, follows
> > > > below. This could provide more realistic statistical behavior in some
> > > > applications.
>
> > > > Dale B. Dalrymplehttp://dbdimages.com
>
> > > Dale
>
> > > Thanks, these are excellent points. I've added an option to the
> > > routine so that a choice can be made between deterministic and
> > > stochastic power.  I've also added an option that normalises the
> > > signal amplitude so that it lies in the range -1 to +1. This may be
> > > useful in certain circumstances. Grab the code from:
>
> > >http://www.eng.ox.ac.uk/samp/powernoise_soft.html
>
> > > Max
>
> > Max
>
> > If you really think a normalization and a peak limiting capability
> > should be included inside the function, I suggest you do so with a
> > gain and a clipping capability instead of putting the mean and scale
> > of the output under the control of the extremal behavior of the noise
> > generated.
>
> > There are test purposes where the original function or my addition
> > could be called multiple times to generate longer sequences of test
> > data. Even these special cases no longer work when gain and offset
> > change from block to block.
>
> > Dale B. Dalrymple
>
> Thanks for the thoughtful suggestions Dale.
>
> I think the normalization to the range [-1, 1] is required for some
> applications (particularly nonlinear time series analysis): at least
> one member of this group has found that some gain is desireable, and I
> reason that it's better to have a signal that uses a standard number
> range. Otherwise, as you imply, there probably isn't much point in
> putting this scaling in the function. If by clipping I understand you
> to mean some kind of nonlinear "clamping" of values outside a certain
> range, I don't agree with this because it will distort the spectrum.
>
> This also leads on to your second paragraph. Unfortunately this
> method, as it stands, can't be called multiple times to generate
> longer sequences: there will be several problems. One of them is that
> the signal is not guaranteed to be continuous across multiple calls.
> The resulting discontinuities will introduce spectral distortions.
>
> Of course, one could imagine a solution to fix this problem that
> matches the endpoints of each successive call, but then you'd also
> need to match as many derivatives of the signal as are required to
> ensure the appropriate differentiability for the chosen "roughness" of
> the signal (which, informally, is inversely proportional to alpha). I
> don't think such an approach would be elegant though, and I'm sure
> there's a better one, although I can't think of it right now.
>
> Max

If you want to rerun a test and generate a new data set to do so, or
run a test against multiple independent data sets, the RMS and DC
values of the next data set will be entirely different from the last
one as you have it. A knee-jerk reaction to clipping seems easier to
some than thinking about the actual effect. If for example, you limit
the amplitude to put the data through a DAC, you should consider how
the error generated by clipping compares with the level of
quantization noise plus the mismatch between the generated data and
'real' 1/f noise. For almost anything except examining the tail ends
of order statistic distributions the clipping of some extremal values
is undetectable in noise. Differences in scaling and offset are much
more obvious to many common processes.

Dale B. Dalrymple

On 26 Mar, 02:02, dbd <d...@ieee.org> wrote:
> On Mar 25, 5:58 pm, maxlittle2...@googlemail.com wrote:
>
>
>
> > On 25 Mar, 14:51, dbd <d...@ieee.org> wrote:
>
> > > On Mar 24, 4:51 pm, maxlittle2...@googlemail.com wrote:
>
> > > > On 23 Mar, 18:33, dbd <d...@ieee.org> wrote:
>
> > > > > On Mar 22, 4:03 pm, maxlittle2...@googlemail.com wrote:
>
> > > > > > ...
> > > > > > Ah, well I did find those resources, but no simple "one pager" of
> > > > > > Matlab code for generating perfect 1/f^alpha noise.
> > > > > > ...
> > > > > > Max
>
> > > > > There can be no page that is 'simple' and 'perfect' for all
> > > > > applications.
>
> > > > > One characteristic of your implementation choice is that the power
> > > > > spectrum does not just have the desired expectation value in any
> > > > > frequency region, but it equals the expectation value. While this is a
> > > > > possible power spectum value, it is not a likely one, and is unlikely
> > > > > to repeat. This can be useful in certain testing situations. It might
> > > > > be inappropriate in others. An alternative is to use independent
> > > > > Gaussian samples for real and imaginary frequency components and
> > > > > weight them with the desired frequency weighting. An example
> > > > > implementation of this, modified from your implementation, follows
> > > > > below. This could provide more realistic statistical behavior in some
> > > > > applications.
>
> > > > > Dale B. Dalrymplehttp://dbdimages.com
>
> > > > Dale
>
> > > > Thanks, these are excellent points. I've added an option to the
> > > > routine so that a choice can be made between deterministic and
> > > > stochastic power.  I've also added an option that normalises the
> > > > signal amplitude so that it lies in the range -1 to +1. This may be
> > > > useful in certain circumstances. Grab the code from:
>
> > > >http://www.eng.ox.ac.uk/samp/powernoise_soft.html
>
> > > > Max
>
> > > Max
>
> > > If you really think a normalization and a peak limiting capability
> > > should be included inside the function, I suggest you do so with a
> > > gain and a clipping capability instead of putting the mean and scale
> > > of the output under the control of the extremal behavior of the noise
> > > generated.
>
> > > There are test purposes where the original function or my addition
> > > could be called multiple times to generate longer sequences of test
> > > data. Even these special cases no longer work when gain and offset
> > > change from block to block.
>
> > > Dale B. Dalrymple
>
> > Thanks for the thoughtful suggestions Dale.
>
> > I think the normalization to the range [-1, 1] is required for some
> > applications (particularly nonlinear time series analysis): at least
> > one member of this group has found that some gain is desireable, and I
> > reason that it's better to have a signal that uses a standard number
> > range. Otherwise, as you imply, there probably isn't much point in
> > putting this scaling in the function. If by clipping I understand you
> > to mean some kind of nonlinear "clamping" of values outside a certain
> > range, I don't agree with this because it will distort the spectrum.
>
> > This also leads on to your second paragraph. Unfortunately this
> > method, as it stands, can't be called multiple times to generate
> > longer sequences: there will be several problems. One of them is that
> > the signal is not guaranteed to be continuous across multiple calls.
> > The resulting discontinuities will introduce spectral distortions.
>
> > Of course, one could imagine a solution to fix this problem that
> > matches the endpoints of each successive call, but then you'd also
> > need to match as many derivatives of the signal as are required to
> > ensure the appropriate differentiability for the chosen "roughness" of
> > the signal (which, informally, is inversely proportional to alpha). I
> > don't think such an approach would be elegant though, and I'm sure
> > there's a better one, although I can't think of it right now.
>
> > Max
>
> If you want to rerun a test and generate a new data set to do so, or
> run a test against multiple independent data sets, the RMS and DC
> values of the next data set will be entirely different from the last
> one as you have it. A knee-jerk reaction to clipping seems easier to
> some than thinking about the actual effect. If for example, you limit
> the amplitude to put the data through a DAC, you should consider how
> the error generated by clipping compares with the level of
> quantization noise plus the mismatch between the generated data and
> 'real' 1/f noise. For almost anything except examining the tail ends
> of order statistic distributions the clipping of some extremal values
> is undetectable in noise. Differences in scaling and offset are much
> more obvious to many common processes.
>
> Dale B. Dalrymple

Thanks Dale, very interesting points as always.

You're right, of course, in that things are nearly always more complex
than they initially appear. I'd certainly be interested to see exact
analytical results detailing the distortion of the 1/f^alpha power
spectrum caused by a given level of clipping.

I do think though that, at least superficially, if you asked around
for a range of opinions (and I mean not just in signals circles), by
"normalization" most would respond that the amplitude is somehow
scaled to some numerical range that is widely understood to be a
"unit" range of some sort. I don't know how many would think of
clipping initially. But that's just my opinion.

Max

On Mar 25, 7:44 pm, maxlittle2...@googlemail.com wrote:

> ...
>
> Thanks Dale, very interesting points as always.
>
> You're right, of course, in that things are nearly always more complex
> than they initially appear. I'd certainly be interested to see exact
> analytical results detailing the distortion of the 1/f^alpha power
> spectrum caused by a given level of clipping.

To first order, each peak clip is equivalent to summing a short
impulse (the clipped region, inverted). This becomes broadband noise.
Is it's level small compared to the main noise signal? That's the
test.

>
> I do think though that, at least superficially, if you asked around
> for a range of opinions (and I mean not just in signals circles), by
> "normalization" most would respond that the amplitude is somehow
> scaled to some numerical range that is widely understood to be a
> "unit" range of some sort. I don't know how many would think of
> clipping initially. But that's just my opinion.
>
> Max

Traditional users of noise-like signals have thought in terms of
amplitude being RMS and not instantaneous. They wouldn't notice
clipping until it showed on the oscilloscope. (There's no clipping
noise from my load resister when I turn up the gain, so no problem!)
New players in the digital world haven't learned how to analyze what
happens or if it matters when they replace calculated values with the
max or min values their digital domain enforces. I'm sure you are
right that neither would think of clipping. But why would you trust
either to tell you how a digital generator should behave?

All noise generators with active gain stages in the noise generation
path have clipping. Theoretically, the noise sources must be capable
of arbitrarily large outputs. The power supplies present limits to the
levels that can be 'represented' even in analog instruments, so
clipping is not new. It's just still not well recognized.

Dale B. Dalrymple

On 25 Mrz., 23:41, Andrew Reilly <andrew-newsp...@areilly.bpc-
users.org> wrote:
> On Tue, 25 Mar 2008 15:19:43 -0700, Andor wrote:
> > The one thing I remember from my statistics class was that the defining
> > characteristic was the asymptotic behaviour of the PSD towards DC.
> > Anything with bounded PSD misses the nature of long-range correlated 1/f
> > noises.
>
> While that's true, and a useful property for lots of estimation and
> theoretical applications, one of the really big uses of "pink" noise is
> as a standard test signal in audio systems. &#4294967295;(And that's because it
> almost matches the PSD of lots of more interesting signals.)

Agreed. For more occurences of 1/f noise see this list:

http://www.nslij-genetics.org/wli/1fnoise/index-by-category.html

> The systems
> are generally known to be band-limited to something like human hearing,
> so it's OK if the signal is too. &#4294967295;(In fact, since they're also usually
> bounded in an absolute sense, it's not too useful to be generating values
> asymptotically approaching infinity.)

If I remember correctly, the random walk (1/f^2) is 0-recurrent,
meaning that it always returns back to its starting point. The
expected length between two 0-crossings is, however, infinite.

Regards,
Andor

On 26 Mrz., 02:08, maxlittle2...@googlemail.com wrote:
> On 25 Mar, 22:19, Andor <andor.bari...@gmail.com> wrote:
>
> > Some interesting algorithms are presented here:
>
> >http://planck.lal.in2p3.fr/IMG/pdf/astro-ph-full.pdf
>
> > The "random midpoint"-algorithm requires only one randn() call per
> > sample! However, it recursively generates a 1/f^alpha sequence by
> > starting with the start and the endpoint, so it's not usable for
> > streaming.
>
> Nice reference, Andor, and a very nice contribution to the discussion
> here.
>
> I particularly like the logistic map algorithm suggested in the paper.

I thought that part in the paper was a bit odd. The author asks how to
simulate 1/f^alpha noise for alpha > 2 (and suggests the logisitc
map). However, the question can be trivially answered. If alpha > 2,
first simulate noise with PSD 1/f^beta, where

beta = alpha mod 2,

then filter the 1/f^beta noise with an integrator of order floor(alpha/
2). Voila.

Regards,
Andor