Exploraciones Matemáticas con Scilab y Linux

Por C.E Pramode

Traducción al español por José Gregorio del Sol Cobos
el día 12 de Abril de 2004, para La Gaceta de Linux

Poco se podría haber imaginado el científico y revolucionario francés Jean Baptiste Joseph Fourier, que las técnicas analíticas que había inventado para estudiar el comportamiento de las funciones matemáticas se convertirían algún día en las herramientas más poderosas en las manos de científicos e ingenieros dedicados a disciplinas tan diversas como la neurofisiología y la comunicación digital.

Dado que me estaba deslizando rápidamente en las profundidades de la ignorancia matemática, pensé que podría refrescar algunos recuerdos del instituto intentando comprender un poco de la matemática de Fourier. Mucho de lo que leí voló bien por encima de mi cabeza - mi único consuelo fue que descubrí en Limux una plataforma ideal no sólo para el hackeo de Sistemas Operativos, sino también para el recreo y la investigación matemáticos.

Me encontré con una gran herramienta llamada Scilab y también con un bonito manual sobre las matemáticas de Fourier de Chris Meyers que demostraba algo de material interesante sobre el análisis y la combinación de ondas en seno empleando código Python. Este artículo demuestra algunos sencillos trucos de Scilab y reimplemnenta el código de Chris en el lenguaje de guiones nativo de Scilab. ¡Se advierte a los lectores que busquen rigor matemático que no confíen demasiado en lo que digo aquí!

¿Qué es Scilab?

Scilab es un entorno de computación matemática poderoso y libre. Proporciona un marco de trabajo extensible para la manipulación general de matrices y "cajas de herramientas" para hacer cosas como diseño de control de sistemas, procesamiento digital de la señal, etc. El código fuente en C/Fortran está disponible para su descarga desde el sitio web del proyecto -yo no tuve ninguna dificultad en construir el sistema - el mágico "configure; make; make install" estándar funcionó perfectamente.

Aquí viene un pantallazo de Scilab ejecutándose en mi caja Linux:

Matemática sencilla

Empecemos haciendo unas sencillas manipulaciones con matrices. Se crea una matriz 3 por 3 simplemente escribiendo, en la consola de Scilab:

-->a = [1,10,20; 5,6,7; 12,11,45]
 a  =
 
!   1.     10.    20. !
!   5.     6.     7.  !
!   12.    11.    45. !
 
-->

Es fácil obetener la matriz transpuesta:

--->a'
ans  =
 
!   1.     5.    12. !
!   10.    6.    11. !
!   20.    7.    45. !
 
-->

Algunas otras funciones:

-->sum(a, 'c')
 ans  =
 
!   31. !
!   18. !
!   68. !
 
-->sum(a, 'r')
 ans  =
 
!   18.    27.    72. !
 
-->diag(a)
 ans  =
 
!   1.  !
!   6.  !
!   45. !
 
-->

Los elementos de las matrices se pueden extraer de muchos modos diferentes - el más sencillo es el procedimiento tradicional del indexado. Escribir a(1,2) daría el elemento en la fila 1 y la columna 2 (observe que el índice empieza en 1). Indexar una matriz más allá de su límite dará un error. Escribir un índice no existente resultará en que la matriz crecera dinámicamente.

-->a(3,4) = 3
 a  =
 
!   1.     10.    20.    0. !
!   5.     6.     7.     0. !
!   12.    11.    45.    3. !
 
-->

El operador "dos puntos"

El operador "dos puntos" es un bonito operador. Podemos crear un vector de números, 1, 2, 3, ... 10 simplemente escribiendo:

-->a = 1:10
 a  =
 
!   1.    2.    3.    4.    5.    6.    7.    8.    9.    10. !
 
-->

Otros muchos trucos son posibles:

-->b
 b  =
 
!   1.    2.    3. !
!   4.    5.    6. !
!   7.    8.    9. !
 

-->b(1:3,2:3)
ans  =
 
!   2.    3. !
!   5.    6. !
!   8.    9. !
 
-->1:2:10
 ans  =
 
!   1.    3.    5.    7.    9. !
 
-->

Observe que 1:2:10 significa crear un vector que empiece en 1, con cada elemento siguiente calculado añadiendo 2, hasta que el valor se hace mayor que 10.

Dibujo sencillo

Miremos un ejemplo del dibujo de una sencilla onda en seno. Queremos un ciclo completo de la curva seno (desde 0 hasta 2*PI) -tomemos 240 puntos en medio, para que cada división sea 2*PI/240. Primero creamos un vector que contenga todos los valores del ángulo en ese rango y después lo representamos (%pi es una constante con el valor de PI):

--> = 0:(2*%pi)/240:2*%pi
 x  =
         column 1 to 5
 
!   0.    0.0261799    0.0523599    0.0785398    0.1047198 !
 
         column 6 to 9
 
!   0.1308997    0.1570796    0.1832596    0.2094395 !
 
         column 10 to 13
 
!   0.2356194    0.2617994    0.2879793    0.3141593 !
 
         column 14 to 17
 
!   0.3403392    0.3665191    0.3926991    0.4188790 !
[More (y or n ) ?] 

Ahora, empleamos una sencilla función de dibujo:

-->plot(x, sin(x))

Escribiendo guiones de Scilab

Escribir guiones de Scilab es sencillo. He aquí un ejemplo de un bucle "for" que se puede ingresar en la propia consola de Scilab:

-->s = 0         
 s  =
 
    0.  
 
-->for i=1:3:10
-->  s = s + i
-->end
 s  =
 
    1.  
 s  =
 
    5.  
 s  =
 
    12.  
 s  =
 
    22.  
[More (y or n ) ?] 
 

Definiendo funciones

La sintaxis de la definición de funciones tiene su aquél. He aquí un ejemplo sencillo:

-->function [r] = my_sqr(x)
-->  r = x * x
-->endfunction
 
-->my_sqr(3)
 ans  =
 
    9.  
-->

Después de la palabra clave "function", proporcionamos una lista de "valores de salida". Cualquier valor escrito a un valor de "salida" será "devuelto" por la función. El argumento "x" es por supuesto el argumento de entrada de la función. La función devuelve el valor "r" que es el cuadrado de "x".

Obviamente la cuestión es qué ocurre si queremos que devuelva dos valores. Ensayamos lo siguiente en la consola de Scilab:

-->function [r1, r2] = foo (x, y)
-->  r1 = x + y
-->  r2 = x - y
-->endfunction
 
-->[p, q] = foo(10, 20)
 q  =
 
  - 10.  
 p  =
 
    30.  
 
-->

Observe el modo especial de llamar a la función. El valor de r1 se transferirá a "p" y el valor de r2 a "q".

Las siguientes invocaciones de "foo" demuestran que el lenguaje está dinámicamente tipado.

-->[p, q] = foo([1,2], 1)
 q  =
 
!   0.    1. !
 p  =
 
!   2.    3. !
 
-->[p, q] = foo([1,2], [3,4,5])
 !--error     8 
inconsistent addition
at line       2 of function foo   called by :  
[p, q] = foo([1,2], [3,4,5])
 
-->

Es posible almacenar definiciones de funciones en un archivo y cargarlas en un instante posterior. Suponga que la definición de la funcion de arriba está almacenada en un archivo llamado "fun.sci". Simplemente necesitamos invocar, en la consola de Scilab:

-->exec('fun.sci')

¡Entre, Fourier!

Nosotros encontramos "señales" por todas partes. El micrófono del PC genera sonido convirtiendo señales eléctricas a vibraciones. Vemos objetos a nuestro alrededor porque esos objetos reflejan las señales de la luz hacia nuestros ojos. Nuestra radio y televisión reciben señales electromagnéticas. ¡Estamos inmersos en un "mar de señales"! El análisis de las señales es pues de central importancia en muchas ramas de la ciencia y la ingeniería.

La filosofía básica de Unix es "Manténlo Sencillo, estúpido". Los físicos (y muchos otros científicos e ingenieros) no pueden a menudo adherirse a ese dicho cuando empiezan a analizar algo, simplemente porque los fenómenos que están estudiando tienen una fea complejidad. Pero parece que las cosas más complejas de este mundo se pueden explicar en base a cosas más sencillas. El pesamiento de Joseph Fourier fue que las señales complejas y variables con el tiempo se pueden expresar como combinación de curvas sen/cos simples de frecuencia y amplitud variables. Verificaremos esta suposición representando algunas simples ecuaciones con la ayuda de Scilab.

Comencemos con una sencilla suma de dos señales seno:

-->delta = (2*%pi)/240
 delta  =
 
    0.0261799  

-->x = 0:delta:2*%pi

-->a = sin(x) - (1/2)*sin(2*x)
-->plot(x, a)
 

Aquí está el dibujo:

Hay pocos indicios aquí de que esté a punto de ocurrir algo interesante. Seguidamente intentamos representar

b = sin(x) - (1/2)*sin(2*x) + (1/3)*sin(3*x)

Si siguiéramos añadiendo términos a la serie, el siguiente sería -(1/4)*sin(4*x), el siguiente +(1/5)*sin(5*x), y así sucesivamente. He aquí lo que obtuve cuando representé esta serie con 200 términos (usted tendrá que escribir una función para que lo haga por usted):

¡Parece magia! ¡La curva seno de ha desvanecido por completo y tenemos una señal completamente nueva! El cómo exactamente el Sr. Fourier "supo" que una serie así nos daría al final algo totalmente diferente de la suma de sus partes se trataría más apropiadamente en una clase de matemáticas (¿escucho vuestro bostezo?¿Tenemos un caso de educación matemática más "práctica" con estudiantes teniendo acceso a cajas de Linux ejecutando Scilab, Python(Numeric) y un montón de otras herramientas educativas libres?).

Determinando las componentes de una señal

Hemos visto que añadiendo senos de frecuencia y amplitud diferentes nos da señales que lucen totalmente diferentes. Ahora la cuestión es, dados algunos númertos que representen una onda compleja, ¿seremos capaces de decir qué combinación (frecuencia, amplitud) de los senos dieron luigar a aquella señal particular? Intentémoslo.

Escribamos primero una función que realice la "integración" numérica sencilla en un rango entre 0 y 2*PI. Dividimos el área bajo nuestra curva en pequeñas tiras, cada una de ancchura, digamos que 2*PI/240. El área de una tira en el punto "x" (0 < x < 2*PI) será su altura multiplicada por su anchura, que será sin(x) * (2*PI/240). Ésta es la idea detrás de la función de integración, que se puede escribir en la consola de Scilab. El argumento a integrar es un vector de valores del seno entre 0 y 2*PI-delta, donde delta es (2*PI)/240. La diferencia entre dos valores sucesivos del vector es "delta".

-->function [r] = integrate(a)
-->  r = sum(a)*(2*%pi)/240
-->endfunction

Intentemos integrar la función seno sencilla, sen(x).

 
-->x = 0:delta:(2*%pi-delta)

-->integrate(sin(x))
 ans  =
 
    3.837E-16  
 

Vemos que la integral es cero. La curva seno tiene igual área por encima y por debajo del punto cero.

Intentemos representar sin(x).*(-sin(x)). (Observe que el operador .* realiza la multiplicación de dos vectores miembro a miembro):

Vemos que toda la función se ha desplazado por debajo del punto cero. Ahora debería tener un área definitivamente distinta de cero.

-->integrate(sin(x).*(-sin(x))) 
 ans  =
 
  - 3.1415927  
 
-->

Scilab nos dice que esto es -PI. Intentemos ahora representar sin(2*x).*(-sin(x)).

La gráfica nos dice que la integral debería ser cero. Lo verificamos:

-->integrate(sin(2*x).*(-sin(x)))  
 ans  =
 
    3.977E-16  

Estamos empezando a coger el "sentido" de la idea que emplearíamos para separar las componentes de nuestra señal compleja. Multiplicando un seno con el opuesto de un seno de diferente frecuencia nos da cero -sólo cuando las frecuencias coinciden obtenemos resultados no nulos. Digamos que nuestra compleja señal es:

sin(x) - (1/2)*sin(2*x) + (1/3)*sin(3*x) - (1/5)*sin(5*x)

Si multiplicamos esto por -sen(x), lo que obtenemos es:

sin(x).*(-sin(x)) - (1/2)*sin(2*x).*(-sin(x)) + 
(1/3)*sin(3*x).*(-sin(x)) - (1/5)*sin(5*x).*(-sin(x))

El primer términos nos da -PI, todos los demás términos se hacen cero. El hecho de que estemos obteniendo un valor no nulo nos dice que sen(x) está presente en la señal. Ahora multiplicamos la señal por -sen(2x). Si obtenemos un resultado no nulo, eso significa que sen(2*x) está presente en la señal. Repetimos este proceso tantas veces como queramos.

¿Cómo obtenemos la amplitud de cada componente? Intentemos otro experimento:

 
-->b = sin(x) - (1/2)*sin(2*x) + (1/3)*sin(3*x) - (1/4)*sin(4*x)
 
-->integrate(b.*(-sin(x)))                                      
 ans  =
 
  - 3.1415927  
 
-->integrate(b.*(-sin(2*x)))
 ans  =
 
    1.5707963  
 
-->integrate(b.*(-sin(3*x)))
 ans  =
 
  - 1.0471976  
 
-->integrate(b.*(-sin(4*x)))
 ans  =
 
    0.7853982  
 

Vemos que dividiendo cada resultado entre -PI nos da la amplitud de cada componente de la señal.

Conclusión

Existen herramientas propietarias de muy alta calidad para hacer matemática simbólica y numérica, pero a veces tienen un precio más allá del alcance del estudiante o el aficionado. Espero que este artículo le haya convencido de que las alternativas dentro del Software Libre existen realmente. Hágame saber amablemente cualquier inexactitud que encuentre en este documento. Se me puede contactar a través de mi página web en pramode.net.

Agradecimientos

Gracias al equipo de Scilab por crear dicha poderosa herramienta y también por documentarla tan ampliamente. Este artículo no se habría escrito sin la ayuda del documento de Chris Meyers explicando la matemática de Fourier -Chris también ha escrito algunos otros interesantes programas en Python que estoy seguro usted disfrutará. ¡Un gran Gracias a él!

 

Soy un profesor trabajando en el IC Software de Kerala, India. Habría adorado convertirme en químico orgánico, pero hago la segunda mejor cosa posible que es ¡jugar con Linux y enseñar a programar!