Transformada Hartley discreta

Transformada matemática relacionada con Fourier

Una transformada Hartley discreta (DHT) es una transformada relacionada con Fourier de datos discretos y periódicos similar a la transformada de Fourier discreta (DFT), con aplicaciones análogas en el procesamiento de señales y campos relacionados. Su principal distinción con la DFT es que transforma entradas reales en salidas reales, sin la participación intrínseca de números complejos . Así como la DFT es el análogo discreto de la transformada de Fourier continua (FT), la DHT es el análogo discreto de la transformada de Hartley continua (HT), introducida por Ralph VL Hartley en 1942. ^[1]

Dado que existen algoritmos rápidos para la DHT análogos a la transformada rápida de Fourier (FFT), la DHT fue propuesta originalmente por Ronald N. Bracewell en 1983 ^[2] como una herramienta computacional más eficiente en el caso común en el que los datos son puramente reales. Sin embargo, posteriormente se argumentó que los algoritmos FFT especializados para entradas o salidas reales se pueden encontrar normalmente con un poco menos de operaciones que cualquier algoritmo correspondiente para la DHT.

Definición

Formalmente, la transformada discreta de Hartley es una función lineal invertible H : R ⁿ → R ⁿ (donde R denota el conjunto de números reales ). Los N números reales x ₀ , ..., x _{N −1} se transforman en los N números reales H ₀ , ..., H _{N −1} según la fórmula

${\displaystyle H_{k}=\sum _{n=0}^{N-1}x_{n}\operatorname {cas} \left({\frac {2\pi }{N}}nk\right)=\sum _{n=0}^{N-1}x_{n}\left[\cos \left({\frac {2\pi }{N}}nk\right)+\sin \left({\frac {2\pi }{N}}nk\right)\right]\quad \quad k=0,\dots ,N-1.}$

La combinación a veces se denota cas( z ) , y no debe confundirse con cis( z ) = e ^iz = cos( z ) + i  sin( z ) , o e ^{− iz} = cis(− z ) que aparece en la definición de DFT (donde i es la unidad imaginaria ). ${\displaystyle \cos(z)+\sin(z)}$ ${\displaystyle ={\sqrt {2}}\cos \left(z-{\frac {\pi }{4}}\right)}$

Al igual que con la DFT, el factor de escala general delante de la transformación y el signo del término seno son una cuestión de convención. Aunque estas convenciones varían ocasionalmente entre autores, no afectan las propiedades esenciales de la transformación.

Propiedades

La transformada puede interpretarse como la multiplicación del vector ( x ₀ , ...., x _{N −1} ) por una matriz N -por- N ; por lo tanto, la transformada discreta de Hartley es un operador lineal . La matriz es invertible; la transformación inversa, que permite recuperar la x _n a partir de la H _k , es simplemente la DHT de H _k multiplicada por 1/ N . Es decir, la DHT es su propia inversa ( involutiva ), hasta un factor de escala global.

La DHT se puede utilizar para calcular la DFT, y viceversa. Para entradas reales x _n , la salida de la DFT X _k tiene una parte real ( H _k + H _N−k )/2 y una parte imaginaria ( H _N−k − H _k )/2. Por el contrario, la DHT es equivalente a calcular la DFT de x _n multiplicada por 1 +  i , y luego tomar la parte real del resultado.

Al igual que con la DFT, una convolución cíclica z = x ∗ y de dos vectores x = ( x _n ) e y = ( y _n ) para producir un vector z = ( z _n ), todos de longitud N , se convierte en una operación simple después de la DHT. En particular, supongamos que los vectores X , Y y Z denotan la DHT de x , y y z respectivamente. Entonces los elementos de Z están dados por:

${\displaystyle {\begin{matrix}Z_{k}&=&\left[X_{k}\left(Y_{k}+Y_{N-k}\right)+X_{N-k}\left(Y_{k}-Y_{N-k}\right)\right]/2\\Z_{N-k}&=&\left[X_{N-k}\left(Y_{k}+Y_{N-k}\right)-X_{k}\left(Y_{k}-Y_{N-k}\right)\right]/2\end{matrix}}}$

donde tomamos todos los vectores como periódicos en N ( X _N = X ₀ , etcétera). Por lo tanto, al igual que la DFT transforma una convolución en una multiplicación puntual de números complejos ( pares de partes reales e imaginarias), la DHT transforma una convolución en una combinación simple de pares de componentes de frecuencia real. La DHT inversa produce entonces el vector deseado z . De esta manera, un algoritmo rápido para la DHT (ver más abajo) produce un algoritmo rápido para convolución. (Esto es ligeramente más caro que el procedimiento correspondiente para la DFT, sin incluir los costos de las transformaciones a continuación, porque la operación por pares anterior requiere 8 operaciones aritméticas reales en comparación con las 6 de una multiplicación compleja. Este recuento no incluye la división por 2, que puede ser absorbida, por ejemplo, en la normalización 1/ N de la DHT inversa).

Algoritmos rápidos

Al igual que para la DFT, evaluar la definición de DHT directamente requeriría O( N ² ) operaciones aritméticas (ver notación Big O ). Sin embargo, existen algoritmos rápidos similares a la FFT que calculan el mismo resultado en solo O( N log N ) operaciones. Casi todos los algoritmos FFT, desde Cooley-Tukey hasta el factor primo , pasando por Winograd (1985) ^[3] y Bruun (1993), ^[4] tienen un análogo directo para la transformada discreta de Hartley. (Sin embargo, algunos de los algoritmos FFT más exóticos, como la QFT, aún no se han investigado en el contexto de la DHT).

En particular, el análogo DHT del algoritmo Cooley–Tukey se conoce comúnmente como el algoritmo de transformada rápida de Hartley (FHT), y fue descrito por primera vez por Bracewell en 1984. ^[5] Este algoritmo FHT, al menos cuando se aplica a tamaños de potencia de dos N , es el tema de la patente de los Estados Unidos número 4.646.256, otorgada en 1987 a la Universidad de Stanford . Stanford colocó esta patente en el dominio público en 1994 (Bracewell, 1995). ^[6]

Como se mencionó anteriormente, los algoritmos DHT son típicamente ligeramente menos eficientes (en términos de la cantidad de operaciones de punto flotante ) que el algoritmo DFT correspondiente (FFT) especializado para entradas (o salidas) reales. Esto fue argumentado por primera vez por Sorensen et al. (1987) ^[7] y Duhamel & Vetterli (1987). ^[8] Los últimos autores obtuvieron lo que parece ser el recuento de operaciones publicado más bajo para el DHT de tamaños de potencia de dos, empleando un algoritmo de base dividida (similar a la FFT de base dividida ) que divide un DHT de longitud N en un DHT de longitud N /2 y dos DFT de entrada real ( no DHT) de longitud N /4. De esta manera, argumentaron que un DHT de longitud de potencia de dos se puede calcular con, en el mejor de los casos, 2 adiciones más que el número correspondiente de operaciones aritméticas para el DFT de entrada real.

En las computadoras actuales, el rendimiento está determinado más por consideraciones de caché y pipeline de CPU que por recuentos estrictos de operaciones, y es poco probable que una ligera diferencia en el costo aritmético sea significativa. Dado que los algoritmos FHT y FFT de entrada real tienen estructuras computacionales similares, ninguno parece tener una ventaja de velocidad sustancial a priori (Popović  [sr] y Šević, 1994). ^[9] Como cuestión práctica, las bibliotecas FFT de entrada real altamente optimizadas están disponibles de muchas fuentes (por ejemplo, de proveedores de CPU como Intel ), mientras que las bibliotecas DHT altamente optimizadas son menos comunes.

Por otra parte, los cálculos redundantes en las FFT debido a las entradas reales son más difíciles de eliminar para números primos grandes N , a pesar de la existencia de algoritmos de datos complejos O( N log N ) para tales casos, porque las redundancias están ocultas detrás de permutaciones intrincadas y/o rotaciones de fase en esos algoritmos. Por el contrario, un algoritmo FFT de tamaño primo estándar, el algoritmo de Rader , se puede aplicar directamente a la DHT de datos reales para aproximadamente un factor de dos menos de cálculo que el de la FFT compleja equivalente (Frigo y Johnson, 2005). ^[10] Por otra parte, también es posible una adaptación no basada en DHT del algoritmo de Rader para DFT de entrada real (Chu y Burrus , 1982). ^[11]

Transformada Hartley discreta multidimensional (MD-DHT)

La rD-DHT (MD-DHT con dimensiones "r") viene dada por

${\displaystyle X(k_{1},k_{2},...,k_{r})=\sum _{n_{1}=0}^{N_{1}-1}\sum _{n_{2}=0}^{N_{2}-1}\dots \sum _{n_{r}=0}^{N_{r}-1}x(n_{1},n_{2},...,n_{r}){\rm {cas}}({\frac {2\pi n_{1}k_{1}}{N_{1}}}+\dots +{\frac {2\pi n_{r}k_{r}}{N_{r}}}),}$

con y donde ${\displaystyle k_{i}=0,1,\ldots ,N_{i}-1}$ ${\displaystyle {\rm {cas}}(x)=\cos(x)+\sin(x).}$

De manera similar al caso 1-D, como transformación real y simétrica, la MD-DHT es más simple que la MD-DFT. Por un lado, la DHT inversa es idéntica a la transformación directa, con la adición de un factor de escala;

y segundo, dado que el núcleo es real, evita la complejidad computacional de los números complejos . Además, la DFT se puede obtener directamente de la DHT mediante una simple operación aditiva (Bracewell, 1983). ^[2]

El MD-DHT se utiliza ampliamente en áreas como el procesamiento de imágenes y señales ópticas. Las aplicaciones específicas incluyen visión artificial, televisión de alta definición y teleconferencias, áreas que procesan o analizan imágenes en movimiento (Zeng, 2000). ^[12]

Algoritmos rápidos para el MD-DHT

A medida que aumenta la velocidad de procesamiento, los problemas multidimensionales más grandes se vuelven factibles computacionalmente, lo que requiere algoritmos multidimensionales rápidos. A continuación se presentan tres de esos algoritmos.

En busca de la separabilidad para la eficiencia, consideramos la siguiente transformación (Bracewell, 1983), ^[2]

${\displaystyle {\hat {X}}(k_{1},k_{2},...,k_{r})=\sum _{n_{1}=0}^{N_{1}-1}\sum _{n_{2}=0}^{N_{2}-1}\dots \sum _{n_{r}=0}^{N_{r}-1}x(n_{1},n_{2},...,n_{r}){\rm {cas}}({\frac {2\pi n_{1}k_{1}}{N_{1}}})\dots {\rm {cas}}({\frac {2\pi n_{r}k_{r}}{N_{r}}}).}$

^{Bortfeld (1995) [13]} demostró que ambos pueden relacionarse mediante algunas adiciones. Por ejemplo, en 3-D,

${\displaystyle X(k_{1},k_{2},k_{3})={\frac {1}{2}}[{\hat {X}}(k_{1},k_{2},-k_{3})+{\hat {X}}(k_{1},-k_{2},k_{3})+{\hat {X}}(-k_{1},k_{2},k_{3})-{\hat {X}}(-k_{1},-k_{2},-k_{3})].}$

Para , se pueden implementar algoritmos de fila-columna. Esta técnica se utiliza comúnmente debido a la simplicidad de dichos algoritmos RC, pero no están optimizados para espacios MD generales. ${\displaystyle {\hat {X}}}$

Se han desarrollado otros algoritmos rápidos, como el de base 2, base 4 y base dividida. Por ejemplo, Boussakta (2000) ^[14] desarrolló el vector de base 3D,

${\displaystyle X(k_{1},k_{2},...,k_{r})=\sum _{n_{1}=0}^{N-1}\sum _{n_{2}=0}^{N-1}\sum _{n_{r}=0}^{N-1}x(n_{1},n_{2},n_{3}){\rm {cas}}({\frac {2\pi }{N}}(n_{1}k_{1}+n_{2}k_{2}+n_{3}k_{3}))}$

${\displaystyle =\sum _{n_{1}:even}\sum _{n_{2}:even}\sum _{n_{3}:even}+\sum _{n_{1}:even}\sum _{n_{2}:even}\sum _{n_{3}:odd}+\sum _{n_{1}:even}\sum _{n_{2}:odd}\sum _{n_{3}:even}}$

${\displaystyle +\sum _{n_{1}:even}\sum _{n_{2}:odd}\sum _{n_{3}:odd}+\sum _{n_{1}:odd}\sum _{n_{2}:even}\sum _{n_{3}:even}+\sum _{n_{1}:odd}\sum _{n_{2}:even}\sum _{n_{3}:odd}}$

${\displaystyle +\sum _{n_{1}:odd}\sum _{n_{2}:odd}\sum _{n_{3}:even}+\sum _{n_{1}:odd}\sum _{n_{2}:odd}\sum _{n_{3}:odd}.}$

^{En Boussakta (2000) [14]} también se presentó que este algoritmo de base vectorial 3D toma multiplicaciones y sumas en comparación con las multiplicaciones y sumas del enfoque fila-columna. El inconveniente es que la implementación de estos algoritmos de base es difícil de generalizar para señales de dimensiones arbitrarias. ${\displaystyle ({\frac {7}{4}})N^{3}\log _{2}N}$ ${\displaystyle ({\frac {31}{8}})N^{3}\log _{2}N}$ ${\displaystyle 3N^{3}\log _{2}N}$ ${\displaystyle ({\frac {9}{2}})N^{3}\log _{2}N+3N^{2}}$

Las transformadas de teoría de números también se han utilizado para resolver la MD-DHT, ya que realizan convoluciones extremadamente rápidas. En Boussakta (1988), ^[15] se mostró cómo descomponer la transformada MD-DHT en una forma que consta de convoluciones:

Para el caso 2-D (el caso 3-D también está cubierto en la referencia indicada),

${\displaystyle X(k,l)=\sum _{n=0}^{N-1}\sum _{m=0}^{M-1}x(n,m){\rm {cas}}({\frac {2\pi nk}{N}}+{\frac {2\pi ml}{M}}),\;}$ ${\displaystyle k=0,1,\ldots ,N-1}$, ${\displaystyle l=0,1,\ldots ,M-1}$

se puede descomponer en convoluciones circulares 1-D y 2-D de la siguiente manera,

${\displaystyle X(k,l)={\begin{cases}X_{1}(k,0)\\X_{2}(0,l)\\X_{3}(k,l)\end{cases}}}$

dónde

${\displaystyle X_{1}(k,0)=\sum _{n=0}^{N-1}(\sum _{m=0}^{M-1}x(n,m)){\rm {cas}}({\frac {2\pi nk}{N}}),\;}$ ${\displaystyle k=0,1,\ldots ,N-1}$

${\displaystyle X_{2}(0,l)=\sum _{m=0}^{M-1}(\sum _{n=0}^{N-1}x(n,m)){\rm {cas}}({\frac {2\pi ml}{M}}),\;}$ ${\displaystyle l=1,2,\dots ,M-1}$

${\displaystyle X_{3}(k,l)=\sum _{n=0}^{N-1}\sum _{m=0}^{M-1}x(n,m){\rm {cas}}({\frac {2\pi nk}{N}}+{\frac {2\pi ml}{M}})\;,}$

${\displaystyle k=1,2,\ldots ,N-1}$

${\displaystyle l=1,2,\ldots ,M-1.}$

Desarrollando más, ${\displaystyle X_{3}}$

${\displaystyle X_{3}(k,l)=\sum _{n=0}^{N-1}x(n,0){\rm {cas}}({\frac {2\pi nk}{N}})+\sum _{m=1}^{M-1}x(0,m){\rm {cas}}({\frac {2\pi ml}{M}})}$

${\displaystyle +\sum _{n=1}^{N-1}\sum _{m=1}^{M-1}x(n,m){\rm {cas}}({\frac {2\pi nk}{N}}+{\frac {2\pi ml}{M}}).}$

En este punto presentamos la transformada del número de Fermat (FNT). El t ^-ésimo número de Fermat está dado por , con . Los números de Fermat bien conocidos son para ( es primo para ), (Boussakta, 1988). ^[15] La transformada del número de Fermat está dada por ${\displaystyle F_{t}=2^{b}+1}$ ${\displaystyle b=2^{t}}$ ${\displaystyle t=0,1,2,3,4,5,6}$ ${\displaystyle F_{t}}$ ${\displaystyle 0\leq t\leq 4}$

${\displaystyle X(k,l)=\sum _{n=0}^{N-1}\sum _{m=0}^{M-1}x(n,m)\alpha _{1}^{nk}\alpha _{2}^{ml}\mod F_{t}}$

con . y son raíces de unidad de orden y respectivamente . ${\displaystyle k=0,\ldots ,N-1,l=0,\ldots ,M-1}$ ${\displaystyle \alpha _{1}}$ ${\displaystyle \alpha _{2}}$ ${\displaystyle N}$ ${\displaystyle M}$ ${\displaystyle (\alpha _{1}^{N}=\alpha _{2}^{M}=1\mod F_{t})}$

Volviendo a la descomposición, el último término para se denotará como , entonces ${\displaystyle X_{3}(k,l)}$ ${\displaystyle X_{4}(k,l)}$

${\displaystyle X_{4}(k,l)=\sum _{n=1}^{N-1}\sum _{m=1}^{M-1}x(n,m){\rm {cas}}({\frac {2\pi nk}{N}}+{\frac {2\pi ml}{M}}),}$

${\displaystyle k=1,2,\ldots ,N-1}$

${\displaystyle l=1,2,\ldots ,M-1.}$

Si y son raíces primitivas de y (cuya existencia está garantizada si y son primos ), entonces y se asignan a . Por lo tanto, al asignar y a y , se obtiene lo siguiente: ${\displaystyle g_{1}}$ ${\displaystyle g_{2}}$ ${\displaystyle N}$ ${\displaystyle M}$ ${\displaystyle M}$ ${\displaystyle N}$ ${\displaystyle g_{1}}$ ${\displaystyle g_{2}}$ ${\displaystyle (n,m)}$ ${\displaystyle (g_{1}^{n}\mod N,g_{2}^{m}\mod M).}$ ${\displaystyle n,m,k}$ ${\displaystyle l}$ ${\displaystyle g_{1}^{-n},g_{2}^{-m},g_{1}^{k}}$ ${\displaystyle g_{2}^{l}}$

${\displaystyle X_{4}(g_{1}^{k},g_{2}^{l})=\sum _{n=0}^{N-2}\sum _{m=0}^{M-2}x(g_{1}^{-n},g_{2}^{-m}){\rm {cas}}({\frac {2\pi g_{1}^{(-n+k)}}{N}}+{\frac {2\pi g_{2}^{(-m+l)}}{M}}),}$

${\displaystyle k=0,1,\ldots ,N-2}$

${\displaystyle l=0,1,\ldots ,M-2}$.

Lo cual es ahora una convolución circular . Con , , y , se tiene ${\displaystyle Y(k,l)=X_{4}(g_{1}^{k},g_{2}^{l})}$ ${\displaystyle y(n,m)=x(g_{1}^{-n},g_{2}^{-m})}$ ${\displaystyle h(n,m)={\rm {cas}}({\frac {2\pi g_{1}^{n}}{N}}+{\frac {2\pi g_{2}^{m}}{M}})}$

${\displaystyle Y(k,l)=\sum _{n=0}^{N-2}\sum _{m=0}^{M-2}y(n,m)h(<k-n>_{N},<l-m>_{M})}$

${\displaystyle Y(k,l)=FNT^{-1}\{FNT[y(n,m)]\otimes FNT[h(n,m)]}$

donde denota la multiplicación término por término. También se afirmó en (Boussakta, 1988) ^[15] que este algoritmo reduce el número de multiplicaciones en un factor de 8 a 20 con respecto a otros algoritmos DHT a costa de un ligero aumento en el número de operaciones de desplazamiento y adición, que se supone que son más simples que las multiplicaciones. La desventaja de este algoritmo es la restricción de que cada dimensión de la transformación tiene una raíz primitiva . ${\displaystyle \otimes }$

Referencias

1. Hartley, Ralph VL (marzo de 1942). "Un análisis de Fourier más simétrico aplicado a problemas de transmisión". Actas del IRE . 30 (3): 144–150. doi :10.1109/JRPROC.1942.234333. S2CID  51644127.
2. Bracewell, Ronald N. (1983). "Transformada Hartley discreta". Revista de la Sociedad Óptica de América . 73 (12): 1832–1835. doi :10.1364/josa.73.001832. S2CID  120611904.
3. Sorensen, Henrik V.; Jones, Douglas L .; Burrus, Charles Sidney ; Heideman, Michael T. (1985). "Sobre el cálculo de la transformada discreta de Hartley". IEEE Transactions on Acoustics, Speech, and Signal Processing . ASSP-33 (4): 1231–1238.
4. Bini, Dario Andrea; Bozzo, Enrico (1993). "Transformada discreta rápida mediante polinomios propios". Computers & Mathematics with Applications . 26 (9): 35–52. doi : 10.1016/0898-1221(93)90004-f .
5. Bracewell, Ronald N. (1984). "La transformada rápida de Hartley". Actas del IEEE . 72 (8): 1010–1018. doi :10.1109/proc.1984.12968. S2CID  21988816.
6. Bracewell, Ronald N. (1995). "Computación con la transformada de Hartley". Computers in Physics . 9 (4): 373–379. Bibcode :1995ComPh...9..373B. doi : 10.1063/1.168534 .
7. Sorensen, Henrik V.; Jones, Douglas L .; Heideman, Michael T.; Burrus, Charles Sidney (1987). "Algoritmos de transformada rápida de Fourier con valores reales". IEEE Transactions on Acoustics, Speech, and Signal Processing . ASSP-35 (6): 849–863.
8. Duhamel, Pierre; Vetterli, Martin (1987). "Algoritmos mejorados de transformada de Fourier y Hartley: aplicación a la convolución cíclica de datos reales". IEEE Transactions on Acoustics, Speech, and Signal Processing . ASSP-35: 818–824.
9. Поповић [Popović], Миодраг [Miodrag] [en serbio] ; Šević, Dragutin (1994). "Una nueva mirada a la comparación de las transformadas rápidas de Hartley y Fourier". Transacciones IEEE sobre procesamiento de señales . 42 (8): 2178–2182. Código bibliográfico : 1994ITSP...42.2178P. doi : 10.1109/78.301854.
10. Frigo, Matteo; Johnson, Steven G. (2005). "El diseño y la implementación de FFTW3" (PDF) . Actas del IEEE . 93 (2): 216–231. CiteSeerX 10.1.1.66.3097 . doi :10.1109/jproc.2004.840301. S2CID  6644892. }
11. Chu, Shuni; Burrus, Charles Sidney (1982). "Un algoritmo de factor primo FTT [ sic ] que utiliza aritmética distribuida". IEEE Transactions on Acoustics, Speech, and Signal Processing . 30 (2): 217–227. doi :10.1109/tassp.1982.1163875.
12. Zeng, Yonghang; Bi, Guoan; Leyman, Abdul Rahim (2000). "Algoritmos de transformada polinómica para transformada Hartley discreta multidimensional". Simposio internacional IEEE sobre circuitos y sistemas (V): 517–520.
13. Bortfeld, Thomas; Dinter, Wolfgang (1995). "Cálculo de transformadas de Hartley multidimensionales utilizando transformadas de Fourier unidimensionales". IEEE Transactions on Signal Processing . 43 (5): 1306–1310. Bibcode :1995ITSP...43.1306B. doi :10.1109/78.382424.
14. Boussakta, Said; Alshibami, Osama (2000). "Algoritmo rápido para la transformada Hartley discreta tridimensional". Conferencia internacional sobre acústica, habla y procesamiento de señales '00 (4): 2302–2305.
15. Boussakta, Said; Holt, Alan GJ (1988). "Transformada de Hartley discreta multidimensional rápida utilizando la transformada de números de Fermat". IEE Proceedings G - Circuitos y sistemas electrónicos . 135 (6): 235–237. doi :10.1049/ip-g-1.1988.0036.

Lectura adicional

• Bracewell, Ronald N. (1986). La Transformada de Hartley (1.ª ed.). Oxford University Press . ISBN 978-0-19503969-6.
• Boussakta, Said; Holt, Alan GJ (1988). "Transformada de Hartley discreta multidimensional rápida utilizando la transformada de números de Fermat". IEE Proceedings G - Circuitos y sistemas electrónicos . 135 (6): 235–237. doi :10.1049/ip-g-1.1988.0036.
• Hong, Jonathan; Vetterli, Martin ; Duhamel, Pierre (1994). "Transformadas de campo base con la propiedad de convolución" (PDF) . Actas del IEEE . 82 (3): 400–412. doi :10.1109/5.272145.
• O'Neill, Mark A. (1988). "Más rápido que Fast Fourier". BYTE . 13 (4): 293–300.
• Olnejniczak, Kraig J.; Heydt, Gerald T. (marzo de 1994). "Escaneo de la sección especial sobre la transformada de Hartley". Actas del IEEE . 82 : 372–380.(NB. Contiene bibliografía extensa.)