Electrónica y Programación en Personal
«Si no se soluciona con un script en Python entonces no es viable»

por Yeison Cardona el 18 de febrero del 2013 a las 2:25 am UTC
La finalidad de este post no es hacer las demostraciones matemáticas que sustentan este concepto, sino, una orientación de cómo hacer una simulación en Python con ayuda de Matplotlib y Numpy

Teoría (si se le puede llamar así)

Es fácil demostrar que la siguiente fórmula para funciones continuas:



La cual corresponde a la representación exponencial de la serie de Fourier para cualquier señal de energía.
Donde cada una de los coeficientes Xn de la serie son determinados por la siguiente expresión:




Para una correcta interpretación teórica recomiendo el libro de mi profesor Análisis de Señales y Sistemas, y como nuevamente no estoy seguro de las licencias, por favor consideren adquirir la versión impresa del libro, y usen esta versión en PDF sólo como sustento temporal.

Objetivo

Entonces nuestra tarea es desarrollar un código en Python que nos permita trabajar con cualquier señal discreta, entonces es necesario reescribir la expresión ya que una sumatoria infinita en sistemas computacionales resulta imposible, también se tienen que reescribir las integrales (y sus límites) en términos de sumatorias, con lo que obtenemos:



Lo que hemos reescrito aquí es todo lo necesario para hacer la simulación:
  • x(t) será nuestra función discreta, un array de numpy.
  • x̂(t) será la señal aproximada
  • N será la dimensión de conjunto ortogonal que evaluaremos.
  • tf, ti será el tiempo final y el inicial de nuestra señal, si suponemos que nuestra señal iniciará en 0 entonces tf será el tamaño del array, la cantidad de datos con la que estamos trabajando.
  • ω corresponde a una variable clave para hacer que conjunto ortogonal de funciones sea ortogonal en el dominio de nuestro muestreo.

Simulación

Dependencias

Serán necesarias algunas librerías extras, las cuales ya hemos usado anteriormente en otras entradas: Numpy y Matplotlib.
sudo apt-get install numpy
sudo apt-get install matplotlib

Aproximación de señal de audio

La señal de audio que usaremos constará de 96 elementos.
#!/usr/bin/env python
#-*- coding: utf-8 -*-

from numpy import e, pi, array, nan
from matplotlib import pyplot as plt
import numpy, random

# Generar el array con la señal a partir de un archivo .wav
signal_input = numpy.fromfile("fairy__.wav", long)

# Tf-Ti
T = len(signal_input)

# Generar el conjunto ortogonal de funciones
Qn = lambda n: lambda t: e**(-1j*n*t*(2*pi)/T)

# Armónicos en términos de n
Xn = lambda n: (1.0/T) * sum(map(lambda t:signal_input[t]*Qn(n)(t), range(T)))

# Función aproximada
Xt = lambda N: map( lambda t: sum(map(lambda n: Xn(n)*e**(1j*n*(2*pi)/T*t), range(-N, N))), range(T))

# Error
E0 = lambda N: (1.0/T) * sum((signal_input-Xt(N))**2)
De éste código hay que resaltar la simplicidad, una expresión para cada componente de análisis, nos hemos valido de funciones anónimas y algunas herramientas para programación funcional que nos provee Python, es exactamente lo que está definido en el último conjunto de expresiones.

Ahora, puede que el código sea un poco engañoso, pero en realidad lo que se tiene es un conjunto de funciones que podemos usar para muchas cosas, por ejemplo, plotear la señal y su aproximación para un numero de funciones dentro del conjunto ortogonal:

#----------------------------------------------------------------------
def plot(N):
    plt.plot(range(T), signal_input)
    plt.plot(range(T), Xt(N))
    plt.legend(("Signal", "Apprx, n=%d"%N))
    plt.grid()
    plt.show()

# 5 funciones ortogonales
plot(5)

# 10 funciones ortogonales
plot(10)

# 20 funciones ortogonales
plot(20)
Y obtenemos:



Como se puede ver, cada vez que usamos un conjunto mayor de funciones ortogonales, obtenemos una mejor aproximación de la señal, ésto para funciones continuas debe tender a infinito, es decir, un conjunto infinito de funciones ortogonales para que la función sea una representación exacta a la original, pero en funciones discretas, la aproximación exacta se alcanza cuando el número de elementos del conjunto ortogonal el igual al número de elementos de la función (o señal).

Aproximación de función cuadrada

Si usamos una función cuadrada es más fácil de ver como la señal se va acercando mas a la función original cada que aumentamos en conjunto ortogonal, para definir una función cuadrada:
#----------------------------------------------------------------------
def rect(t):
    if t < 50: return 11
    else: return -11
        
entonces debemos cambiar la señal de entrada por:
signal_input = map(rect, range(100))
Ahora, plotearemos varias señales en una misma gráfica, y así veremos de lo que estamos hablando:
#----------------------------------------------------------------------
def multiPlot(min, max):
    plt.plot(range(T), signal_input)
    legend = []
    for i in range(min, max+1):
        plt.plot(range(T), Xt(i))
        legend.append("n=%d"%i)
        print i
    plt.legend(legend)
    plt.grid()
    plt.show()

# Plotearemos desde 2 hasta 6 elementos por conjunto ortogonales
multiPlot(2,6)

y obtenemos:
y si se sigue aumentando:

Error

Dado que el error es acumulativo, entonces se necesita de una forma para calcular el error, este error se define como el error cuadrático medio y está dada por la siguiente expresión:

Esta función ya se ha definido al inicio del código (E0), y para gráficar una una relación error contra conjunto ortogonal:

E = 46
plt.plot(range(1, E), map(lambda n: E0(n), range(1, E)))
plt.legend(("Error", ))
plt.grid()
plt.show()

y vemos algo como:
Si se sigue calculando veremos como el error crece cuando se usan conjuntos ortogonales mayores a la cantidad de elementos de la señal (se demorará su tiempo):

Descarga

El código python completo, incluyendo el archivo de audio que se ha usado en esta entrada se puede descargar desde éste link.

Conclusiones

¿Quien necesita de Matlab cuando se tiene Python?

También podría interesarte:

Añadir un comentario:
Si desean una respuesta para su comentario sólo deben agregarme en G+ y hacer una mención a Yeison Cardona, así les podré responder lo antes posible.