domingo, 7 de agosto de 2016

Tonantzintla

Epicuro afirmaba que lo único que hace verdaderamente feliz al ser humano es vivir y trabajar en comunidad con sus amigos. Creo que cualquiera que haya pasado una escuela, taller o verano de investigación con buena compañía le puede dar la razón. Lo único malo, como le oí decir a alguien durante la última semana del taller, es el "síndrome de abstinencia post-escuela". El forjado de experiencias de vida es una de las cosas más personalmente enriquecedores que puede ofrecernos la ciencia. No podría dedicarme a otra cosa.

Debo admitir que cuando mandé mi solicitud estaba seguro con toda mi arrogancia que iba quedar admitido. Si bien no me equivoqué, me di cuenta al investigar a mis compañeros que no era mejor candidato de latinomérica que digamos, algo que en realidad me animó mucho. Por primera vez me tocó un evento con una selección impecable. Si alguien ha seguido esta serie de escritos ya habrá leído alguna que otra desventura. Pero esta es, desde aquel TCJ 2007, fue una de las mejores escuelas de verano en las que he estado. Había comentado alguna vez que no parece apetecerme mucho escribir cuando las cosas andan muy bien. Creo que ha sido el caso estos últimos días, pero vale la pena contar estas cosas.

Después de S. llegué a pensar que sería muy difícil sentir algo por alguien con la misma fuerza. Hubo otras chicas el año pasado, pero siendo sincero, sólo fueron para mi una salida para sacármela de la mente. Este año no tenía ningún plan no académico. Hasta que conocí a una de las mujeres más perfectas que podría encontrar. A primera vista me pareció muy guapa, pero eso jamás me llevado a sentir nada más que atracción. Supongo que fue a la mitad de la primer semana cuando me di cuanta de que estaba perdido por ella. Traté de evitarlo lo más que pude. Es una extranjera que difícilmente volvería a ver con la frecuencia que quisiera. Pero terminé completamente loco por ella, como no me hacia sentido en años. ¿Por qué? Porque simplemente era la mezcla de todas las mujeres que me habían gustado hasta el presente y mucho más. Un chica extrovertida, sin ser fastidiosa, inteligente y más lógica de lo que había visto antes. Apasionada y comprometida con su carrera. El tipo de mujer que te convence de que la perfección en las personas se puede lograr. ¿Cómo iba a poder resistirme? Su nombre es N.

La segunda semana del taller consistió en trabajar con un proyecto individual. Fue muy interesante y ha sido lo más cercano a lo que he estado a la investigación real. Trabajé con datos de Herschel y GTM Aprendí bastante sobre el WCS y sobre como convertir entre indicies y coordenadas ecuatoriales. Esa semana vi muy poco a N. Todos estuvimos muy ocupados con nuestro proyectos en las oficinas de nuestros asesores. A diferencia de Morelia, nadie se quejaba por la carga de trabajo. Fue todo lo contrario, como la Legión Extranjera, eramos casi los primeros en llegar y los últimos en irnos del instituto. Just the way I like it. Hice muy buenos amigos de los que espero ser colegas en un futuro no muy lejano.

La escuela terminó un viernes, y el sábado todo el mundo partió del instituto. N. regresaba a su país el lunes y aún tenía planes de dar la vuelta por la CDMX. Yo no tenía más que unos 200 pesos en la cuenta, pero sin pensarlo dos veces, acepté dar un curso en la capital dos días antes. Era mi última oportunidad de verla de nuevo en mucho tiempo. La vi el domingo junto con uno de sus amigos, que me cayó muy bien y terminó siendo mi compañero de cuarto en la aventura que escribiré en la siguiente entrada (espero). La tarde se pasó muy rápido para mi. No pude decirle todo lo que me hubiera gustado... Sin embargo, mientras regresaba, solo y suspirando en el metro me decidí a despedirme y declararme en el aeropuerto con en las películas. Y así fue. Después de hacer cálculos y revisar el mapa de la ciudad, ahí estaba yo a las 5:00 am en metro Universidad, en el primer tren del día. Casi una hora y media después me encontraba corriendo en la Av. FAM hacia la terminal 2. Por la hora ya había dado por perdido mi cometido pero no me detuve hasta llegar. Para mi sorpresa, ¡aún seguía ahí!, estaba dormida en las bancas. Tome algunos minutos para recobrar el aliento, ordenar mi mente y esperar a que despertara. Cuando lo hizo, me acerque a saludarla. Me hizo la pregunta que estaba esperando, que era lo que hacía ahí. Todo el discurso que tenía se me borró por completo. Supongo que notó mi sonrisa estúpida sin palabras y agregó sonriendo "vino a despedirse, ¿verdad? Vamos a sentarnos, aún tengo tiempo". Unos 20 minutos después la acompañe a la entrada de la sala de abordar. Tan pronto la vi cruzar la puerta hacia la aduana sentí un extraño nudo en la garganta. En ese momento le envié un mensaje que hasta ahora sigo preguntádome si leyó, debido a un problema que ocurrió con su vuelo que pudo haberlo perdido en medio del apuro. Sólo he llorado 3 veces en un aeropuerto. Ver partir a alguien que me hizo sentir algo en tan poco tiempo fue motivo de esa tercera.

miércoles, 20 de julio de 2016

Concepto de frecuecia en tiempo discreto

Este es un tema que si no se entiende bien al principio de un curso de procesamiento digital de señales puede acarrear muchos problemas al ir avanzando en temas más complicados. En esta entrada trataré de explicar lo más claro posible lo que es la frecuencia de una señal discreta. Primero, comencemos describiendo a una función cosenoidal en tiempo continuo con la notación del Proakis [1]:
El subindice a se utiliza para denotar que x_a(t) es una señal analógica y no confundirla con su contraparte discreta x[n]. Ahora, tomemos esta misma señal discreta y multipliquémosla por tren de pulsos unitarios con un periodo de muestreo Ts (proceso conocido como muestreo).  Gráficamente, el resulto obtenido es el siguiente:



























Tenemos ahora una nueva función que depende de una variable entera n multiplicada por el periodo de muestreo Ts. Es este escalamiento por Ts lo que permite que la señal discreta x[n] coincida en el tiempo con la señal x_a(t). x[n] es una nueva señal que está definida en función de una secuencia de números enteros que representan una posición en un arreglo de datos (estas posiciones pueden tomar valores negativos si se toman con respecto a un origen arbitrario en el arreglo):
Como consecuencia del teorema de muestreo de Nyquist-Shannon, sólo es posible obtener información sobre una señal o componente a una frecuencia dada sin riesgo de ambigüedades (conflictos con sus alias) mientras su frecuencia sea menor o igual a la mitad de la frecuencia de muestreo (que es simplemente el inverso de su periodo de muestreo Fs = 1/Ts). Por lo tanto, la frecuencia relativa (también llamada frecuencia normalizada) esta limitada al intervalo:
Si se grafica un ejemplo numérico, resulta bastante intuitivo visualizar que cuando un coseno discreto tiene una frecuencia relativa de pi radianes por muestra se alcanza el periodo mínimo distinguible para una función periódica discreta. Los sinusoides en tiempo discreto muestran periodicidad en frecuencia. Esto significa que los sinusoides discretos que guarden una diferencia en frecuencia de 2*pi son idénticos. Esto puede verse más claramente en esta animación que hice en Matlab:
Naturalmente, lo que nos interesan son los parámetros físicos a los que las señales en tiempo discreto hacen referencia. El periodo de muestreo Ts es el puente clave para darle significado en tiempo continuo a estas señales. Resumiendo las transformaciones de frecuencia:
 Referencias:
Sec. 1.3 "Concepto de frecuencia en señales continuas y discretas en el tiempo",  "Tratamiento Digital de Señales", John G. Proakis
Teorema del Muestreo, L. Javier Morales Mendoza, UGto

martes, 12 de julio de 2016

¿Qué tan lejos están las estrellas que podemos ver a simple vista?


Hace unos años se había abierto una discusión en reddit (también en Quora)a causa de la tira de xkcd de arriba. Hoy me encontré con una discusión similar ahora por un post en Facebook. Como ni las respuestas en reddit ni Quora me convencieron, decidí encontrar la respuesta yo mismo analizando algunos datos.

Primero, descargué la base de datos HYG (Hipparcos, Yale Bright Star, and Gliese) desde aquí, que nos da la información combinada de los 3 catálogos estelares mencionados, una muestra de cerca de 120,000 estrellas. Teniendo el archivo, escribí un programa en Python usando la librería Pandas [pueden bajarlo aquí] para analizar los datos. Y esto fue lo que encontré para la todas las estrellas resueltas [que pueden distinguirse como fuentes puntuales individuales] con una magnitud menor a 6:

De acuerdo al catálogo, el 98% de las estrellas que podemos ver a simple vista en la Tierra (muy lejos de las ciudades y sin luna) se encuentran a menos de 1000 años luz. Aún más interesante, las distancias de estas estrellas tienen una media de sólo 145 años luz con una desviación estándar de 136 años luz . Considerando los porcentajes en una distribución normal:
estamos hablando de que el 95% de esas estrellas se encuentran a menos de 417 años luz. Debe notarse también de que el 2% de las estrellas que se encuentran a más de 1000 años luz en este catalogo tienen una incertidumbre considerable y todas ellas se redondean a 10,000 años luz. En conclusión, cuando uno mira una estrella en la noche, lo más probable es que esté a menos de 400 años luz de distancia.

miércoles, 29 de junio de 2016

Convolución de dos señales en Matlab

La convolución continua de dos funciones esta dada por:
De esta forma, lo que esperaríamos ver al convolucionar dos pulsos cuadrados de amplitud y ancho unitario se muestra en la siguiente animación:

 Para implementar la convolución de dos señales en Matlab se debe considerar que la función conv() [ver documentación] está realizando la operación de convolución discreta que está expresada como:
Si uno utiliza esta función esperando los resultados de una convolución continua, se lleva la sorpresa de que los valores de la convolución resultante son muchos más altos que los valores correctos. Esto se ocurre porque hace falta considerar que una convolución continua se aproxima como el producto del periodo de muestreo por la convolución discreta de las mismas señales evaluadas en tiempo discreto. Esto es:
Así, el programa en Matlab nos queda:

Ts = 0.01;
t = -2:Ts:2;
f = (t>-0.5) - (t>0.5);
g = f;

%Convolución
cnv = Ts*conv(f,g,'same');

subplot(3,1,1),plot(t,f), ylim([0 1.5]),title('f(t)'),grid on
subplot(3,1,2),plot(t,g,'r'), ylim([0 1.5]),title('g(t)'),grid on
subplot(3,1,3),plot(t,cnv,'g'), ylim([0 1.5]),title('conv(f(t),g(t))'),grid on

Referencias:
Using MATLAB with the Convolution Method, California State University Northridge

miércoles, 1 de junio de 2016

Language Barriers

I hadn't written in English here before. Now, regarding that I need to take the GRE this year, it's something that I enormously regret. Mostly what I post here are technical notes of problems I solved and that there is no documentation online about it in Spanish. Someone commented once that I should write my technical notes in English so people all over the world can use them. What I always say is that there's plenty of resources written in English while good ones written in Spanish are scarce. That's why I prefer to write in my mother language. I feel I'm contributing to the develop of my region, not only Mexico but all the Hispanic world. But this blog is not only about engineering and science stuff (things I really passionate about), although it seems to be the case. It's about my life too, the whole of it. I think writing about ourself and our lives is not narcissism. We need to know other peoples thoughts and passions in order to make a better view of the world. So I know that writing always in Spanish have trapped me in bubble too. I haven't done blogger friends outside Hispanic America. But the blogosphere is kind of dead these day anyway, isn't it? Or may be it's just in the Spanish Internet! I don't even know it for the same reason. Well, I'll give it shot.

sábado, 14 de mayo de 2016

Duloxetina

A Quay at Caen - Stanislas Lépine
Algo muy extraño que ocurrió cuando tomaba Cymbalta fue perder sólo durante esos días el gusto por las pinturas oscuras de Stanislas Lépine, Christian Dahl y Atkinson Grimshaw. No es un suceso que me haya caído en completa sorpresa, pero me parece muy curioso que el gusto en arte esté tan directamente condicionado por la química cerebral. Otro efecto interesante fue perder temporalmente la capacidad de acceder a memorias asociadas con ansiedad o tristeza. ¿Es este mismo mecanismo el mismo que me alejaba del gusto por las pinturas oscuras? Tan pronto dejé el antidepresivo, hace ya algunas semanas, mi gusto volvió. No lo dejé por orden médica ni porque no me sintiera cómodo con los resultados, sino porque me incomodaba el hecho de de tener que vivir con una tranquilidad artificial.

Atkinson Grimshaw
Debo aclarar que no padezco de depresión sino que he sufrido de un desorden de ansiedad desde que tengo memoria. Un problema que mis psicólogos no supieron identificar hasta que no pude soportarlo más y tuve que visitar a un psiquiatra. Ahora lo único que lamento es no haberme tratado este problema antes. Pienso en todo el tiempo que la ansiedad me arrebató de mi vida sin poder impedirlo a pesar de mis intentos. En el fondo sabía que en mí había problema más grande que la timidez. Había sido como una batalla en un cuarto competentemente oscuro. Tarde, si, pero finalmente me siento verdaderamente bien conmigo mismo, como no lo sentía desde niño. Sé ahora que esa presión en el pecho y esa angustia constante que me envolvía no eran parte de mí. He aprendido tanto en estas últimas semanas.  

miércoles, 27 de abril de 2016

Leer imagenes en formato .img en Matlab

Es muy común que datos científicos como imágenes medicas o astronómicas se encuentren con la extensión .img. Al menos hasta la versión de 2016, Matlab no tiene una función para extraerlas pero en esta entrada veremos como hacerlo leyendo directamente los datos binarios con fopen(). Para este ejemplo usaremos los datos topográficos del altímetro láser de la sonda espacial Mars global Surveyor que pueden ser descargados aquí. Es muy importante saber de antemano lo siguiente::
  • Dimensiones de la imagen 
  • Formato y precisión de los datos
Esta información es necesaria para poder redimensionar y obtener correctamente el array de datos leídos con la función fread(). En este caso usaremos el mapa topográfico de la superficie de Marte a una resolución de 4 pixeles por grado (archivo megt90n000cb.img). Estos datos tienen un archivo de documentación asociado llamado megt90n000cb.lbl donde podemos leer que la imagen tiene dimensiones de 1440 por 720 px con datos en formato big-endian int16. Sabiendo esto escribimos el siguiente programa:

%Abrir archivo .img
fid = fopen('megt90n000cb.img','rb');
dim = [1440 720];
img = fread(fid,dim,'int16','ieee-be');
img = img';

%Gráfica 2D

figure(1)
imagesc(img)
xticklabels = 0:30:360;
xticks = linspace(1, size(img, 2), numel(xticklabels));
set(gca, 'XTick', xticks, 'XTickLabel', xticklabels)
yticklabels = -90:30:90;
yticks = linspace(1, size(img, 1), numel(yticklabels));
set(gca, 'YTick', yticks, 'YTickLabel', flipud(yticklabels(:)))
colormap(hot), grid('on'), title('Mars Orbiter Laser Altimeter Data')
xlabel('Longitud(E)'),ylabel('Latitud')

%Gráfica 3D (Valle Marineris)
figure(2)
marineris = img(325:478,1079:1322);
surf(marineris,'EdgeColor','none')
xlim([0 200]), ylim([0 150])
title('Valles Marineris')


Plots:


martes, 26 de abril de 2016

Máquina de Estados Finitos en Python

Una máquina de estado o autómata finito (Finite State Machine en inglés) es un modelo computacional bastante útil en robótica y automatización. Para este ejemplo haremos una implementación por software de una máquina de estados en Python que puede ser fácilmente adaptada para algún proyecto en una Raspberry Pi. Esta implementación es bastante sencilla, para FSM's más robustas basadas en clases recomiendo checar aquí y aquí. El diagrama de estados para FSM de esté ejemplo será el siguiente:
Este tipo particular de FSM es llamada Máquina de Moore ya que la transición de estado sólo depende del estado actual y del valor de la entrada. Esta implemntación en Python aprovechará una característica peculiar de los diccionarios que permite asociar una función con una llave [ref]. De esta manera los estados quedan definidos como funciones dónde podemos incluir todas las acciones que no interese realizar durante dicho estado (o durante las transiciones). El programa completo:

# -*- coding: utf-8 -*-
"""
Created on Thu Feb  4 23:36:28 2016

@author: Rodolfo
"""

from time import sleep
from random import randint

#Variable global
estado = 'i'

#Estados
def EDOi(entrada):
    global estado
    print('Estado Inicial')
    #Transiciónes
    sleep(2)
    if entrada == 0:
        estado = 'i'
    if entrada == 1:
        estado = 0
        print(u'Transición hacia 0...')

def EDO0(entrada):
    global estado
    print('Estado 0')
    #Transiciones
    sleep(2)
    if entrada == 0:
        estado = 1
        print(u'Transición hacia 1...')
    if entrada == 1:
        estado = 2
        print(u'Transición hacia 2...')
        
def EDO1(entrada):
    global estado
    print('Estado 1')
    #Transiciones
    sleep(2)
    if entrada == 0:
        estado = 0
        print(u'Transición hacia 0...')
    if entrada == 1:
        estado = 2
        print(u'Transición hacia 2...')

def EDO2(entrada):
    global estado
    print('Estado 2')
    #Transiciones
    sleep(2)
    if entrada == 0:
        estado = 1
        print(u'Transición hacia 1...')
    if entrada == 1:
        estado = 0
        print(u'Transición hacia 0...')

#Finite State Machine (FSM)   
def FSM(entrada):
    global estado
    switch = {
       'i':EDOi,
        0 :EDO0,
        1 :EDO1,
        2 :EDO2,
    }
    func = switch.get(estado, lambda: None)
    return func(entrada)

#Programa Principal
while True:    
 FSM(randint(0,1))
 sleep(2)


Salida del programa en la consola en Raspberry:

domingo, 10 de abril de 2016

Streaming de datos por TCP en Raspberry Pi + LabVIEW

 Si se desea hacer un streaming de datos desde un sistema de instrumentación remota como una estación meteorológica, sismológica o un telescopio robótico, puede resultar muy útil y económico tener un servidor TCP en un Raspberry Pi y visualizar los datos en tiempo real con LabVIEW desde cualquier lugar. Esto puede hacerse muy fácilmente utilizando el modulo socket para Python.

Un socket es una abstracción de programación para la representación de conexiones y permiten realizar una comunicación bidireccional. En Python se crea de la siguiente manera:

import socket
s = socket.socket()    #Instaciamiento para el objeto 's'
host = ''                      #Host local por defecto
port = 8006                #Puerto
s.bind((host,port))     #Asignación de dirección y puerto a la instancia
s.listen(5)                  #Número máximo de conexiones entrantes 

Para este ejemplo simularemos 2 señales que usaremos para enviar por TCP. En una aplicación real estas señales serían adquiridas mediate ADC's u otros dispositivos conectados a los pines GPIO de la tarjeta:

t = np.linspace(0,2*np.pi,100)
y0 = np.float16(np.sin(t))
y1 = np.float16(0.5*np.sin(t)*np.sin(2*t))

Se puede observar que estoy usando la función float16() para reducir la precisión de los datos para ahorrarle un poco de trabajo al cliente en LabVIEW pero pueden omitirla si requieren más dígitos significativos. Ahora, para parte principal del programa tenemos:

print u"Esperando conexión"
s, addr = s.accept()
print u"Conexión desde: ", str(addr)
i = 0
try:
 while True:
     s.send('CH0'+str(y0[i])+'\n')
     s.send('CH1'+str(y1[i])+'\n')
     i += 1
     if i>=99:
         i = 0
     time.sleep(0.01)
except KeyboardInterrupt:
s.close()

 Observen que en el argumento del método send() estoy concatenando las etiquetas de 'CH0' y 'CH1' para poder separar las señales después con LabVIEW. (recuerden que '\n' es una secuencia de escape que indica nueva linea).

El diagrama a bloques del cliente para LabVIEW es bastante sencillo:
(Click para agrandar)
En el bloque TCP Open Conection colocamos la dirección IP* de la Raspberry y el puerto que establecimos para el socket en el programa de Python. Después de hacer la lectura utilizamos los bloques Match Patern para identificar las etiquetas de los canales y separar sus valores. Finalemnte se convierten las cadenas aisladas a numeros de tipo double para enviarlos al graficador. El sistema funcionando se ve así:


*Aquí estoy usando una ip local. En la práctica se requerirá utilizar una IP pública. Hay varias formas de hacerlo: #1, #2 , #3

Referencias:
Python Advanced Tutorial 6, DrapsTV
Raspberry Pi + Arduino + LabVIEW, Alfredo Cruz

domingo, 3 de abril de 2016

Modelo lineal de un motor DC

El modelo lineal de un motor eléctrico de DC consiste en 2 ecuaciones diferenciales acopladas: el modelo eléctrico y el modelo mecánico. Debido a que ambos modelos están relacionados podemos escribir un modelo general el cual nos permitirá obtener una función de transferencia. En este ejemplo encontraremos la función de transferencia que relacione voltaje (entrada) con posición angular (salida). Primero, observemos el diagrama del circuito equivalente de la armadura del motor y el diagrama de cuerpo libre de rotor:
Para obtener la ecuación diferencial para el modelo eléctrico consideramos la ley de voltajes de Kirchoff:
 La fuerza contra-electromotriz se genera al iniciar el movimiento del rotor debido a que el campo magnético fijo del estator induce un voltaje en el devanado de la armadura (este voltaje es negativo con respecto al voltaje de entrada). La fcem es proporcional a la velocidad angular del rotor por lo que la constante Kb puede determinarse experimentalmente graficando el voltaje en las terminales del motor contra la velocidad angular del rotor (se verá que la relación no es realmente lineal en un intervalo grande pero puede usarse sólo la región lineal como una aproximación para el modelo).

Para obtener el modelo mecánico consideramos la segunda ley de Newton para movimiento angular:

Dónde J y b son el momento de inercia del rotor y el coeficiente de amortiguamiento por fricción respectivamente. Vemos que el torque es proporcional a la corriente en el motor. Asumiendo que no hay perdidas electromagnéticas ni por calor:
Por lo que el valor para el factor de fcem encontrado experimentalmente puede ser usado como valor para Kt. Observamos ahora que las ecuaciones (1) y (2) están relacionadas por la función de corriente.Para facilitar la sustitución obtendremos primero la transformada de Laplace de ambas ecuaciones:
Despejando I(s) de (4) y sustituyendo en (3) obtendremos el modelo unificado para el motor DC. Teniendo ya una sola ecuación obtenemos la función de transferencia posición/voltaje:


Ejecutando el siguiente código en Matlab asignando algunos valores a las constantes obtenemos su respuesta al escalón unitario en lazo cerrado:

J = 3.2284E-6;
b = 3.5077E-6;
K = 0.0274;
R = 4;
L = 2.75E-6;
s = tf('s');
P_motor = K/(s*((J*s+b)*(L*s+R)+K^2))
sys_cl = feedback(P_motor,1)
step(sys_cl)




Referencias:
DC Motor Position: System Modeling, University of Michigan
DC Motor, MathWorks
DC Motor, How It Works? [video subtitulado en español], Learn Engineering