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,inf,'int16','ieee-be');
img = reshape(img,dim)';

%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 Mealy 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 rotor debido a que su campo magnético fijo 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 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

martes, 22 de marzo de 2016

Configuración inicial en Debian

[Esta entrada es sólo un apunte personal para recordarme lo que debo de hacer depués de instalar Debian para configurarlo a mi gusto.]

1. Instalación y configuración de WiFi:

https://wiki.debian.org/es/WiFi

2. Configuración de GNOME:
Escribir "Retoques" o "Tweak" en el buscador de aplicaciones. Abrimos la aplicación para configurar GNOME. En la pestaña de "Ventanas" podemos agregar los botones de maximizar y minimizar

3. Instalación de códecs multimedia
apt-get install libavcodec-extra 
4. Habilitación de sudo

En  Debian no vienen configurados sudoers por defecto pero pueden encontrar un tutorial para instalarlos aquí.

5. Instalación fuentes de Windows

apt-get install ttf-mscorefonts-installer

6. Instalación de Cairo Dock
Para instalar el dock solo tenemos que abrir Synaptic y escribir "Cairo Dock" en el buscador. Instalamos doble click en el paquete y seguimos el procedimiento de instalación dando click en "Aplicar". Una vez instalado, presionamos Alt+F2 y escribimos "cairo dock" para iniciarlo. Después abrimos la aplicación de retoques como (el mismo del punto 2) y abrimos la pestaña "Aplicaciones al inicio". Damos click en el signo "+" y damos click en Cairo Dock. De esta manera el dock se iniciara automáticamente al iniciar sesión.


martes, 1 de marzo de 2016

Programa de ejemplo en C para MicroBlaze

 
(Tarjeta de desarrollo Symbhia)

#include <xparameters.h>
#include <xgpio.h>

XGpio LEDS;
XGpio SW;

u8 sws; // Variable tipo uint8
int i;

int main(){
//Inicialización
  XGpio_Initialize(&LEDS,XPAR_GPIO_1_DEVICE_ID);
  XGpio_Initialize(&SW,XPAR_GPIO_0_DEVICE_ID);

while(1){
//Lectura de switches
  sws = XGpio_DiscreteRead(&SW,1); //(puntero de instancia,canal)

//Escritura de la variable sws en LEDS
  XGpio_DiscreteWrite(&LEDS,1,sws); // (puntero de instancia,canal,dato)
   }

return 0;

}

sábado, 20 de febrero de 2016

OpenCV en Python

Instalación de Anaconda en Linux:

Descargar instalador aquí y ejecutar con:

$ bash nombredelinstallador.sh

Instalación de OpenCV usando conda (en la terminal):

$ conda install opencv 

Abrir Spyder:

$ spyder

Ejemplo. Cargar una imagen y mostrarla con matplotlib:

import cv2 
from matplotlib import pyplot as plt 

img = cv2.imread('eco.jpg',0) #Cargar imagen 
plt.imshow(img,cmap = 'gray') # Mostrar imagen 
plt.xticks([]), plt.yticks([]) #Ocultar numeración en ejes 
plt.title('Umberto Eco (1932-2016)'


jueves, 14 de enero de 2016

8 de Enero

Apenas corre la primera semana de este año y ya da mucho de qué hablar. Esta noche, después de un curioso suceso, me he propuesto  volver a escribir frecuentemente pues me he terminado de convencer finalmente de que olvidamos lentamente la vida de una forma silenciosa. Pequeños detalles se pierden para siempre. Un beso en al lluvia, una carta hecha pedazos o un mensaje de texto hiriente. Olvidamos, convenientemente, que también hemos sido unos hijos de puta. Pero hablemos un poco del año pasado. Inicié 2014 con expectativas en el suelo, quebrado por dentro, emocionalmente abatido. No sabía que esperar de 2015 y fuera lo que fuera sabía que no estaba listo (y tenía razón). Doce meses después esperaba la llegada de otro año sentado en un viejo sillón en la terraza de la casa mi abuela con una extraña tranquilidad. Sabía que había sido un año memorable por todo lo bueno y lo malo. Académicamente ha sido quizá mi mejor año hasta ahora, pero emocionalmente fue un año que me dejó algunas cuantas heridas de guerra. Lastimé y fui lastimado, conocí verdaderamente a la filosofía y definí mi ética personal. Veo con todo esto a 2016 como el incierto año que es pero esta vez sintiéndome listo y armado para el  mal tiempo. Si el hombre sólo puede definirse en la medida que se realiza, como decía Sartre, tengo el ánimo de ver cada día de este año como una página en blanco que no quiero dejar que pase en vano. O lo que quizá sea aún peor: olvidarlas después haberlas escrito como tantas hojas que el tiempo arranca de mi memoria.

viernes, 4 de diciembre de 2015

Implementación de un controlador PID discreto en un PIC18F2550

En esta entrada se utilizará el método desarrollado en la entrada anterior, por lo que saltarán algunos pasos durante el proceso. En esa entrada se implementó un controlador PID para una planta descrita por un doble integrador. En esta ocasión haremos una simulación de un controlador PID implementado en un PIC para controlar una planta construida con amplificadores operacionales. Si bien es posible construir un doble integrador físicamente, resulta impráctico debido a la saturación provocada por la acumulación continua. Para resolver esto, esta planta estará caracterizada por dos integradores con reset resistivo o integradores de ganancia controlada (que una forma elegante de llamar a un filtro pasa-bajas activo de un solo polo) [1].  Eligiendo los valores de resistencia y capacitancia de tal manera que cada integrador tenga ganancia unitaria, la función de transferencia de los 2 integradores nos queda:
 Para controlar esta planta utilizaremos el siguiente controlador discreto a partir de la transformación backward Euler:
Dónde:

Obteniendo transformada inversa, la salida de control descrita como ecuación de diferencias nos queda:
Hacemos la siguiente sustitución de coeficientes:

Estos coeficientes pueden encontrarse manualmente siguiendo el procedimiento descrito aquí, o ejecutando éste programa de Matlab: [PID_Proteus.m]. Teniendo estos valores podemos realizar una verificación de la ecuación en Simulink:

Algo muy importante que se debe notar es que en esta simulación se está considerando la limitación física de la saturación de la salida del controlador real. Sin esta limitación la respuesta del controlador sería más rápida y con un sobreimpulso mucho menor. Sabiendo ahora que nuestra ecuación de diferencias es correcta podemos proceder a escribir el programa en el PIC en lenguaje C con el compilador MikroC. El programa completo se adjuntará más adelante. La función de control principal queda descrita como:

void pid(void) // Función para el controlador
{
   e2=e1; e1=e0; u2=u1; u1=u0; // update de variables
   y = ADC_Read(0);   // Lectura analógica por AN0
   yd = (double) y;  // Conversión a double
   vin = (yd/1023*5);
   y = ADC_Read(1);   // Lectura analógica por AN1
   yd = (double) y;  // Conversión a double
   reff = (yd/1023*5);
   e0 = reff - vin;  // Error
   u0 = ke0*e0 + ke1*e1 + ke2*e2 - ku1*u1 - ku2*u2; // Ecuación de diferencias
   salida = (int) floor(((u0+5)/10)*1023);
   if (salida > 1023) salida = 1023;  // Saturación de DAC
   if (salida <= 0) salida = 0;
   DAC_Output(salida); // Enviar dato
}


Es muy importante que ésta función sea ejecutada dentro de una rutina de interrupción de manera que se tenga control del tiempo de muestreo Ts en exactamente 0.01 segundos para que el controlador funcione correctamente. Para la salida analógica del controlador se esta utilizando un DAC LTC1660 [2][3]. Compilando el programa y simulando en Proteus:


 Nota: El relay en la simulación lo único que hace es evitar que la salida bipolar del DAC envíe un voltaje de saturación negativo al inicio [ver este programa].

lunes, 16 de noviembre de 2015

Implementación de un controlador PID discreto

Supongamos que queremos controlar el mismo modelo de una masa móvil sin fricción que habíamos tratado en esta entrada cuya función de transferencia es la siguiente:

Dados las ganancias del controlador PID que se habían considerado (recomiendo utilizar el app PID Tuning de Matlab si no tienes estos valores para la planta de tu interés) :

La función de transferencia en tiempo continuo de nuestro controlador nos queda:

En la práctica es recomendable utilizar un LPF en la etapa derivativa a modo de reducir el ruido. En este caso la FT es:

Dónde N representa la frecuencia de corte del filtro que suele ser de 100 radianes/s. Teniendo la función de transferencia de nuestro controlador, el proceso de implementación de este controlador en un sistema digital (PIC, Arduino, FPGA, etc) es el siguiente:

La ecuación de diferencias es la que describirá el algoritmo que deberemos implementar en algún lenguaje de programación o un lenguaje de descripción de hardware. Para realizar la transformación entre la FT en dominio de Laplace al dominio de Z debemos considerar las siguientes aproximaciones numéricas para la derivada:


De aquí es fácil notar que podemos representar a la variable s en terminos de z (dónde Ts es el periodo de muestreo) de las siguientes maneras:


Estas 2 transformaciones de la variable s tienen propiedades distintas. Mientras la transformación backward Euler es estable siempre que la FT en tiempo continuo sea estable, la transformación forward Euler no lo es en todos los casos. (Para más detalles consultar la sección 19.2.3.1 de "Pasive, Active and Digital Filters", Wai-Kai Chen [1]). Para este ejemplo utilizaremos la transformación forward Euler, por lo que nuestro controlador queda representado en tiempo discreto por la TF:
Para este ejemplo se sabe de antemano que este controlador es estable pero si se arriesga a utilizar foward Euler se recomienda verificar estabilidad (puede usarse sisotool de Matlab).  Verificamos en Simulink que nuestro controlador funcione correctamente:


Ahora podemos proceder a la conversión de formatos de la FT. Para eso utilizamos el este programa de Matlab: PIDdiscreto.mat. Usando los datos que nos muestra el programa la FT nos queda:


Reordenando (click para agrandar):


Calculando transformada inversa:


Finalmente, verificamos que la ecuación de diferencias sea correcta en Simulink:


En una próxima entrada se implementará este controlador en un código en C para PIC y Arduino.