*From*: Eamon Caddigan <ecaddiga@xxxxxxxx>*To*: visionegg@xxxxxxxxxxxxx*Date*: Mon, 8 Dec 2008 13:24:00 -0600

Hi all,

-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()

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