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(o bject):

updateTable = [0,1,0,2,0,1,0,3 ,0,1,0,2,0,1,0, 4]

updateTable.ext end(updateTable[:-1])

del updateTable[-1]

def __init__(self, rng=stats.norm) :

self.key = 0

self.rng = rng

self.whiteValue s = self.rng.rvs(si ze=5)

def getNextValue(se lf):

self.key += 1

self.key = self.key % len(self.update Table)

self.whiteValue s[self.updateTabl e[self.key]] = self.rng.rvs()[0]

return (sum(self.white Values) + self.rng.rvs()[0])/6

def getManyValues(s elf, size):

data = zeros((size,), Float)

for i in range(size):

data[i] = self.getNextVal ue()

return data

def sampleData(self , size=1024):

data = self.getManyVal ues(size)

p = power(absolute( fftshift(fft(da ta))), 2)/size

f = fftshift(fftfre q(size))

return data, f, p

--

Robert Kern

rk***@ucsd.edu
"In the fields of hell where the grass grows high

Are the graves of dreams allowed to die."

-- Richard Harter