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)+1), numpy.arange(-numpy.ceil(0.5*DIM)+1, 0)) / DIM # Reproduce these frequencies along every row u = numpy.tile(u, [DIM, 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), numpy.arange(-numpy.ceil(0.5*DIM)+1, 0)) / DIM # Reproduce these frequencies along every column v = numpy.tile(v, [DIM, 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/2, screen.size/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()