Ruido blanco

Uno de los elementos fundamentales para la generación de contenido procedural es el ruido. Empezaremos por la forma de ruido más simple, el ruido blanco.

white_noise

Una imagen bidimensional compuesta por ruido blanco

como oirse (bajar el volumen antes de escuchar el sonido!).

Se denomina de esta forma a una señal que posee las distintas componentes espectrales “en igual proporción”. Dese el punto de vista de las señales continuas es solo una idealización ya que, si tuviera una potencia no nula en una banda espectral dada, tendría que tener una potencia total infinita. Pero en las señales discretas, al estar su frecuencia acotada superiormente, es posible tener un “espectro plano”. (También es posible definir ruido blanco en una banda de frecuencias dadas para señales continuas, pero eso no es el tema que trataremos en este post.)

Una forma común de generar ruido blanco es mediante el uso de muestras aleatorias independientes repartidas uniformemente en un intervalo. Por ejemplo, la imagen y el archivo WAV anteriormente mostrados fueron generadas mediante los siguientes scripts:

# requires PIL 1.1.6
# http://www.pythonware.com/products/pil/
import Image, random

def make_white_noise(fname, size):
    im = Image.new('L', size)
    pix = im.load()
    for i in range(size[0]):
        for j in range(size[1]):
            pix[i, j] = random.randint(0, 255)
    im.save(fname)

if __name__ == '__main__':
    make_white_noise('white_noise.png', (256, 256))

import random, struct, wave

def make_white_noise(fname, time):
wfp = wave.open(fname, ‘w’)
wfp.setparams((1, 2, 44100, 0, ‘NONE’, None))
for i in xrange(44100*time):
wfp.writeframesraw(
struct.pack(‘h’,
random.randint(-1<<15, 1<<15-1))) wfp.close() if __name__ == '__main__': make_white_noise('white_noise.wav', 1) [/code] En lo que resta de este post analizaremos porqué esto funciona…

Varianza y correlación

La varianza de una variable aleatoria es una medida de la dispersión de sus valores y se define, para el caso de una variable aleatoria real X, como:

\sigma^2(X) = E[X^2] - E[X]^2

donde E[.] representa al valor esperado. Utilizando la propiedad de linealidad del valor esperado podemos llegar a la siguiente identidad:

\sigma^2(X + Y) = E[(X+Y)^2] - E[X+Y]^2

\sigma^2(X + Y) = E[(X+Y)^2] - (E[X]+E[Y])^2

\sigma^2(X + Y) = E[X^2+2XY+Y^2] - (E[X]^2+2E[X]E[Y]+E[Y]^2)

\sigma^2(X + Y) = (E[X^2]-E[X]^2)+2(E[XY]-E[X]E[Y])+(E[Y^2]-E[Y]^2)

\sigma^2(X + Y) = \sigma^2(X)+\sigma^2(Y)+2(E[XY]-E[X]E[Y]).

El término que aparece sumado a las varianzas se denomina (bueno, en realidad la mitad de ese término se denomina…) covarianza y se anula si las variables aleatorias son independientes ya que:

E[XY] = \int_{\mathbb{R}^2} xy h(x, y) dx dy

E[XY] = \int_{\mathbb{R}^2} xy f(x) g(y) dx dy

E[XY] = \int_{\mathbb{R}} x f(x) dx \int_{\mathbb{R}} y g(y) dy

E[XY] = \int_{\mathbb{R}} x f(x) dx \int_{\mathbb{R}} y g(y) dy

E[XY] = E[X]E[Y],

donde llamamos f(x) a la densidad de probabilidad de X, g(y) a la densidad de probabilidad de Y y h(x, y) a la densidad de probabilidad conjunta, además de suponer que X e Y son independientes.

Por lo tanto, en caso de ser X e Y variables aleatorias independientes tendremos:

\sigma^2(X + Y) = \sigma^2(X)+\sigma^2(Y).

Espectro de una señal compuesta por muestras independientes y distribuidas uniformemente

Llamemos x_n a la señal compuesta por N muestras independientes distribuidas uniformemente en [-1, 1] y X_k a su DFT. Entonces tendremos

X_k = \sum_{n=0}^{N-1} x_n e^{-\frac{2 \pi i}{N} k n}.

Como E[x_n] es nulo, cada término de la suma tendrá un valor esperado nulo y, por lo tanto, E[X_k] es nulo. Pero estamos interesados en la distribución de potencia por intervalo espectral, por lo que debemos observar a E[|X_k|^2] en lugar de a E[X_k].

Aplicando directamente la definición y operando algebraicamente tenemos:

E[|X_k|^2] = E\left[\left|\sum_{n=0}^{N-1} x_n e^{-\frac{2 \pi i}{N} k n}\right|^2\right]

E[|X_k|^2] = E\left[\left|\sum_{n=0}^{N-1} x_n \left(\cos\left(-\frac{2 \pi}{N} k n\right) + i \sin\left(-\frac{2 \pi}{N} k n\right) \right) \right|^2\right]

E[|X_k|^2] = E\left[\left|\sum_{n=0}^{N-1} x_n \cos\left(-\frac{2 \pi}{N} k n\right) + i \sum_{n=0}^{N-1} x_n \sin\left(-\frac{2 \pi}{N} k n\right) \right|^2\right]

E[|X_k|^2] = E\left[\left(\sum_{n=0}^{N-1} x_n \cos\left(-\frac{2 \pi}{N} k n\right)\right)^2 + \left(\sum_{n=0}^{N-1} x_n \sin\left(-\frac{2 \pi}{N} k n\right)\right)^2 \right].

Ahora, como los valores de x_n son independientes, tenemos que la covarianza entre dos elementos distintos cualesquiera para cada una estas sumas será nula. Eso nos permite, junto con el hecho de que E[x_n] = 0, transformar la expresión anterior a la siguiente:

E[|X_k|^2] = E\left[\sum_{n=0}^{N-1} x_n^2 \cos^2\left(-\frac{2 \pi}{N} k n\right) + \sum_{n=0}^{N-1} x_n^2 \sin^2\left(-\frac{2 \pi}{N} k n\right) \right].

Aplicando una identidad bien conocida 🙂 y distribuyendo a E[.] aprovechando su linealidad, llegamos a:

E[|X_k|^2] = \sum_{n=0}^{N-1} E[x_n^2].

Aunque podríamos encontrar el valor exacto, no es de mucha importancia ya que podemos observar algo mucho más importante: no depende de k. Eso nos indica que el espectro “es plano”, el resultado que buscábamos. (En realidad lo que es plano no es “el” espectro sino su valor esperado; pero eso es lo que suele interesar en el caso de una señal aleatoria.)

FFT: Qué es realmente?

Cursando en la facultad, escuché varias veces hablar sobre la FFT e incluso la utilicé en una materia. Si bien tuvimos una preparación razonable en cuanto a las series y transformadas de Fourier, en ningún momento nos explicaron nada acerca de este algoritmo.

Esto me llevó naturalmente a suponer que debía ser un algoritmo complejo y, en consecuencia, a sorprenderme al ver que la idea básica tras él es muy simple y bien conocida por los que estudiamos Informática: “divide & conquer”.

Para empezar consideremos primero la definición de la transformada discreta de Fourier:

X_k = \sum_{n=0}^{N-1} x_n e^{-\frac{2 \pi i}{N} k n} \quad \quad k = 0, \dots, N-1

Si lo implementáramos directamente obtendríamos un algoritmo como el siguiente:

def naive_dft(X):
    n = len(X)
    basic_root = math.e ** (-2.0 * math.pi * 1j / n)
    K = [sum(x * root ** i for i, x in enumerate(X))
         for root in (basic_root ** j for j in range(n))]
    return K

Como tenemos esencialmente una iteración completa sobre los N elementos de x_n para calcular cada uno de los N elementos X_k, la complejidad del algoritmo será \Theta(N^2) (puede recortarse el número de operaciones aprovechando ciertas identidades pero el número de operaciones seguirá siendo cuadrático).

Pero la expresión dada anteriormente para la DFT puede reescribirse del siguiente modo (asumiendo que N es par y tomando M = N/2):

X_k = \sum_{n=0}^{N-1} x_n e^{-\frac{2 \pi i}{N} k n}

X_k = \sum_{m=0}^{M-1} (x_{2m} e^{-\frac{2 \pi i}{2M} k (2m)} + x_{2m+1} e^{-\frac{2 \pi i}{2M} k (2m + 1)})

X_k = \sum_{m=0}^{M-1} x_{2m} e^{-\frac{2 \pi i}{M} k m} + \sum_{m=0}^{M-1} x_{2m+1} e^{-\frac{2 \pi i}{M} k m} e^{-\frac{2 \pi i}{2M} k}

X_k = \sum_{m=0}^{M-1} x_{2m} e^{-\frac{2 \pi i}{M} k m} + e^{-\frac{2 \pi i}{N} k} \left( \sum_{m=0}^{M-1} x_{2m+1} e^{-\frac{2 \pi i}{M} k m} \right)

Puede observarse que obtuvimos dos expresiones muy similares a la expresión original para la DFT, pero tratando con secuencias con una longitud de un medio de la original. La única diferencia es que en la secuencia original podíamos garantizar que k < N, mientras que k puede ser mayor que M. Pero como

e^{-\frac{2 \pi i}{M} k m} = e^{-\frac{2 \pi i}{M} (k - a M) m} e^{-\frac{2 \pi i}{M} a M m}

e^{-\frac{2 \pi i}{M} k m} = e^{-\frac{2 \pi i}{M} (k - a M) m} e^{-2 \pi i a m}

e^{-\frac{2 \pi i}{M} k m} = e^{-\frac{2 \pi i}{M} (k - a M) m},

los valores de las DFT parciales simplemente se repetirán para K = M, \dots, N-1.

Expresando esta definición en Python tendremos:

def fft(X):
    n = len(X)
    assert n & (n - 1) == 0 # n must be a power of two
    if n == 1:
        return X
    even_fft = fft(X[0::2])
    odd_fft = fft(X[1::2])
    w = math.e ** (-2.0 * math.pi * 1j / n)
    first_half = [ke + (w ** i ) * ko
                  for i, ke, ko in zip(itertools.count(),
                                       even_fft, odd_fft)]
    second_half = [ke - (w ** i ) * ko
                   for i, ke, ko in zip(itertools.count(),
                                        even_fft, odd_fft)]
    return first_half + second_half

donde el signo “-” entre las partes par e impar en la segunda mitad es debido a que

e^{-\frac{2 \pi i}{N} k} = e^{-\frac{2 \pi i}{2 M} (k - M)} e^{-\frac{2 \pi i}{2 M} M} = e^{-\frac{2 \pi i}{2 M} (k - M)} e^{-\pi i} = -e^{-\frac{2 \pi i}{2 M} (k - M)}.

Obviamente esta implementación en Python no es competitiva con FFTW 😀 pero se observan claramente las diferencias que el comportamiento asintótico hace en la performance de los algoritmos. En este caso tenemos una complejidad total lineal en cada nivel de recursión, dando una complejidad general \Theta(n \log n).

Sería interesante comparar la performance entre una implementación en C del algoritmo DFT básico y esta implementación en Python, de modo de ver para que valor de n el comportamiento asintótico pasa a dominar sobre las constantes…

Terrenos multifractales – estado actual

Este es un screenshot del estado actual de la demo:

Screenshot de la demo de terrenos multifractales en su estado actual.

Screenshot de la demo de terrenos multifractales en su estado actual.

El terreno se ve bastante poco atractivo por la carencia absoluta de texturas e iluminación y existen algunos problemas en la forma de las costas, pero parece un camino prometedor. prácticamente no sufrió cambios en el aspecto desde hace un par de semanas, ya que me concentré en migrar desde SDL y GLEW a GLFW y GLEE (en busca de mayor simplicidad y menos DLLs).

El número de la barra de títulos es el número de FPS, pero no es tan bajo por la ineficiencia de mi código (que es mucha, sin lugar a dudas, pero compensada por la placa de video) sino porque no logro descubrir como desactivar el vsync :-S

“Portal” y otros temas relacionados

Hace una semana terminé a Portal, uno de los pocos juegos de los que puedo decir eso (la mayoría me aburren rápido y nunca los termino).

Puede verse el concepto claramente en el trailer: consta de utilizar un “arma” que crea portales, que son básicamente conexiones entre puntos del espacio. Por razones de practicidad y jugabilidad solo se pueden crear en ciertas superficies, pero debo decir es la primera vez que veo un juego basado en espacios que no sean simplemente conexos.  (Nota: bueno, en realidad el Asteroids se jugaba en un espacio toroidal… digamos que es el primer juego que conozco basado en espacios tridimensionales múltiplemente conexos.)

En la realidad se desconoce si la creación de portales o “agujeros de gusano” (que siempre me pareció una traducción poco feliz del término wormhole) es posible en principio y, como es obvio, no tenemos la menor idea de como crearlos con tecnología remotamente a nuestro alcance. Pero se han realizado estudios teóricos bastante detallados de sus propiedades tales como Lorentzian Wormholes: From Einstein to Hawking, llevando a propuestas de buscar objetos con esas características (aunque desconozco si fue realizada alguna búsqueda).

La existencia de estos objetos en nuestro universo me parece algo dudosa, pero ciertamente ayudarían a realizar algunos de los sueños de la ciencia ficción 😀

Extra: para los que no piensan jugar al Portal o no les molestan los spoilers, los credits del juego en HQ (la voz pertenece a una computadora émula de HAL 9000).