# [visionegg] Pink Noise in VisionEgg

• To: visionegg@xxxxxxxxxxxxx
• Date: Mon, 8 Dec 2008 13:24:00 -0600

```Hi all,

```
I took my first crack at writing a pink noise stimulus for vision egg. Most of the work was done by Jon Yearsley, who wrote the original matlab code upon which this is based; I just transcribed his routine into python and wrapped it into a subclass of VisionEgg.Textures.Texture.
```
```
I'm posting here in case anybody else finds it useful. Please let me know if you spot any bugs; I'm particularly suspicious of the apparent need to rescale the data.
```

```
```#!/usr/bin/env python
"""Code for 2D 1/f noise. This module provides the spatialPattern function,
used to generate such noise, and the PinkNoiseTexture class, a subclass of
VisionEgg.Textures.Texture, for use in vision experiments.

Eamon Caddigan, University of Illinois, 2008-12-08"""

import numpy
import VisionEgg.Textures

def spatialPattern(DIM, BETA=-1.0):
"""This function generates 1/f spatial noise, with a normal error
distribution (the grid must be at least 10x10 for the errors to be normal).
1/f noise is scale invariant, there is no spatial scale for which the
variance plateaus out, so the process is non-stationary.

DIM is a two component vector that sets the size of the spatial pattern
(DIM=[10,5] is a 10x5 spatial grid)
BETA defines the spectral distribution.
Spectral density S(f) = N f^BETA
(f is the frequency, N is normalisation coeff).
BETA = 0 is random white noise.
BETA  -1 is pink noise
BETA = -2 is Brownian noise
The fractal dimension is related to BETA by, D = (6+BETA)/2

Note that the spatial pattern is periodic.  If this is not wanted the
grid size should be doubled and only the first quadrant used.

Time series can be generated by setting one component of DIM to 1

The method is briefly descirbed in Lennon, J.L. "Red-shifts and red
herrings in geographical ecology", Ecography, Vol. 23, p101-113 (2000)

Many natural systems look very similar to 1/f processes, so generating
1/f noise is a useful null model for natural systems.

The errors are normally distributed because of the central
limit theorem.  The phases of each frequency component are randomly
assigned with a uniform distribution from 0 to 2*pi. By summing up the
frequency components the error distribution approaches a normal
distribution.

Written by Jon Yearsley  1 May 2004
j.yearsley@xxxxxxxxxxxxxx

S_f corrected to be S_f = (u.^2 + v.^2).^(BETA/2);  2/10/05

Transcribed to Python/numpy by Eamon Caddigan, 2008"""

# Type checking
if len(DIM) == 1:
DIM = (DIM, DIM)
elif len(DIM) != 2:
raise ValueError("Size must be a one- or two-element sequence")

# Generate the grid of frequencies. u is the set of frequencies along the
# first dimension
# First quadrant are positive frequencies.  Zero frequency is at u[0, 0].
u = numpy.append(numpy.arange(0, numpy.floor(0.5*DIM[0])+1),
numpy.arange(-numpy.ceil(0.5*DIM[0])+1, 0)) / DIM[0]
# Reproduce these frequencies along every row
u = numpy.tile(u, [DIM[1], 1]).transpose()
# v is the set of frequencies along the second dimension.  For a square
# region it will be the transpose of u
v = numpy.append(numpy.arange(0, numpy.floor(0.5*DIM[1])+1),
numpy.arange(-numpy.ceil(0.5*DIM[1])+1, 0)) / DIM[1]
# Reproduce these frequencies along every column
v = numpy.tile(v, [DIM[0], 1])

# Generate the power spectrum
S_f = numpy.power(u**2 + v**2, float(BETA)/2)

# Set any infinities to zero
S_f[S_f==numpy.inf] = 0

# Generate a grid of random phase shifts
phi = numpy.random.random(DIM)

# Inverse Fourier transform to obtain the the spatial pattern
x = numpy.fft.ifft2(numpy.sqrt(S_f) *
(numpy.cos(2*numpy.pi*phi) + 1j*numpy.sin(2*numpy.pi*phi)))

# Pick just the real component
x = numpy.real(x)

return(x)

class PinkNoiseTexture(VisionEgg.Textures.Texture):
"""1/f^alpha Noise texture for VisionEgg"""
def __init__(self, size=(256,256), alpha=1):
"""Create a 1/f^alpha noise texture."""

# Generate noise, scale to [0, 1]
texels = spatialPattern(size, -alpha)
texels = texels - numpy.min(texels)
texels = 1.0/numpy.max(texels) * texels
texels = texels.transpose()

# Create Texture object
VisionEgg.Textures.Texture.__init__(self, texels, size)

if __name__ == "__main__":
import VisionEgg
VisionEgg.start_default_logging(); VisionEgg.watch_exceptions()
import VisionEgg.FlowControl

screen = VisionEgg.Core.get_default_screen()

stim = VisionEgg.Textures.TextureStimulus(texture=PinkNoiseTexture(),
position=(screen.size[0]/2, screen.size[1]/2), anchor="center",
mipmaps_enabled=False)

viewport = VisionEgg.Core.Viewport(screen=screen, stimuli=[stim])

p = VisionEgg.FlowControl.Presentation(viewports=[viewport],
go_duration=(5.0, 'seconds'))
p.go()

```

### Other related posts:

• » [visionegg] Pink Noise in VisionEgg - Eamon Caddigan