Lo que nunca me enseñaron: Los filtros de Cauer (parte I)

Introducción

Según cuenta Van Valkenburg (ver "Analog Filter Design" 1982, página 379), sobre 1935 los ingenieros de Bell Laboratories se llevaron una desagradable sorpresa. Su competencia alemana había sacado un teléfono al mercado que igualaba la calidad de sonido del suyo, pero usaba un inductor menos en los filtros.
En 1935 los inductores no eran lo que hoy. Un inductor pesaba su cuarto de kilo y costaba bastante. El que los alemanes hubieran conseguido eliminarlo sin deteriorar las características del dispositivo podría llevar la empresa americana a la ruina y no era época para tonterías con los alemanes.

Teléfono de 1930

Pero esa ruina no ocurrió. El inventor del método de diseño que hacía capaz ese ahorro, Wilhelm Cauer, estaba necesitado de dinero (supongo) y quiso vender sus patentes. Para eso dio unas conferencias y los muy avispados matemáticos de Bell Laboratories tomaron nota. Una familia de funciones matemáticas se mencionó varias veces: las funciones elípticas de Jacobi.

Cuenta la leyenda que durante las dos semanas siguientes el departamento entero de matemáticas de Bell Laboratories se encerró en la biblioteca pública de Nueva York, estudiando tales funciones. Si mal no recuerda Darlington, que lo sufrió en sus propias carnes, allí encontró el artículo original de Jacobi de 1829 escrito en latín. En él estaba todo: tablas, transformaciones, aproximaciones… ¡todo!

A fin de cuentas, Cauer fue estudiante de Hilbert en Göttingen. Dada la prominencia del maestro, seguramente había visto de sobra estas funciones elípticas y, simplemente, les había encontrado otra aplicación. Por ese motivo, a este tipo de filtros se los denomina filtros elípticos o filtros de Cauer.

David Hilbert y Wilhelm Cauer

Los filtros

Un filtro eléctrico no es más que un dispositivo que toma una señal eléctrica como la del teléfono y elimina algunas frecuencias. De esta manera podemos oír la voz sin los ruidos raros del ADSL introduciendo un filtro que elimine las frecuencias que usa el ADSL. Lógicamente, podremos usar otro filtro para quedarnos con el ADSL y quitar la voz. Así que, por un mismo cable y gracias a un par de filtros, tengo voz y datos.

Splitter de ADSL que no es más que un filtro para que no se oiga el ruido de los datos de ADSL en el audio del teléfono.

Esta no es más que una de las aplicaciones de los filtros. Hay muchísimas más: multiplexar, eliminar el aliasing, hacer efectos de sonido y ecualizar, compensar desfases, modular, demodular… y un sinfín más que no vamos a enumerar aquí porque lo que nos interesa es el trasfondo teórico de los filtros.

La función de transferencia y la respuesta en frecuencia

La forma que tienen los filtros de procesar las señales es mediante una función de transferencia [$H(s)$]. Esta función está definida en los números complejos y, para que pueda ser construida físicamente, ha de ser racional (P y Q son polinomios). [$$H(s)=\frac{P(s)}{Q(s)}$$] La magnitud de la función de transferencia en el eje imaginario es la respuesta en frecuencia. Esta respuesta es lo que nuestro filtro multiplicará a la amplitud cada frecuencia. Si su valor a una frecuencia es mayor que uno, esa frecuencia se amplifica. Si es menor que uno, se atenúa. Escribiremos esta magnitud así:[$$|H(j\omega)|$$] Donde [$\omega$] no es más que [$2\pi$] veces la frecuencia ([$\omega=2\pi f$]) para no tener que escribir [$2\pi$] cada vez que vayamos a usar un seno, coseno o exponencial compleja. La [$j$] es la constante imaginaria que, en ingeniería, se usa para no confundirla con la intensidad de corriente [$i$].

Es importante señalar aquí que para realizar el filtro necesitamos la función [$H(s)$] con [$s$] en todo el plano complejo, pero la respuesta en frecuencia que deseamos es sólo esa función en la recta imaginaria [$|H(j\omega)|\ \ $]. Gran parte del problema estriba en pasar de [$|H(j\omega)|\ \ $], lo que queremos, a [$H(s)$], cómo hacerlo.

Pero primero veamos qué queremos.

El filtro ideal y los filtros reales

Nos interesaría que la respuesta en frecuencia fuera ideal: todo lo que haya por encima de una frecuencia de corte [$\omega_c$] se elimina (se multiplica por 0 en la banda de rechazo) y todo lo que haya por debajo se preserva (se multiplica por 1 en la banda de paso). Eso significa que el [$|H(j\omega)|$] ideal tendría esta forma:

Un filtro ideal que elimina las frecuencias de la banda de rechazo y deja intactas las frecuencias de la banda de paso.

Sin embargo, esto no es posible porque las únicas funciones de transferencia que podemos realizar eléctricamente son funciones racionales (recordemos: cocientes de dos polinomios) y son funciones continuas excepto en asíntotas verticales. No es posible hacer el escalón discontinuo que vemos en la gráfica de arriba. Hay que transigir y permitir no tener exactamente un 1 o un 0 (a esto se le denomina rizado) y dejar una banda de transición entre las bandas de paso y de rechazo.

Forma de un filtro real donde se han especificado los requisitos del mismo: frecuencia de paso [$\omega_p$], frecuencia de rechazo [$\omega_s$], rizado de paso [$r_p$] y rizado de rechazo [$r_s$].

Ya no tenemos una frecuencia de corte, sino una frecuencia de paso [$\omega_p$] por debajo de la cual las frecuencias no se van a tocar mucho (entre [$1$] y [$r_p$]) y una frecuencia de rechazo [$\omega_s$] por encima de la cual no van a sobrevivir muchas frecuencias (como mucho se multiplica por [$r_s$] que será todo lo cercano a cero posible).

Nota: El subíndice s viene del inglés stopband.

Complicando para simplificar

Como vemos, todo el rango de trabajo está entre el [$0$] y el [$1$]. Si elevamos al cuadrado, seguimos estando entre [$0$] y [$1$]. Esto permite trabajar con la magnitud al cuadrado que es más sencilla de calcular y siempre es real positiva. [$$ |H(j\omega)|^2=H(j\omega) H^*(j\omega)=H(j \omega )H(-j \omega )$$] El asterisco indica complejo conjugado. La última igualdad saca provecho a que [$H(s)$] es racional con coeficientes reales y, por tanto, su parte imaginaria ha de ser hermítica.

Trabajar entre [$0$] e infinito es también más sencillo. ¿Por qué? Si el límite es el [$1$], nos podemos pasar de [$1$] y eso sería erróneo. Nos tendríamos que dedicar a comprobar que nuestros cálculos nunca superasen el uno y los complicaría. En cambio, si el límite es infinito, ¡no hay límite! Los cálculos no tienen necesidad de asegurar que haya un límite y se simplifican.

Forzaremos que la función [$|H(j\omega)|^2$] tenga esta forma: [$$ |H(j\omega)|^2=\frac{1}{1+F^2 ( \omega )}$$]Ahora, la función [$F^2 ( \omega )$] no tiene límite superior. A cambio, hemos de recalcular los rizados.

La misma respuesta en frecuencia de antes, pero con la transformación mencionada para trabajar entre [$0$] e infinito. Hay que hacer notar cómo se han modificado los parámetros que especificaban los rizados.

Estos nuevos parámetros [$\epsilon_s$] y [$\epsilon_p$] se relacionan con los rizados originales así: [$$r_p=\sqrt{\frac{1}{1+\epsilon_p^2}},\ \ \ \ r_s=\sqrt{\frac{\epsilon_s^2}{1+ \epsilon _s^2}}$$] Llamaremos a [$\omega_p,\epsilon_p,\omega_s$] y [$\epsilon_s$] los requisitos del filtro.

Hay que destacar que la función [$F^2 (\omega)$] sigue siendo una función racional y deberá ser el cociente de dos polinomios.

Continúa en parte dos.

miniSL parte 14 - La Evaluación

Retomamos el lenguaje de script miniSL. Hasta ahora hemos definido los datos sobre los que vamos a tratar (las celdas) y hemos implementado las operaciones básicas sobre ellos: creación y destrucción. También hemos destacado algunos datos como el entorno global.

Ahora vamos a implementar la principal operación a realizar sobre los datos: tomar una celda que represente una expresión y calcular su valor (que será otra celda). Esto es lo que se llama la evaluación.

Existen varios modelos de evaluación. Aquí nos centraremos en la evaluación por entornos (se explica aquí y aquí). Este tipo de evaluación es realmente sencilla: si queremos saber el valor de una variable, hemos de buscarla en el entorno actual. Si no, buscamos en el entorno padre y así sucesivamente. Esto lo hace la función FindName() que se vio en la parte 10.

La evaluación de literales es inmediata: representan su propio valor (el 5 es el 5, la cadena “hola” es la cadena “hola”). La evaluación de combinaciones (aplicaciones de una función a argumentos) procede de la siguiente manera.

1. Se evalúa el operador (la función).
2. Se envían los operandos sin evaluar a la función.

Como hay dos tipos de funciones: la definida por el usuario y la nativa, tendremos que contemplar ambos casos. En principio, la definida por el usuario evaluará los operandos; pero veremos que las nativas necesitan sus operandos sin evaluar (en las partes 29,30 y siguientes).

La evaluación de una lista es la lista de sus miembros evaluados. Así, [1+2,5] es la lista [3,5]. La lista vacía es precisamente un literal porque no hay que evaluar nada para calcular su valor.

Cualquier otra celda no es considerada expresión y genera un error de ejecución.

El código de la función de evaluación es el siguiente:

CELL& Script::Evaluate(CELL& c, CELL& envir)
{
 CELL* aux;
 switch(c.type)
 {
  //Non evaluable
 case UNUSED: throw L"Evaluating an unused cell";
 default:  throw L"Evaluating an unknown cell";
 case LAMBDA_VAL:throw L"Evaluating a lambda value";
 case NATIVE_VAL:throw L"Evaluating a native";
 case ENVIR_VAL: throw L"Evaluating an environment";

  //Literals
 case INT_LIT: case BOOL_VAL: case STRING_LIT: case EMPTY_LIT:
  return c;

  //Code constructors
 case CONS_CTOR:  return CreateCell(CONS_CTOR, &Evaluate(*c.head, envir), &Evaluate(*c.tail, envir));

 case NAME_CODE: 
  if((aux=FindName(c, envir))==NULL)
   throw L"Unknown name";
  return *aux;

 case COMBINE_CODE:
  switch((aux=&Evaluate(*c.op, envir))->type)
  {
  case LAMBDA_VAL: return ApplyLambda(*aux->code, *c.operands, *aux->closure, envir);
  case NATIVE_VAL: return aux->native(*this, *c.operands, envir);
  default:   throw L"Non-combinable value";
  }
 }
}

Lo más complejo es la evaluación de la combinación. Esto es debido a los casos antes mencionados. Cuando tenemos una operación nativa, directamente llamamos a la función que la implementa. Si tenemos una operación definida por el usuario (lambda) usaremos una función auxiliar ApplyLambda() que veremos en la siguiente parte de esta serie.

The vessel with the pestle

Como me he tomado algunas vacaciones, estoy buscando y viendo algunas de las películas que me gustaron de niño. Las veo en versión original para practicar el inglés. Realmente me he encontrado con algunas sorpresas muy agradables. Por ejemplo, he encontrado esta joya.

Activad los subtítulos donde pone CC.

Algunas palabrejas poco usuales del inglés que se usan en este contexto son:

• Vessel: vasija
• Pestle: mortero
• Pellet: pastilla
• Flagon: jarra
• Brew: brebaje

Números binarios en C++

ATENCIÓN ACTUALIZACIÓN (2014-12-13): El estándar de C++ de 2014 incluye la posibilidad de usar literales binarios directamente con la notación 0b1001.


En C++ (y en C y en otros muchos lenguajes de programación) hay tres formas de escribir números enteros literales.

1. En base 10 (decimal) cuando escribo algo como 1042348
2. En base 16 (hexadecimal) cuando prefijo 0x como en 0x3FFF
3. En base 8 (octal) cuando prefijo sólo 0 como en 0663

Cuando se programa en bajo nivel muchas veces es más directo escribir directamente en base 2 (binario), ¿pero cómo? Me he topado con la siguiente macro que lo consigue

#define Ob(x)  ((unsigned)Ob_(0 ## x ## uL))
#define Ob_(x) (x & 1 | x >> 2 & 2 | x >> 4 & 4 | x >> 6 & 8 |  \
 x >> 8 & 16 | x >> 10 & 32 | x >> 12 & 64 | x >> 14 & 128)

El truco se basa en dos pasos, el primero Ob (no es un cero, es una o) pone delante un 0 (que sí es un cero) y detrás un "uL". El "uL" es la forma que tenemos de decirle al C++ que el número es sin signo y largo. El cero, por otra parte, sí que es importante ya que nos va a convertir el número a octal. Lo que tenemos entonces es un número binario interpretado como un número octal. ¿Cómo afecta eso a su valor?

Un número binario se expresa de la siguiente forma [$$ n=\cdots + d_3 2^3 + d_2 2^2 + d_1 2^1 + d_0 2^0 $$] donde los [$d_i$] son los dígitos. Los dígitos serán o cero o uno. Al escribir el número binario como octal lo que hacemos es  [$$ n=\cdots + d_3 8^3 + d_2 8^2 + d_1 8^1 + d_0 8^0 $$] pero aún con los dígitos en binario. Así que lo único que necesitamos es leer los bits de tres en tres ya que [$8=2^3$].  Justamente eso es lo que hace Ob_.

Por ejemplo:
Ob(1001) //Como binario es 9 en decimal
Ob_(01001uL) //Como octal es 513 en decimal y 1_000_000_001 en binario
9 //Tomando un bit de cada 3 en el 513 

Clasificación de los bugs

Leyendo a Walter Bright me he encontrado con la palabreja "heisenbug" que según él es un tipo de bug. Ya me imaginaba yo que tendría que haber una taxonomía de los bugs. El resumen de lo que he encontrado es este:

• Bohrbug: Error sistemático que se da siempre cuando las condiciones lo propician.
• Mandelbug: Error sistemático tan complejo que no parece sistemático, pero lo es.
• Heisenbug: Error sistemático que desaparece cuando intentamos depurarlo. Ya sea por que estamos usando el depurador y cambia la estructura de la memoria, porque el código está compilado en modo de depuración o por otras razones.
• Schrödinbug: Error oculto que, una vez descubierto, todo el mundo se tropieza con él.
• Estadístico: Error que no se detecta en una única ejecución del programa, sino que hay que hacer una estadística de los resultados hasta darnos cuenta que no son como deberían ser. (Esto me recuerda a la minería de datos en los MMORPG)
• Alfa bug: Error que ves una vez y no vuelves a ver. Lo que ocurría antiguamente cuando un rayo cósmico le daba a una celda de memoria.

Supongo que luego están los errores de condiciones de carrera, que aparecen cuando les da la gana si no se han sincronizado bien las hebras. Ya inventarán una palabreja para eso algún día.

El lío de las composiciones

La composición en las relaciones

Una relación [$\le$] entre [$A$] y [$B$] es un subconjunto del conjunto cartesiano.[$$\le \subseteq A\times B$$]Cuando dos elementos [$a\in A$] y [$b\in B$] están relacionados escribimos [$$ a \le b \Leftrightarrow (a,b)\in \le$$] La composición de dos relaciones [$\le_1$] entre [$A$] y [$B$]; y [$\le_2$] entre [$B$] y [$C$] es [$\le_1 \circ \le_2$] definida así [$$ a \le_1\circ \le_2 c \Leftrightarrow \exists b\in B.a\le_1 b\le_2 c $$] Es interesante ver el parecido natural entre [$$  a \le_1\circ \le_2 c$$] y [$$ a\le_1 b\le_2 c $$] como si el [$\circ$] representase un elemento indefinido a buscar.

La composición en las funciones

Una función [$f:A\to B$] no es más que una relación en la cual [$$\forall a\in A.\exists ! b\in B. a f b $$] Para los que no lo hayan visto nunca [$\exists !$] significa "existe un único".
Bien, los problemas empiezan porque si existe un único [$b$] lo natural es representarlo por [$f$] y [$a$]. La forma de hacerlo es esta: [$$f(a)=b \Leftrightarrow a f b$$] ¿Por qué es problemático? Porque si pienso que una función es una relación y defino la composición de funciones como una composición de relaciones obtenemos [$$ (f\circ g)(a)=g(f(a)) $$] Invirtiéndose el orden de [$f$] y [$g$], lo que lía muchísimo. Así que lo que se hace es definir al revés la composición de funciones. [$$ (f\circ g)(a)=f(g(a)) $$] Por lo que deja de ser igual que la composición de relaciones.

Soluciones

Las que se me ocurren así a bote pronto son:

1. Usar una notación postfija para las funciones [$(a)f$] en vez de [$f(a)$]. De esta manera [$ (a)(f\circ g)  =((a)f)g$] y el orden coincide con la composición de relaciones.
2. Invertir una de las dos composiciones. O bien [$ (f\circ g)(a)=g(f(a)) $] o bien [$ a \le_2\circ \le_1 c \Leftrightarrow \exists b\in B.a\le_1 b\le_2 c $]. Hagas lo que hagas, estás invirtiendo el orden natural de la notación y vas a liarte (el elemento [$a$] va primero a [$f$] y a [$\le_1$] que están lejos de él).
3. Usar dos notaciones. Por ejemplo, usar [$\le_1 ; \le_2 = \le_2 \circ \le_1$] y [$ f \circ g = g ; f $]. Esta solución la suele usar alguno que otro autor.
4. Definir las funciones al revés, [$\forall a\in A.\exists ! b\in B. b f a $] de forma que   [$b=f(a) \Leftrightarrow b f a$]. Claro que también así vamos en contra de unos cuantos siglos de historia y habría que darle la vuelta a muchas definiciones como relaciones inyectivas.

Exponenciando la derivación

Realmente para definir la exponencial (o cualquier otra función analítica) mediante una serie de potencias hacen falta pocas cosas.

• Una operación de multiplicación por escalar
• Una operación de suma
• Una operación de potencia

Con estas tres operaciones podemos definir
[$$ e^x = \sum_{k=0}^\infty \frac{x^k}{k!} $$]
Los operadores lineales cumplen todas estas condiciones tomando la potencia como composición iterada.
[$$f^0=Id;  f^{n+1}=f\circ f^n$$]
Un caso clásico es cuando las operaciones lineales están representadas por matrices.
Otro de estos operadores lineales es la derivación, que escribiremos [$D$]. Se puede realizar la exponencial de la derivación.
[$$e^D = \sum_{k=0}^\infty \frac{D^k}{k!} $$]
Para ver mejor sus efectos, aplicaremos el resultado de [$\alpha$] veces esta exponencial sobre [$x^n$].
[$$e^{\alpha D} x^n = \sum_{k=0}^\infty \frac{\alpha^k D^k}{k!} x^n $$]
Podemos derivar hasta [$n$] veces [$x^n$]
[$$\sum_{k=0}^n \frac{n!  \alpha^k x^{n-k} }{k! (n-k)!}$$]
Es interesantísimo ver que lo obtenido es justamente el binomio de Newton.
[$$ \sum_{k=0}^n \frac{n!  \alpha^k x^{n-k} }{k! (n-k)!} = (x+\alpha)^n $$]
Debido a que estamos trabajando con operadores lineales, podemos volver a usar una serie de potencias para aplicar la exponencial de la derivación sobre una función analítica arbitraria.
[$$ (e^{\alpha D} f)(x) = e^{\alpha D} \sum_{i=0}^\infty {a_i x^i} =   \sum_{i=0}^\infty {a_i   e^{\alpha D} x^i} =  \sum_{i=0}^\infty {a_i  (x+\alpha)^i} = f(x+\alpha) $$]
Así que la exponencial de la derivación no es más que una traslación.
Los físicos me dicen que acabo de descubrir que el operador de momento genera las traslaciones. ¡Vaya! Pues no lo sabía.