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.

domingo, 15 de noviembre de 2015

Temas de actualidad

Yo que soy un derechairo-liberal declarado tengo la sospecha que el atentado en París ha sido un evento de falsa bandera. De ser verdad podría ayudarnos a nosotros, ratas de ese laboratorio sociológico llamado Facebook, a decidirnos de que lado estar en esta batalla virtual. Pero ahora, ¿hasta que punto una bandera representa al Leviatán o representa al pueblo? Aún resolviendo eso queda otro problema, ¿a que pueblo se le debe mostrar empatía? A todos en el caso ideal, por supuesto. Pero esta solución genera otro problema; por más utopías que quieran sacarse de la manga, los seres humanos no podemos vivir sin fronteras (en todos los contextos posibles). Siempre he tenido como lema personal esta frase de Huygens, "El mundo es mi país, la ciencia mi religión", pero me he dado cuenta que en un contexto no político, esto no es verdad para mí. Existen muchos grupos sociales que respeto (siendo terriblemente sincero, otros no) más no puedo identificarme. Cargo en mi mente, como todos, fronteras invisibles definidas por reglas de pertenencia social. Y aquí una triste verdad. Aunque lo correcto sea mostrar empatía por la humanidad entera, la paradoja de la elección viene a arruinaros la paz. Elegir todo no es diferente a elegir nada.

miércoles, 11 de noviembre de 2015

Guanajuato

Hace unos días estuve una escuela-taller en el CIMAT, lo que fue el pretexto perfecto para conocer el instituto y la ciudad. Me gustó bastante y fue una lastima tener tan pocos días para recorrerla. Si bien el taller me dejó un poco frío me llevé buena idea de las áreas en las que trabaja el CIMAT y bastante información del posgrado. Durante la estancia me quedé en la casa del asesores de un amigo mío, que por cierto conocí en un congreso en Guadalajara el año pasado (aquí aquel breve episodio). Su asesor es una persona realmente agradable, un astrofísico brasileño de la UG. El primer día me recibió dando una clase de tango a la que no pude rechazar la invitación de unirme pues era algo que había querido hacer desde hace mucho. Conocí también al housemate de mi amigo, uno de sus mejores amigos que estudió junto con el en la BUAP que igual me cayo muy bien. En las tardes, después de que regresaba de las sesiones taller, bajábamos caminado al centro. Algunas partes de la ciudad me recordaron bastante a Venecia. El penúltimo día dimos un recorrido por los cerros que rodean a la ciudad, algo que suelen hacer. En algún momento tuvieron que irme a rescatar cuando me quedé aferrado como gato asustado a una pared de piedra sin poderme mover. Lo más gracioso fue que en aquel momento sentí que estaba en risco de 300 metros cuando en realidad no estaba tan alto. Ya en lo alto la vista era increíble. Debo decir que mi amigo lleva una vida de estudiante envidiable. Para mi fue un corta pero buena aventura.

jueves, 29 de octubre de 2015

Existencialismo

Si me preguntaran cuales han sido los 3 libros que marcaron mi vida, simplemente diría que sólo hay un ensayo que realmente lo ha hecho, y lo ha hecho más tarde de lo que hubiera querido. Me refiero a "El Existencialismo es un Humanismo" de Sartre. Pasé casi toda mi vida culpando al mundo de mis fallas, pasando por alto una idea tan fundamental; la existencia del hombre precede a su esencia y somos nosotros mismos quienes nos definimos. Quizá en fondo lo sabemos y nos aterra la idea de haber nacido en un mundo indiferente que existe y debe definirse desde cero. Que tarde lo he conocido señor Sartre.

martes, 27 de octubre de 2015

Convertir tiempo UTC a GMST (Tiempo Sideral Medio de Greenwich) en Python

 -*- coding: utf-8 -*-
"""
Created on Tue Oct 27 18:37:49 2015

@author: Rodolfo Escobar et al
"""

from datetime import datetime
from astropy.time import Time
from numpy import floor

# Reductor de rango 0-24h
def to24(x):
 while x >= 24:
  if x < 0:   
   x = x + 24
  if x > 24:
   x = x - 24
 return x
#


# Conversor de formato de hora
def tohms(x):
 time_hours = x
 time_minutes = time_hours * 60
 time_seconds = time_minutes * 60

 horas    = int(floor(time_hours))
 minutos  = int(floor(time_minutes % 60))
 segundos = int(floor(time_seconds % 60))
 return (horas,minutos,segundos)
#

#Programa principal
ut = Time(datetime.utcnow(), scale='utc')
JD= ut.jd # Fecha Juliana
D = JD - 2451545.0
GMST = 18.697374558 + 24.06570982441908*D
GMST = to24(GMST)
(horas,minutos,segundos) = tohms(GMST)
#

print "Hora Sideral:","{h}:{m}:{s}".format(h=horas, m=minutos, s=segundos)


Referencias:
Approximate Sidereal Time, USNO
Local Sidereal Time, Durham University
Astronomy With Your Personal Computer, Peter Duffett-Smith
Astropy.org

viernes, 23 de octubre de 2015

winvid en Matlab 2014/2015

Por alguna razón no viene instalada por defecto ninguna interfaz de video en las ultimas versiones de Matlab. Para resolver esto solo hay que escribir en la consola supportPackageInstaller e instalar el siguiente paquete:


Con esto podrás adquirir la imagen de tu webcam sin problemas después de reiniciar Matlab.

lunes, 28 de septiembre de 2015

Interrupciones externas en Raspberry Pi con Python

A partir de la versión 0.5.6 de la librería RPi.GPIO es posible utilizar rutinas de interrupción externas. Como los usuarios de microcontroladores experimentados deben saber, una rutina de interrupción permite ejecutar una o una serie de instrucciones sin perder tiempo en bucles de lectura. Es decir que si el procesador se encuentra ejecutando una tarea, éste pausará el proceso en curso y ejecutará la rutina de interrupción ante un flanco de subida o bajada ocurrido en un pin GPIO. Si no has comprado tu Raspberry recientemente y no has actualizado su software es recomendable hacerlo:

sudo apt-get update 

sudo apt-get upgrade #(Esto podría tardar hasta más de una hora)

El siguiente código es un ejemplo sencillo sobre el uso de interrupciones. Tenemos dos rutinas de interrupción CuentaA() y CuentaB() asociadas a GPIO23 y GPIO24 respectivamente. Cada una incrementan una variable de conteo y muestran el resultado en consola. El programa se detiene cuando el contador A es mayor o igual a 5. Para más detalles sobre interrupciones pueden revisar la documentación de la librería RPi.GPIO. Para evitar el el efecto de rebote se utilizan capacitores de entre 22 o 10 uF conectados en paralelo con las resistencias pulldown de 10K como se muestra en el circuito [aunque también puede utilizarse un debounce vía software agregando el argumento bouncetime=200 a la función add_event_detect()]:
import RPi.GPIO as GPIO
import os

contaA = 0
contaB = 0

#Setup
GPIO.setmode(GPIO.BCM)
GPIO.setup(23,GPIO.IN)
GPIO.setup(24,GPIO.IN)

#Callbacks
def CuentaA(channel):
    global contaA
    contaA += 1
    os.system("clear")
    print "Contador A: ", contaA
    print "Contador B: ", contaB

def CuentaB(channel):
    global contaB
    contaB += 1
    os.system("clear")
    print "Contador A: ", contaA
    print "Contador B: ", contaB

#Interrupciones
GPIO.add_event_detect(23, GPIO.RISING, callback = CuentaA)
GPIO.add_event_detect(24, GPIO.RISING, callback = CuentaB)

print "Contador A: ", contaA
print "Contador B: ", contaB

#Bucle principal
while(contaA < 5):
    pass

GPIO.cleanup()

jueves, 17 de septiembre de 2015

Filtro de Notch de 60 Hz

Debido a la restricción de los valores comerciales para los resistores y capacitores, este circuito tiene una frecuencia central de 58.94 Hz. Sin embargo, como se puede ver en una simulación, se tiene una atenuación de -34 dB (un factor de de 0.0004) para la componente de 60 Hz.


Para más detalles sobre este tipo de filtros pueden consultar éste artículo.

domingo, 13 de septiembre de 2015

Circuito para electrocardiografía (ECG/EKG)

Este circuito está basado en el que se puede encontrar en la hoja de datos del INA114. En el circuito original se asume que el rango de voltaje de la señal de entrada es de 1 mV ~ 5 mV más un offset de DC de 300 mV. Sin embargo, en la práctica estos valores resultan ser aún más pequeños (posiblemente por las condiciones poco controladas en un laboratorio de ingeniería universitario o el uso de caimanes y protoboard) así que es necesario incrementar la ganancia del amplificador de instrumentación. Para esta modificación se tiene una ganancia de 2501 y estos fueron los resultados:


En este video utilizan un circuito con ganancia de 100 con buenos resultados (el cual no funcionó en mi caso), lo que me hace sospechar que existen varios factores que afectan el voltaje de entrada. Posiblemente este circuito no funcionará en condiciones diferentes por lo que recomiendo experimentar con diferentes valores de ganancia.

Un segundo circuito que también funciona muy bien es una modificación del que pueden encontrar en The Biosignal How-To:
 

miércoles, 29 de julio de 2015

“Que nadie trate de deducir quién fui 
de todo lo que hice y todo lo que dije. 
Había un obstáculo que deformaba 
mis acciones y mi modo de vivir. 
Había un obstáculo que me detenía 
muchas veces cuando iba a hablar. 
Por medio de mis acciones más inadvertidas 
y mis escritos más velados, 
sólo por medio de estas cosas podré ser comprendido. 
Pero quizá no valga la pena dedicar 
tanto interés y tantos esfuerzos a descubrir quién soy. 
Más adelante –en una sociedad más perfecta- 
otro, hecho exactamente como yo, 
sin duda aparecerá y actuará con libertad.” 
 — Cosas escondidas, Constantino Cavafis.

lunes, 27 de julio de 2015

Morelia

A mediados de mayo recibí dos muy buenas noticias el mismo día. Mi universidad aprobó el financiamiento de mi proyecto de tesis (un radiotelescopio de 1m para la banda Ku) y la noticia de que había sido aceptado en una escuela de verano que había estado esperando con muchas ansias. Sin embargo había un pequeño detalle en esto último. Como había temido, iba a toparme durante dos semanas con mi ex nuevamente. Después de la mala experiencia de Monterrey había decido no permitir que volviera a ocurrir así que tuve que contactarla, después de meses sin hablarnos, para negociar una reconciliación amistosa. Esto ocurrió una semana antes de la escuela. Para mi sorpresa aquella tarde terminó muy bien. Aclaramos muchas cosas sobre el congreso de Monterrey, hablamos como buenos amigos. Aquel día también escuché lo que tanto había temido. Ya salía con alguien más. Debo decir que en aquel momento lo tomé muy bien, e incluso pasé el resto de la tarde sintiéndome orgulloso de mi madura respuesta, como si subir un nivel de la vida se tratara. Pero la realidad me golpeó con toda su fuerza al siguiente día al verlos besándose en mi propia facultad, pues por fatalidades de la vida el susodicho era compañero mío. Vislumbrado un infierno para las siguientes semanas, pero no fue así. Por mucho, fue una de las mejores experiencias que viví en este año. Académicamente me abrió todo un panorama que desconocía. Resolví muchas de mis eternas dudas sobre síntesis de radioimágen y conocí en persona a toda una leyenda de la radioastronomía en México, Stan Kurtz. Manejé remotamente una antena de 35 metros en España (Proyecto PARTNeR), procesé datos crudos de ALMA, datos de rayos X de Chandra y aprendí un poco sobre software del NRAO. No sentía una euforia tan intensa por el conocimiento desde hacía un buen tiempo. Stan, una persona muy amable a quién visité frecuentemente mostrándole mis proyectos y exponiendo mis dudas, acepto mi propuesta de trabajar con el como voluntario o becarlo el próximo año y apenas lo podía creer ¡Tanto tiempo que lo había creído inaccesible para mi! En cuanto mi ex... si hubiera escrito en esos días de junio estaría hablando enérgicamente sobre el "mito" de la imposibilidad de la amistad después del amor y demás. Nos pasamos hablando durante el largo viaje desde Puebla hasta Morelia en autobús. A lo largo de la estancia tuvimos algunas peleas, pero como en los buenos tiempos, aclarábamos las cosas escuchándonos. Hubo un día en el que no pude evitar recordar a un bloguero moreliano que contaba sus aventuras y desamores en aquella ciudad. Ahora era yo quien estaba en lagrimas en una noche en zócalo de Morelia. Fue  toda una sopresa llevarme tan bien con ella en esas dos semanas. Aunque nunca lo aceptó, notaba como se le llenaban los ojos de algunas lagrimas en los momentos que reíamos o cuando le contaba apasionadamente sobre mis platicas con Stan, como solíamos hacerlo alguna vez... Es una lastima que se alejó de mi tan pronto volvimos del viaje.

domingo, 19 de julio de 2015

Comunicación serial con Arduino y Scilab

 Recuerda instalar Serial Communication Toolbox primero con Module Manager.

Programa en Scilab

h=openserial("/dev/ttyACM0","9600,n,8,1"); // *
xpause(1e06); // retardo (aprox. 1 sec)
writeserial(h,char(1)) // envía 1(uint8)
xpause(1e06) // retardo (aprox. 1 sec)
writeserial(h,char(0)) // envía 0 (uint8)
closeserial(h);

* En Windows se debe sustituir "/dev/ttyACM0" por el número (sin comillas) de COM. Por ejemplo, para el COM5 quedaría: h=openserial(5,"9600,n,8,1");

Probamos la comunicación con un programa simple en Arduino dónde el led del puerto 13 encienda enviando un 1 y apagándolo con un 0.

Programa de prueba en Arduino

void setup() {
  Serial.begin(9600);
  pinMode(13, OUTPUT); 
  digitalWrite(13,LOW);  
}

void loop() {
    while (Serial.available() > 0) {
     int lectura = Serial.read();    
     if (lectura == 1){digitalWrite(13,HIGH);}
     else if(lectura == 0){digitalWrite(13,LOW);}   
      
    }    
} 

Como ejercicio, se puede escribir un programa dónde se envíe un valor PWM (0-255) desde Scilab hacía un led conectado a la tarjeta. 

miércoles, 15 de julio de 2015

Comunicación SPI con un DAC en VHDL

Para esta implementación se está utilizando un DAC LTC1661. Ya había escrito una entrada sobre comunicación SPI con este mismo DAC con Arduino. La ventaja didáctica de una implementación hardware, como en este caso, es que permite un mejor entendimiento de este protocolo de comunicación. Si más preámbulo, entremos al tema. Como quizá ya lo haya mencionado antes, las hojas de especificaciones lo son todo. Antes de siquiera pensar en escribir nuestro programa debemos entender que queremos hacer. Asumiré que ya se tiene cierto conocimiento del protocolo y me centraré en lo necesario para mantener breve la descripción del programa.

De la hoja de especificaciones, tenemos el siguiente diagrama de tiempos:
 Los periodos que nos interesan son t5, que es el periodo mínimo entre datos ( el tiempo mínimo que el FPGA debe esperar para enviar el siguiente dato de manera que el DAC pueda reaccionar)  y t3 y t4 que nos indican que SCK debe tener un duty cycle de 50% con un periodo T = t3+t4. En la página 4 de la hoja podemos encontrar que los tiempos mínimos requeridos son: t5 = 100 ns y una frecuencia para SCK de 10 Mhz. Los periodos t1 y t2 nos indican una especificación igual de importante. Indican que el flanco de subida de SCK debe ocurrir justo a la mitad del bit enviado. Entendiendo estas restricciones, podemos ahora comenzar a escribir nuestro programa. Nuestra entidad deberá tener 2 entradas: clk y un reset, y 3 salidas: SCK, CD y datos:

La implementación estará formada por los siguientes bloques elementales:
  • Divisor de frecuencia (50 a 10 MHz)
  • Divisor de frecuencia para reloj de CS
  • Memoria ROM
  • Registro PISO (Parallel Input Serial Output)
  • Generador de reloj SCK
 Sólo describiré los bloques más relevantes:

Memoria ROM

Los datos almacenados por esta memoria deben tener la estructura especificada por el fabricante (pag. 8), que en este caso es la siguiente:

 A3A2A1A0D9D8D7D6D5D4D3D2D1D0X1X0
Control Datos Don't Care

Para este ejemplo debemos enviar la mitad del voltaje de referencia al canal A y actualizar su salida, por lo que los valores para el código de control deben ser "1001" (ver tabla de la pag. 9). Se crean entonces las siguientes señales que serán usadas en la memoria ROM:

type rom is array(0 to 2) of std_logic_vector(15 downto 0);
constant control: std_logic_vector(3 downto 0):="1001";    
constant dc: std_logic_vector(1 downto 0):= "00"; -- Bits don't care constant myrom: rom
    :=( control&"0111111111"&dc,
        control&"0111111111"&dc,
        control&"0111111111"&dc);


La razón de implementar una memoria ROM como un arreglo de varios elementos para un ejemplo tan simple es para que pueda ser modificada fácilmente para  hacer cosas más interesantes como enviar una función triangular o cosenoidal (en cuyo caso habrá que incrementar el tamaño del array a entre 50 o 100 elementos). Otro detalle que debo mencionar es que, como se verá en el programa completo, el incremento de la dirección se hace con el reloj CS. Esto se hace a modo de tener el bit-rate máximo. Se puede agregar un nuevo divisor a una frecuencia más baja para cambiar el dato si así lo desean.

Registro PISO

Este registro serializará los datos recibidos de la salida de la memoria ROM y los enviará de forma síncrona en cada flanco de bajada* del reloj secundario de 10 MHz. La serialización únicamente se efectuará después de recibir un flanco de subida del reloj CS.

*Esto de modo que se cumpla la especificación de los tiempos t1 y t2, así SCK deberá estar en fase con el reloj secundario.

Generador de reloj SCK

Este reloj deberá tener una frecuencia de 10 MHz pero requiere cumplir  una condición: debe activarse durante el siguiente flanco de subida del reloj secundario (10 MHz) después del flanco de bajada del pulso CS. Esto se logra mediante una implementación if síncrona normal. Sin embargo, dado que no es posible implementar físicamente una instrucción en flanco de bajada inmediatamente después de un flanco de subida ("bad synchronous description") se recuré a un truco combinacional para generar el reloj SCK que puede verse en el programa completo.

Teniendo el programa listo, así luce en simulación con ISim:

(Click para agrandar)

La implantación física se realizó en una tarjeta de desarollo Basys2 y funcionó perfectamente. Desafortunadamente, como había comentado en alguna entrada, sufrí un robo durante un congreso en Monterrey, por lo que no habrá fotos esta vez.

Programa completo

LIBRARY IEEE;
USE IEEE.STD_LOGIC_1164.ALL;
USE IEEE.STD_LOGIC_ARITH.ALL;
USE IEEE.STD_LOGIC_UNSIGNED.ALL;

entity DAC_SPI is

    port(
        clk         : in    std_logic; -- 50 MHz          
        reset       : in    std_logic;
        CS          : out   std_logic;
        SCK         : out   std_logic-- 10 MHz
        datos       : out   std_logic); -- Salida serial
   
end entity;

architecture rtl of DAC_SPI is
-- Señales
    -- Divisor

    signal count: integer :=1;
    signal clk_10: std_logic:='0';
   
    -- ROM
    signal address: integer range 0 to 2;
    type rom is array(0 to 2) of std_logic_vector(15 downto 0);
    constant control: std_logic_vector(3 downto 0):="1001"; -- Cargar y actualizar canal A
    constant dc: std_logic_vector(1 downto 0):= "00"; -- Bits don't care
    constant myrom: rom
    :=( control&"0111111111"&dc,
        control&"0111111111"&dc,
         control&"0111111111"&dc);
    signal rom_out: std_logic_vector(15 downto 0):=(others=>'0');
    signal clk_cs:  std_logic:='0';
   
    -- PISO
    signal temp: std_logic_vector(15 downto 0):= (others =>'0');
    signal t:    std_logic;
   
    -- Otros
    signal a: integer range 0 to 18 :=1;
    signal sck_t: std_logic;
--   
begin

-- Divisor de frecuencia (10 MHz)
process(clk) begin
   if(clk'event and clk='1') then
      count <=count+1;
      if(count = 5) then -- dónde count = frec de reloj / frec deseada
         clk_10 <= not clk_10;
         count <=1;
end if;
end if;
end process;
--

-- PISO

   process (clk_10,rom_out,clk_cs)begin
     if (CLK_10'event and CLK_10='0') then
       if (clk_cs='1') then
             temp(15 downto 0) <= rom_out(15 downto 0);
       else
             t <= temp(15);
             temp(15 downto 1) <= temp(14 downto 0);
             temp(0) <= '0';
       end if;
     end if;   
    end process;

    datos <= t;
--
-- SCK

    process(clk_10,clk_cs) begin
        if (clk_10'event and clk_10 = '1') then
                if(clk_cs='0')then
                    SCK_t <= '1';
                else
                    SCK_t <= '0';
                end if;
        end if;       
    end process;

  SCK <= (sck_t) and (clk_10) and not(clk_cs);   
--               

-- Divisor de frecuencia para CLK_CS

process(clk_10) begin
   if(clk_10'event and clk_10='1') then
      a <=a+1;
      if(a = 17) then
         clk_cs <= '1';
         a <=0;
        else
            clk_cs <= '0';
        end if;
    end if;
end process;
--

-- ROM

process(clk_cs,reset) begin
    if reset = '1' then
        address <= 0;
    elsif( clk_cs'event and clk_cs = '1' ) then
        if address < 2 then
            address <= address + 1;
        else
            address <= 0;
        end if;
    rom_out <= myrom(address);
    end if;
end process;

-- CS
CS <= clk_cs;
--
end rtl;

martes, 14 de julio de 2015

New Horizons finalmente en Plutón

Imagen de Plutón a 766000 kilómetros de distancia tomada por la New Horizons 
16 horas antes del encuentro (NASA/APL-JHU/SwRI).
  
9 años había esperado para ver esta fotografía y ha valido totalmente la espera. Un símbolo de lo lejos que puede llegar la inventiva humana. 

 "A la espera de los datos del encuentro, hoy ya vivimos en un mundo que sabe qué aspecto tiene Plutón. A partir de ahora, cuando queramos hablar de Plutón o de un objeto del cinturón de Kuiper ya no usaremos borrosas imágenes del telescopio Hubble o pinturas de artistas, porque al fin tenemos imágenes reales obtenidas por un artefacto humano que hemos mandado a los límites del sistema solar."

lunes, 13 de julio de 2015

Simulando circuitos con ngspice

Ngspice es un simulador de circuitos gratuito y de código abierto desarrollado por el CAD Group de la Universidad de California, Berkeley. Si ya están familiarizados con otros programas basados en SPICE no debe ser mayor problema aprender a usarlo. Por otro lado, si están acostumbrado a las comodidades de programas como Multisim o Proteus hay algunas cosas que deben saber para empezar. La forma de trabajar con ngspice es la siguiente:
 Es necesario crear un directorio en donde estén guardados los archivos .mod (SPICE models) que contienen los modelos de simulación de los dispositivos electrónicos que van a utilizarse. Estos archivos no están incluidos en la instalación por defecto de ngspice, pero son proveídos por los mismos fabricantes de componentes. Estos modelos deben ser enlazados a los componentes desde gschem.

Instalación en Debian/Ubuntu/Linux Mint:

$ sudo apt-get install ngspice

Para instalar gschem y gnetlist instalamos el paquede gEDA:

$ apt-get update && apt-get install geda pcb gerbv

Abrimos gschem:

$ gschem

, buscamos "añadir componente" en la barra de herramientas (su ícono es una compuerta and) y usamos los siguientes componentes:

SPICE simulation elements > vsin-1.sym
Diodes (generic) > diode-1.sym
Basic devices > resistor.sym
Power rails > gnd-1.sym

La interfaz de gschem es bastante intuitiva. Para girar 90° un componente sólo hay que seleccionarlo y presionar las teclas E y R (una después de la otra). Una vez terminado el circuito se deben configurar los componentes para poder generar un netlist que pueda ser simulado correctamente por ngspice. Dando doble click en el diodo se abrirá una ventana dónde debemos añadir una propiedad value cuyo valor será 1N4001 y damos click en "añadir":
Hacemos lo mismo con la fuente de AC. En este caso, value ya aparecerá en la tabla por lo que solo debemos modificar su valor directamente escribiendo "sin 0 5 60Hz". Para agregar un valor a la resistencia, le damos click derecho a su símbolo y buscamos "añadir propiedad", se abrirá una ventana dónde en "nombre"  escribimos "value" y en "valor" escribimos su resistencia en Ohms (100 para este ejemplo), en la última opción seleccionamos "Mostrar solo valor". Después de esto ahora sólo hace falta enlazar el archivo de modelo de SPICE al símbolo del diodo. Para esto regresamos a "agregar componente" y buscamos SPICE simulation elements > spice-model-1.sym. Para este modelo en particular las propiedades deben ser las siguientes:

Observen que en file he escrito la ruta al archivo de modelo para el diodo 1N4007 (pueden descargar un modelo de ese diodo aquí, cambiando la extensión de .txt a .mod) . Guardamos el esquema con el nombre diodo.sch. El esquema terminado debe verse algo así:

Para generar el netlist escribimos:

$ gnetlist -g spice-sdb -o diodo.net diodo.sch

[NOTA: consultar sintaxis y opciones de gnetlist]

Ejecutamos con ngspice:

$ ngspice diodo.net

Ya dentro de la consola de ngspice escribimos:

tran 100us 100ms #dónde se indica 100us de step simulando de 0 a 100ms
plot v(2)

Se abrirá una ventana con la gráfica del voltaje en la resistencia:


NOTA: Deben tomar siempre en cuenta el nombre del modelo que está declarado dentro del archivo .mod. El modelo del 1N4007 que estoy usando es diferente al que pueden descargar en diodes.com. Si usan ese archivo, deben asigar valores para la propiedad value de DI_1N4007 en vez de 1N4007 (o simplemente cambiar el nombre en el script del modelo). Este es el contenido del archivo de diodes.com, señalo el lugar donde va el nombre del modelo:

*SRC=1N4007;DI_1N4007;Diodes;Si;  1.00kV  1.00A  3.00us   Diodes, Inc. diode
.MODEL DI_1N4007 D  ( IS=76.9p RS=42.0m BV=1.00k IBV=5.00u
+ CJO=26.5p  M=0.333 N=1.45 TT=4.32u )

viernes, 19 de junio de 2015

With our eyes so full of evening

domingo, 7 de junio de 2015

Monterrey

En marzo de este año presenté un póster sobre un trabajo de radioastronomía amateur en un congreso en Monterrey. El trabajo consistía en una modernización del proyecto RadioJOVE de la NASA utilizando sintonizadores de televisión HD (SDR-RTL) de bajo costo. En aquellos días estaba demasiado sentimental como para hacer una reseña más o menos objetiva de aquel viaje. Alguna vez un amigo, estudiante de física, me advirtió del error que es andar con una mujer de tu misma área profesional, un consejo que me llegó muy tarde (y que meses después mi amigo tuvo que tragarse). Yo no diría que es un error pero si una apuesta arriesgada. Si las cosas funcionan por años, debe ser una experiencia muy bonita. Congresos, escuelas de verano, escuelas de invierno, estancias de investigación, viajar por el país y por el mundo juntos como pareja haciendo lo que los les gusta. Uno de esos romances dignos de una pantalla grande. Pero no, yo no tuve esa suerte. Cuando uno pierde esa apuesta se paga con una experiencia que puede llegar a ser bastante incomoda. Congresos, escuelas de verano, escuelas de invierno, estancias de investigación, viajar por el país y por el mundo topándote con exnovia, quizá por el resto de tu vida profesional. Por supuesto, si uno queda en buenos términos de amistad puede ser divertido, pero aquí es donde vuelvo a perder también.  Y como si el destino gozara de ver en desgracia a los mortales, mi cámara fue robada en el hotel en una situación tan extraña que parecía un relato de Sherlock Holmes. C'est la vie. No fue el viaje perfecto pero a pesar de todo lo disfruté mucho, y creo que eso es algo que puedo decir hasta ahora. Iba acompañado de buenos amigos que siempre supieron como animarme, aprendí mucho en muchos aspectos y, como siempre ocurre en estos eventos, conocí gente interesante y muy agradable. El drama de todos los días de un estudiante de ciencias.

sábado, 6 de junio de 2015

Escribir y compilar un programa en C con gcc desde terminal

Abrimos una consola, asegurándonos de que estamos en nuestra carpeta personal (puedes hacerlo con el comando pwd) y escribimos:
$ mkdir Programas  # Creamos una nueva carpeta
$ cd Programas       # Nos movemos a la carpeta Programas
~/Programas$ nano ejemplo.c # Creamos nuevo archivo .c
~/Programas$ gcc -o ejemplo ejemplo.c # Compilamos y creamos código objeto con el nombre ejemplo.o
~/Programas$ ./ejemplo # Ejecutamos código objeto

viernes, 5 de junio de 2015

Meetup

Esta semana hice algo de lo que tenía mucha curiosidad: ir a un meetup. Este en particular fue muy interesante ya que la temática era hacer un intercambio cultural entre extranjeros radicando en Puebla y mexicanos. Compartí una buena tarde con gente de Turquía, India (un estudiante de doctorado de astrofísica en INAOE) y EU (Ian, el fundador del grupo Foreing Exapts, Mexicanos and Poblanos). Meetup no es una plataforma muy conocida en México, sumémosle que la seguridad aquí no es nuestro fuerte, y es difícil encontrar eventos en la provincia. Espero que eso pronto pueda cambiar. Es justo el tipo de sitios donde el internet se usa para romper tu burbuja de confort y vivir algo nuevo en la vida real.

viernes, 10 de abril de 2015

sábado, 21 de marzo de 2015

02:25 am

Me arrepiento de perder la costumbre, si alguna vez la tuve, de escribir con regularidad. Este blog sería como un registro al que podría acudir para auxiliar a mi memoria o analizar mi vida de forma más ordenada. Una libreta de apuntes técnicos de mi vida. Pero siendo honesto conmigo mismo, esa no era la verdadera razón por la que yo escribía. Era y es una necesidad de dejar una huella en este mar de la existencia. Puedo presumir que a mis 24 años he llevado una vida de novela literaria. ¿No es esto un ya buen motivo para escribir sobre ella? Quizá temo palidecer los detalles sobre una mala redacción. But who knows? nunca es tarde para recuperar o construir un hábito.

jueves, 19 de marzo de 2015

Gráfica de un Patrón de Radiación 3D con Python

Para este programa se están utilizando datos obtenidos de una simulación de un dipolo de media onda en HFSS en un archivo csv.

# -*- coding: utf-8 -*-
"""
@author: Rodolfo Escobar
"""


import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

data = np.loadtxt('dipolo_16.csv',delimiter=',',skiprows=1)

theta = data[:,0]
phi = data[:,1]
r = data[:,2]

a = 0.0174532925
phi_rad = a*phi;
theta_rad = a*theta;

#Conversión de coordenadas
x = r * np.sin(phi_rad) * np.cos(theta_rad)
y = r * np.sin(phi_rad) * np.sin(theta_rad)
z = r * np.cos(phi_rad);

x = np.reshape(x,(73, 37))
y = np.reshape(y,(73, 37))
z = np.reshape(z,(73,37))
R = np.reshape(r,(73,37))

N = R/R.max()

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(x, y, z, rstride=1, cstride=2,facecolors=cm.coolwarm(N),cmap = 'coolwarm',linewidth=0.5 )
surf.set_edgecolor('k')

# Equivalente a la instrucción axis equal en Matlab
max_range = np.array([x.max()-x.min(), y.max()-y.min(), z.max()-z.min()]).max() / 2.0

mean_x = x.mean()
mean_y = y.mean()
mean_z = z.mean()

ax.set_xlim(mean_x - max_range, mean_x + max_range)
ax.set_ylim(mean_y - max_range, mean_y+ max_range)
ax.set_zlim(mean_z - max_range, mean_z + max_range)
#

plt.title(u'Patrón de Radiación')
m = cm.ScalarMappable(cmap=cm.coolwarm)
m.set_array(R)
cbar = fig.colorbar(m,shrink=0.8, aspect=9)
cbar.set_label('dB')
plt.rc('ytick', labelsize=7)
plt.show()




lunes, 16 de marzo de 2015

Graficar una esfera con coordenadas esféricas en Matlab

Matlab utiliza la convención alta-azimutal para coordenadas esféricas por lo que se tiene que considerar a la hora de traducir funciones expresadas con otras convenciones.

%Esfera de radio unitario en coordenadas esféricas
theta = linspace(0,2*pi,50);
phi = linspace(-pi,pi,50);
r = ones(1,50);

%Construcción de malla
[phi,theta] = meshgrid(phi,theta);
r = meshgrid(r);

%La conversión se realiza después de construir la malla
[X,Y,Z] = sph2cart(theta,phi,r);

mesh(X,Y,Z), axis equal