[visionegg] Pink Noise in VisionEgg

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. 

-Eamon Caddigan


#!/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: