# Fast 2d brown/pink noise generation?

Discussion in 'Python' started by Mathias, Mar 4, 2005.

1. ### MathiasGuest

Dear NG,

I'm looking for a fast way to produce 2d-noise images with 1/f or 1/f^2
spectrum. I currently generate the noise via inverse FFT, but since I
need lots of data (~10^7 for a monte carlo simulation) it needs to be
really fast. Does someone know a faster way than my approach?

- Dimensionality is between 20x20 and 100x100
- The spectrum doesn't need to be exactly pink/brown, an approximation
is fine.
- Implementation in either matlab or scientific python (LAPACK anyway)

Thanks a lot for hints, literature, algorithms!
Mathias

Mathias, Mar 4, 2005

2. ### Robert KernGuest

Mathias wrote:
> Dear NG,
>
> I'm looking for a fast way to produce 2d-noise images with 1/f or 1/f^2
> spectrum. I currently generate the noise via inverse FFT, but since I
> need lots of data (~10^7 for a monte carlo simulation) it needs to be
> really fast. Does someone know a faster way than my approach?
>
> - Dimensionality is between 20x20 and 100x100
> - The spectrum doesn't need to be exactly pink/brown, an approximation
> is fine.
> - Implementation in either matlab or scientific python (LAPACK anyway)

This is a 1D version that I have using scipy. It's naive, so I'm sure
that it is slower. However, I believe the general technique can be
implemented on a larger scale.

The basic idea is to sum up a bunch of white time series with different
time steps. The first level is white noise at every time step. The
second level changes at every second time step. The third changes at
every fourth, etc.

I think you can replicate this by generating a few white noise arrays of
the appropriate sizes, judiciously using repeat(), and summing them
together. I got this scheme from an article I found by googling for pink
noise algorithms, I believe.

from scipy import *

class PinkGenerator(object):
updateTable = [0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4]
updateTable.extend(updateTable[:-1])
del updateTable[-1]

def __init__(self, rng=stats.norm):
self.key = 0
self.rng = rng
self.whiteValues = self.rng.rvs(size=5)

def getNextValue(self):
self.key += 1
self.key = self.key % len(self.updateTable)
self.whiteValues[self.updateTable[self.key]] = self.rng.rvs()[0]
return (sum(self.whiteValues) + self.rng.rvs()[0])/6

def getManyValues(self, size):
data = zeros((size,), Float)
for i in range(size):
data = self.getNextValue()
return data

def sampleData(self, size=1024):
data = self.getManyValues(size)
p = power(absolute(fftshift(fft(data))), 2)/size
f = fftshift(fftfreq(size))
return data, f, p

--
Robert Kern

"In the fields of hell where the grass grows high
Are the graves of dreams allowed to die."
-- Richard Harter

Robert Kern, Mar 5, 2005

## Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.