Generación de una señal coseno discreta

Podemos generar una función senoidal o cosenoidal de forma discreta mediante una ecuación de diferencias que es obtenida a partir de la transformada Z de una función senoidal o cosenoidal según sea el caso. Esta ecuación de diferencias es fácilmente implementable en un microcontrolador o FPGA. En este ejemplo crearemos un generador de señal cosenoidal. Partimos de la transformada Z de una función coseno discreto considerándola como la respuesta al impulso de la función de transferencia de un sistema Y(z)/X(z):

Dónde W0 es la frecuencia angular en tiempo continuo y Ts es el tiempo de muestreo [ver mi entrada sobre el concepto de frecuencia en tiempo discreto]. Multiplicamos ahora por z^-2/z^-2 para obtener un modelo físicamente sintetizable (debido a que los valores positivos en los exponentes de z implican muestras desde el futuro):

Multiplicamos ambos lados por el denominador de la función de transferencia y despejamos para Y(z):

Aplicamos transformada Z inversa para obtener finalmente nuestra ecuación de diferencias:

Para simular el generador utilizaremos Xcos de Scilab que es una alternativa de código abierto de Simulink [si no están familiarizados con esta herramienta aquí hay un video curso en español]. Antes de colocar los bloques debemos ir a menú superior en Simulation>Set Context y definimos las siguientes variables:

Ts = 0.1
W0 = 2*%pi/10
k0 = cos(W0*Ts)
k1 = 2*k0
 
El diagrama en Xcos es el siguiente:

Recordemos que hemos supuesto desde el principio que x[n] = delta[n]. La función Delta de Kronecker (impulso discreto) es un super bloque formado por el siguiente subsistema:

Dónde el retardo 1/z está inicializado a 1. Obsérvese que en el diagrama principal estamos enviando los resultados de la simulación al workspace. Personalmente, las gráficas nativas Scilab me parecen feas. Para hacer una grafica en Python exportaremos los datos de la simulación a un archivo CSV con el siguiente comando en la consola de Scilab después de ejecutar la simulación:

csvWrite([A.time,A.values],'cos_data.csv',',')

Teniendo el archivo podemos generar una gráfica en Python con el siguiente código:

import numpy as np
import matplotlib.pyplot as plt

A = np.loadtxt('cos_data.csv',delimiter = ',')

plt.stem(A[:,0],A[:,1])
plt.ylim([-1.5,1.5])
plt.grid(True)
plt.xlabel('Ts*n')
plt.ylabel('y[n]')
plt.show()

  Se deja como ejercicio construir un generador para función seno.

Segmentación de imágenes por umbral con scikit-image

Segmentar una imagen a través de umbrales es probablemente la forma más sencilla de identificar regiones de interés en una imagen. Esto hace que este método sea vulnerable a ruido, sombras y otros factores más a la hora de utilizarlo en la práctica. Sin embargo es una forma sencilla de introducir a las funciones de segmentación del módulo scikit-image. Primero, para este ejemplo tomaremos la siguiente pintura de Kazimir Malevich:

 Separaremos la imagen en secciones para las cuales obtendremos un array con los índices de cada pixel para cada una de ellas. Para ello cargamos primero la imagen:

#Cargar imagen
im = skimage.io.imread("malevich-map-style.png")
imGS = skimage.color.rgb2gray(im) #Escala de grises

 y creamos una imagen binarizada a partir de un umbral en la imagen en escala de grises:

#Umbral
thresh = 0.7
bw = binary_erosion(closing(imGS < thresh, np.ones((3,3))))

La variable bw es la imagen binarizada corregida y es obtenida aplicando la función closing seguida de binary_erosion (ambas contenidas en skimage.morphology). Closing() realiza una operación de dilatación seguida de una erosión, lo que permite cerrar "agujeros" oscuros en una imagen. El primer argumento que toma en este ejemplo es la imagen binarizada resultante de aplicar el umbral imGS < 0.7 y el segundo es una matriz que determina el tamaño de la dilatación/erosión. Binary_erosion es solo una operación de erosión que elimina las falsas regiones debidas a las zonas brillantes del fondo. Teniendo la imagen binarizada limpia, obtenermos los segmentos:

#Segmentación
labeled_im = skimage.measure.label(bw)
lab_img = label2rgb(labeled_im, image=imGS,bg_label=0)

dónde labeled_im es una nueva imagen obtenida con la función label()en dónde cada pixel de la imagen en escala de grises original esta clasificado con su número de región . Si esto no es muy claro, pueden visualizar este arreglo de datos con el explorador de variables de Spyder dando doble click en el nombre de la variable:

La imagen lab_img es una visualización gráfica de las etiquetas de región sobre la imagen original que al mostrarla con Matplotlib y utilizando la función plt.text() junto con los datos contenidos en labeled_im, el resultado final es el siguiente:

 Código completo:

# -*- coding: utf-8 -*-
"""
Created on Tue Jun 27 13:27:53 2017

@author: rodolfo
"""
import numpy as np
from matplotlib import pyplot as plt
import skimage
from skimage.morphology import closing, binary_erosion
from skimage.color import label2rgb

#Cargar imagen
im = skimage.io.imread("malevich-map-style.png")
imR = im[:,:,0]
imG = im[:,:,1]
imB = im[:,:,2]
imGS = skimage.color.rgb2gray(im)

#Umbral
thresh = 0.7
bw = binary_erosion(closing(imGS < thresh, np.ones((3,3))))

#Segmentación
labeled_im = skimage.measure.label(bw)
lab_img = label2rgb(labeled_im, image=imGS,bg_label=0)

#Array de arrays de indices por segmento
reg_index = []
for i in range(0, np.max(labeled_im)):
    u = np.where(labeled_im == i)
    reg_index.append(zip(u[0],u[1]))

#Plots
R = []
cr = []
for i in range(0,np.max(labeled_im)):
    R.append(np.where(labeled_im==i+1))
    cr.append((np.mean(R[i][0]),np.mean(R[i][1])))

plt.figure(1)
plt.imshow(im)

plt.figure(2)
plt.imshow(lab_img)
for i in range(0,np.max(labeled_im)):
    plt.text(cr[i][1],cr[i][0],str(i+1),color='w')

Muestras I/Q: Un ejemplo con GNU Radio

En la entrada anterior resumía todo lo relevante sobre las señales analíticas (I/Q). Esta vez haremos un ejemplo que nos ayudará a visualizarlas mejor y hacernos de una intuición sobre su naturaleza. El ejemplo consiste en obtener las componentes I/Q de la siguiente señal:

Al ser una señal real, esperamos que tenga un espectro simétrico tal como podemos observar en el sink de frecuencia:

Para obtener las componentes I/Q de esta señal requerimos implementar un modulador en cuadratura con un filtro pasabajas para cada componente y posteriormente usar un bloque de conversión Float to Complex. No hará falta agregar simular un ADC ya que la señal generada ya es discreta. Si aún no están familiarizados con GNU Radio tengo un entrada introductoria [Asumo también cierto conocimiento en filtros digitales, aunque los bloques de estos son fáciles de usar y pueden consultar su documentación aquí] y tampoco estará de más tener muy claro el concepto de frecuencia en tiempo discreto. El flowgraph nos queda de la siguiente manera:

Click para agrandar

He agregado un slider gráfico que nos permite cambiar de forma interactiva la frecuencia de los generadores de señal que representan al oscilador local en cuadratura. Esto me parece un recurso didáctico muy bueno porque ayuda visualizar la relación que existe entre la frecuencia del oscilador local y la frecuencia que vemos en el centro (0 Hz) del espectro de la señal analítica. Basta asignar una frecuencia de 14.5 KHz para sintonizar a la componente de esa misma frecuencia que existe en nuestra señal de prueba:

Descargar flowgraph

Entendiendo a las señales analíticas (muestras I/Q)

Las muestras I/Q son una representación matemática discreta de una señal de radiofrecuencia cuya característica más útil es tener un espectro de frecuencia negativa nulo (un espectro compacto y asimétrico). Dado que el procedimiento de obtención de estas muestras involucra una reducción de frecuencia, su uso esta muy extendido en en las aplicaciones de radio definido por software y comunicaciones digitales en general. Las dos características anteriores nos dan dos ventajas: la primera elimina la redundancia en el espectro y la segunda permite reducir en el ancho de banda.

¿Cómo se obtiene una representación compleja de una señal real? Para el caso continuo el proceso es el siguiente:

 Simplemente se multiplica a la señal de radiofrecuencia por un sinusoide complejo y se aplica un filtro pasa-bajas para eliminar la redundancia en el espectro. Esto es lo que se conoce como señal analítica. Este procedimiento se vuelve más claro en el dominio de la frecuencia:

Recordemos que una multiplicación en el dominio del tiempo se manifiesta como una convolución en el dominio de la frecuencia. La convolución entre el espectro asimétrico de un sinuosoide complejo produce un corrimiento en frecuencia, el cual está dado por la frecuencia del sinusoide complejo (que lleva la etiqueta de LO por Local Oscillator) .

Ahora bien, el proceso anterior es puramente matemático. ¿Cómo se hace en el mundo real? El problema se resuelve implementando un modulador en cuadratura como se muestra en el siguiente diagrama:

 Tomando en cuenta a la identidad de Euler, una multiplicación por un sinusoide complejo se emula como la multiplicación y filtrado separado de dos sinusoides reales en cuadratura (heterodino en cuadratura). Posteriormente, al digitalizar independientemente a cada una de las componentes, se consideran a las muestras del ADC en cuadratura como la parte imaginaria multiplica por -1. La I viene de In-phase y la Q de Quadrature. Prácticamente todos los dongles SDR que pueden conseguirse fácilmente en el mercado como los RTL's y el FunCube funcionan de esta manera. En la siguiente entrada haremos un ejemplo con GNU Radio para clarificar aún más las cosas.

Referencias:
"RTL-SDR: Inexpensive Software Defined Radio", EE123 University of California Berkeley
"I/Q Data for Dummies", Mikael Q Kuisma
"Complex Sinusoids", PDS Related

Iniciando con GNU Radio

GNU Radio es un proyecto que provee módulos de Python (escritos en C++) y herramientas de software para el procesamiento de señales en aplicaciones de radio definido por software (SDR). Puede ser instalado en Windows aunque aún no lo he probado ahí. Instalarlo en Debían y Fedora es muy sencillo. Para este tutorial voy a asumir que se ha tomado un curso básico de procesamiento digital de señales (pero tanto para aprender como repasar se puede recurrir a la vieja confiable serie Schaum's).

Hay dos maneras de utilizar GNU Radio: En modo gráfico creando diagramas a bloque con GNU Radio Companion (GRC) o creando scripts de Python desde cero. La ventaja de usar GRC es que permite desarrollar aplicaciones muy rápidamente y con un entorno amigable. Sin embargo, para aplicaciones más complejas, escribir scrips provee una flexibilidad más adecuada. Naturalmente iniciaremos usando GRC.

GRC es visualmente similar a Simulink o LabVIEW. Sin embargo trabaja de la siguiente manera: GRC crea un script de python (cuyo nombre por defecto es top_block.py)  a partir del diagrama a bloques creado. Este programa se encarga de ejecutar los bloques en el orden correcto y llamar a la librería gráfica para visualizar interactivamente los datos si es el caso. Gracias al GNU Block API, cada bloque se ejecuta en un thread lo da muy buen desempeño si tienes un procesador multinúcleo.

Bien, al abrir GRC (buscadolo en el menú o desde terminal ejecutando "gnuradio-companion") aparecerá el espacio de trabajo, la ventana de bloques a la derecha y una ventana de mensajes en la parte inferior:

 En todo nuevo proyecto nos aparecerán dos bloques por efecto en el espacio de trabajo: el bloque "Options" y un bloque de variable samp_rate. El bloque "Options" nos permite cambiar al nombre del script principal (del que hablaba más arriba), agregar autor, descripción, elegir la librería gráfica entre otras opciones. Las variables son muy útiles para facilitar el desarrollo de cualquier programa. La tasa de muestreo es tan importante en procesamiento digital de señales que es muy útil agregarla al espacio de trabajo por defecto.

Para nuestro primer ejemplo generaremos una señal cosenoidal de 10 KHz con una amplitud unitaria más ruido gaussiano de amplitud 0.1 y la visualizaremos con un instrumento virtual. Usaremos los siguientes bloques:

• Waveform Generators > Signal Source
• Wavefoem Generator > Noise Source
• Math Operators > Add
• Misc > Throttle
• Instrumentation > QT > QT GUI Sink

 y los conectaremos de la siguiente manera (la conexión se hace dando click en la salida del primer bloque y después en la entrada del segundo bloque):

Se notará que el color de los puertos será azul por defecto. Este color denota el tipo de dato (el azul representa a los complejos flotantes). El código de colores completo se puede ver yendo al menú superior en help > types. En este ejemplo debemos cambiar a tipo float en la configuración de todos los bloques (naranja).

El bloque "Throttle" se encarga de reducir la carga de trabajo del procesador limitando el flujo de datos a través del bloque de manera que no haya retrasos o congelamientos en la visualización gráfica de los datos. Es muy importante aclarar que este bloque sólo debe usarse para visualizar datos y nunca para trabajar con salidas físicas como audio o instrumentos USRP.

Finalmente, lo que queda es configurar los bloques "Signal Source" y "Noise Source" para que nos den los valores de salida deseados. Damos click en "Execute flow graph":

 Veremos ejecutarse la interfaz gráfica del instrumento virtual QT que incluye displays en dominio de la frecuencia y del tiempo así como un espectrograma y un diagrama de constelación:

 Referencias
 "Core concepts of GNU Radio", GNURadio.org

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