CÁLCULO DE LA PROBABILIDAD DE DETECCIÓN Y FALSA ALARMA DE UN DETECTOR DIGITAL UNIFICADO

FALSE ALARM AND DETECTION PROBABILITIES CALCULATION OF A DIGITAL UNIFIED DETECTOR

 

JUAN C. VÉLEZ
Ph.D. Universidad del Norte, Profesor TC, jcvelez@uninorte.edu.co

MIJAIL Y. LISHAK
Ph.D. Moscow Power Electrical Institute. Profesor, LishakMY@mpei.ru

YURY N. ANTONOV-ANTIPOV
Ph.D. Moscow Power Electrical Institute. Profesor, AntonYN@mpei.ru

 

Recibido para revisar Enero 22 de 2009, aceptado Julio 9 de 2009, versión final Julio 27 de 2009

 

---

RESUMEN: El objetivo central del presente trabajo es encontrar las expresiones analíticas para el cálculo de la probabilidad de detección correcta y de falsa alarma de un detector digital no lineal, usado como elemento central en la etapa de detección, de un receptor que procesa señales de comunicación y navegación para bajos niveles de relación señal a ruido y altos corrimientos Doppler. Para ello se encontrará el espectro energético del cuadrado de la envolvente compleja del proceso entrante y la p.d.f resultante en la salida del procesador de Fourier. Dichos resultados permiten graficar las curvas de detección de la señal, para una probabilidad de alarma falsa dada, y comparar con resultados obtenidos experimentalmente en simulación.

PALABRAS CLAVE: detección, probabilidad de detección, probabilidad de falsa alarma.

ABSTRACT: The central objective of the present work is to find analytical expressions for the detection and false alarm probabilities computation in a nonlinear digital detector, used as the main detection device of a receiver that processes communication and navigation signals with low signal to noise ratios and high Doppler shifts. In order to accomplish this task, the energy spectrum of the squared complex envelope of the incoming process and the resulting p.d.f out of the Fourier processor are calculated. These results allow to find the curves of detection probabilities of the signal, for a given false alarm probability, and compare them with experimentally curves, found by means of simulation.

KEYWORDS: detection, detection probability, false alarm probability.

---

 

1. INTRODUCCIÓN

Entre las exigencias requeridas para los dispositivos electrónicos portátiles modernos, usados en los sistemas aeronáuticos y estaciones satelitales livianas, se encuentra la de proveer servicios integrados de comunicación, navegación e identificación [1-4]. La base para el diseño y la implementación unificada de estos sistemas son los esquemas de software-radio [5,6], gracias a su versatilidad, confiabilidad y economía. La detección de la señal y estimación de sus parámetros en la etapa inicial de entrada en servicio en este tipo de receptores unificados, que pueden operar en condiciones de muy baja relación señal a ruido, del orden de 1,5-5 dB, y con grandes rangos de indeterminación en frecuencia, 4 kHz, mejora considerablemente el desempeño de los mismos [7]. En las referencias [8,9] se presenta un detector digital que permite identificar la presencia de la señal y realizar una primera estimación gruesa de la frecuencia de la misma, el cual está basado en una transformación no lineal de la mezcla de señal y ruido, seguido luego del análisis espectral de dicha señal. Este detector es la parte básica inicial del módulo de entrada en servicio de un receptor que ha sido implementado en un sistema unificado de procesamiento de señales de comunicación y navegación satelital. El sistema opera para 4 tazas de transmisión: 2,4 Kbps, 4,8 Kbps, 9,6 Kbps y 19,2 kbps.

La idea del presente trabajo es encontrar expresiones analíticas para el cálculo de la probabilidad de detección correcta D y de falsa alarma F de dicho detector. Para ello se encontrará el espectro energético del cuadrado de la envolvente compleja del proceso de entrada, conformado por una señal BPSK y ruido blanco, así como la función de densidad de probabilidad (p.d.f.) del módulo de los valores en la salida del procesador de Fourier, para una entrada de solo ruido y señal más ruido respectivamente.

El método de análisis del proceso aleatorio en la salida del bloque no lineal se realizará utilizando el método directo [10,11]. Con la ayuda de estos resultados se graficarán las curvas de detección de la señal, para una probabilidad de alarma falsa dada, y se compararán con los valores encontrados experimentalmente por medio de simulación. Los resultados mostrados se presentarán para las señales de comunicación, debido a que las señales de navegación necesitan un procesamiento secundario que no se describe en este artículo.

El diagrama de bloques simplificado del detector se muestra en la figura 1. La mezcla inicial de la señal modulada y el ruido , da origen al proceso . En el bloque formador de la señal en cuadratura (FC) se forman dos componentes digitales de baja frecuencia en cuadratura. Estas componentes, descritas en adelante como la parte real e imaginaria de la envolvente compleja del proceso , pasan por el filtro pasa bajas digital (DLPF) y luego entran al bloque de elevación al cuadrado de la envolvente compleja , en donde se elimina la modulación de fase. Posteriormente, en el procesador de Fourier, bloque FFT, se encuentra el espectro complejo del cuadrado de la envolvente compleja .

Figura 1. Diagrama de bloques simplificado del detector digital unificado
Figure 1. Simplified block diagram of the unified digital detector

En la última etapa se busca el módulo de la componente máxima del espectro Sp_max, cuyo valor se compara con un umbral escogido previamente según el criterio de Neyman-Pearson [12] y una técnica de verificación similar a las presentadas en [13]. Si Sp_max es mayor que dicho umbral se establece una bandera que indica la presencia de señal y se forma la estimación de la frecuencia de la portadora. De lo contrario, se repite la rutina de búsqueda hasta el momento de detección de la señal.

 

2. CARACTERÍSTICAS ESTADÍSTICAS DEL PROCESO EN LA SALIDA DEL BLOQUE ELEVADOR AL CUADRADO

Para simplificar el análisis del proceso aleatorio, en la salida del bloque de elevación al cuadrado, tendremos en cuenta lo siguiente:

1. Se supondrá que el bloque FC forma dos componentes cuadráticas del proceso sin distorsión.
2. Los procesos observados serán considerados funciones de tiempo continuo t.

Escribamos la envolvente compleja de la señal BPSK como

donde es la amplitud de la señal, es la fase inicial de la señal, es la desviación de la frecuencia de la señal con respecto a su valor nominal, es la función de modulación, que toma valores equiprobables ±1.

El proceso aleatorio complejo en la salida del bloque cuadrático es igual a:

en donde y son la envolvente compleja de la señal y el ruido en la salida del bloque DLPF correspondientemente, obtenidas mediante la convolución de y la respuesta al impulso .

Encontremos la función de autocorrelación (AF) del proceso . Teniendo en cuenta que los procesos y tienen valor promedio cero, que la esperanza matemática de la envolvente compleja también es cero y utilizando los resultados presentados en [10,14], es posible demostrar que la AF del cuadrado de la envolvente compleja es

En donde y son las AF de la envolvente compleja de la señal y de su cuadrado respectivamente. De manera similar, y son las AF de la envolvente compleja del

ruido y de su cuadrado correspondientemente. De esta forma, la AF del cuadrado de la envolvente compleja del proceso tiene tres componentes: la primera de ellas ocasionada solo por la señal; la segunda es producto de la interacción entre la señal y el ruido como consecuencia de la transformación no lineal y la tercera es causada solo por el ruido.

La dependencia de la AF con respecto al tiempo es causada por los procesos transitorios en el filtro. Teniendo en cuenta que las características estocásticas del detector, se determinan a partir de la varianza del ruido en el canal del analizador de espectros FFT y que la duración de la realización del proceso , para la cual se calcula la FFT, es mucho mayor que la duración de los procesos transitorios en el filtro, entonces es posible despreciar la no estacionariedad para el cálculo de la AF del proceso . Adicionalmente, la simulación ha mostrado que, incluso, para la máxima velocidad observada (19,2 Kbps), cuando el ancho de banda de la señal es aproximadamente 1,5 veces menor que el ancho de banda del filtro, la primera componente en la expresión (3), causada por las distorsiones en el filtro, es menor en un orden de 2, que el resto de componentes. Por estas razones, en adelante no se tendrá en cuenta esta componente para el cálculo. De esta forma, la AF del proceso se encuentra aproximadamente como

Es posible demostrar que, si el ruido es Gaussiano, entonces, en correspondencia con [10]

Teniendo en cuenta esto, finalmente obtenemos:

Aplicando el teorema de Einstein-Wiener-Jinchin [10,15] a (6) y utilizando las propiedades de las transformadas de Fourier encontramos la Densidad Espectral de Potencia (PSD) del cuadrado de la envolvente compleja para la mezcla de señal más ruido:

en donde y son los espectros energéticos de la envolvente compleja de la señal y ruido respectivamente.

Veamos un caso particular, cuando el DLPF tiene una respuesta en frecuencia ideal cuadrada con un ancho de banda de paso igual al ancho de banda de la señal y no afecta la misma. En este caso, la envolvente compleja del ruido tiene un espectro energético simétrico uniforme igual a

en donde es la frecuencia de corte del filtro. El valor f_n se determina a través de la mitad del espectro de la señal y la desviación de frecuencia máxima f_smax como

El nivel del espectro energético del ruido es igual a

En donde es la varianza de las componentes en cuadratura del ruido.

La AF de la envolvente compleja de la señal BPSK es igual a

en donde es la AF de la función moduladora, la cual se define, para una duración del símbolo T_c, como [13]

Ahora, teniendo en cuenta la expresión (7), encontremos el espectro energético del cuadrado de la envolvente compleja de la mezcla de señal y ruido. Para ello denotemos

El espectro energético de la envolvente compleja de la señal es

Utilizando las expresiones (8) y (14) encontramos la convolución (13):

Haciendo el cambio de variable , llevamos la expresión (15) a la forma

y Si(x) es la función seno integral.

La integral de convolución (13), que determina , es igual a

Reemplazando en (7) las expresiones (16) y (17) obtenemos la expresión final para el espectro energético en la salida antes de la FFT

De la última expresión (18) se observa que, debido a la presencia de la señal, el espectro del proceso, en la salida del bloque encargado de elevar la señal al cuadrado, tiene una componente adicional de ruido, concentrada cerca de la frecuencia central de la señal de entrada.

En la figura 2,a se muestra los resultados del cálculo del espectro energético obtenido con la ayuda de las expresiones (16-18). Se escogieron los siguientes valores de parámetros para el cálculo: tasa de transmisión de 9,6 kbps, relación señal a ruido en la entrada del detector igual a 1,5 dB y una desviación de frecuencia de la señal de 1,875 kHz. En la figura 2,b se presenta los resultados obtenidos, teniendo en cuenta la respuesta en frecuencia total del sistema de filtrado del detector, tomada experimentalmente. Para ello se combinaron parte de los resultados obtenidos analíticamente con técnicas numéricas en la solución de la ecuación (7). Los resultados del cálculo se compararon con los resultados obtenidos en un modelo de simulación. Se observa que las curvas teóricas y experimentales prácticamente concuerdan.

Figura 2. Espectros energéticos del proceso para una tasa de transmisión de 9,6 kbps.
Figure 2. Energy spectrum of the process for a transmission rate of 9,6 kbps

 

3. CÁLCULO DE LAS CARACTERÍSTICAS ESTADÍSTICAS DE LOS PROCESOS EN LA SALIDA DEL PROCESADOR FFT

En virtud del teorema del límite central es posible considerar que cuando el tamaño de la transformada de Fourier es grande, la función de densidad de probabilidad de las muestras del espectro, en la salida del procesador de Fourier, es normal. Por esta razón la p.d.f w_s+n(z) de una muestra del espectro en amplitud en el canal del procesador de Fourier, donde cae la señal, se describe con la distribución de Rice [10]:

donde А_m и s_v² son la amplitud de la señal y la varianza del ruido en la salida del canal de la FFT, I₀ es la función de Bessel modificada de orden cero. Para encontrar la varianza del proceso en la salida del canal de la FFT que contiene la señal en un momento dado (este canal corresponde al doble de la frecuencia de la señal, puesto que ha pasado por un bloque que eleva dicha señal al cuadrado) tenemos en cuenta que en los límites de un canal del analizador de espectros la densidad espectral del ruido puede considerarse constante, por lo tanto

en donde es el ancho de banda del ruido del canal de la FFT. Es posible demostrar que este valor se determina, a través de la longitud de la realización discreta entrante N_s y la frecuencia de muestreo ó discretización , como

Encontremos la amplitud А_m de la señal como el valor promedio del proceso en la salida del canal de la FFT que contiene la señal. Este valor promedio se encuentra como la esperanza matemática del espectro del cuadrado de la envolvente compleja de la señal

en donde el operador F{ } representa la transformada de Fourier. Además

Encontremos la esperanza matemática del cuadrado de la envolvente compleja (22)

Reemplazando (24) en (22), encontramos

Veamos la expresión . Debido a que la transformada de Fourier de una multiplicación de funciones es igual a la convolución de sus transformadas de Fourier, entonces

en donde es la función de transferencia del filtro. Teniendo en cuenta (27) y (28) la expresión (26) toma la forma

en donde . Si escribimos

entonces la expresión (29) queda como

Encontremos ahora la esperanza matemática de la multiplicación de los valores de la envolvente compleja de la señal en los momentos de tiempo t₁ y t₂ , la cual entra en la expresión (25). Debido a que el espectro en amplitud no depende de la fase inicial de la señal, entonces para simplificar el tratamiento consideraremos que j_о = 0. Por consiguiente,

Reemplazando las expresiones obtenidas (30) y (31) en la fórmula para el cálculo de la media del espectro (25) y pasando de la integral en a la integral en q, obtenemos

La integral de t₁ da una función delta ubicada en el doble de la frecuencia de la señal:

. (33)

Transformemos la primera integral de (32), teniendo en cuenta la ecuación (11)

Reemplacemos en (34) la expresión definida para la función B(w,q) y cambiemos el orden de la integración:

La integral en los corchetes es la transformada de Fourier de la AF de la envolvente de la señal, es decir, el valor del espectro energético en la frecuencia . De esta manera,

La integral obtenida es la convolución de las funciones y . Teniendo en cuenta esto obtenemos finalmente la siguiente expresión para calcular la esperanza matemática del espectro complejo de la señal en la salida del bloque de elevación al cuadrado:

La amplitud de la señal А_m, para el caso de concordancia entre y la frecuencia de uno de los canales del analizador de espectros, se encuentra como

en donde N_fft es el volumen ó tamaño de la FFT. La expresión (37) se convierte en (36) cuando .

En la figura 3 se muestra la función de la relación señal a ruido en la salida de la FFT (S/N)_out = con respecto a la relación (S/N)_in en la entrada del detector para cuatro velocidades de transmisión. En este caso se escogió un desplazamiento de frecuencia de la señal de 1875 Hz, el tamaño de la secuencia de la señal fue N_s=2048 y el volumen de la FFT N_fft = 4096.

Figura 3. Función de la relación (S/N)_out en el canal de señal, en la salida de la FFT, con respecto a la relación (S/N)_in en la entrada del detector
Figure 3. (S/N)_out rate in the signal's channel, in the FFT output, as a function of (S/N)_in rate in the detector input

De dicha figura se observa que a partir de la relación S/N de entrada de 1,5 dB el valor de relación S/N en la salida es mayor a 18 dB, de esta manera se alcanza el valor relación señal a ruido exigido para una detección confiable de la señal [16]. De la gráfica también se observa que a medida que aumenta la relación S/N en la entrada del dispositivo los valores del a S/N en la salida de la FFT se unifican para todas las tasas de transmisión.

 

4. PROBABILIDAD DE ALARMA FALSA Y DE DETECCIÓN CORRECTA DE LA SEÑAL

Encontremos inicialmente la probabilidad de alarma falsa F. Para ello se debe tener en cuenta la distribución del máximo de los módulos de las componentes espectrales del ruido en N canales de la FFT y el valor del umbral de detección V_n de la siguiente manera:

en donde es el máximo de los módulos de las componentes espectrales del ruido; es la función de densidad de probabilidad acumulativa (c.d.f.) del máximo.

Teniendo en cuenta que la detección de la señal se realiza en el intervalo de desviaciones de frecuencia ±df_s, entonces el análisis debe hacerse en N canales, etiquetados con los números 1,2, ., N_dfs +1 (para las desviaciones de frecuencia positivas) y N_fft - N_dfs, N_fft - N_dfs +1, ., N_fft (para las desviaciones negativas), en donde

la operación [.] representa el redondeo de la parte entera del número.

Para encontrar la función la secuencia de muestras del espectro complejo se analiza como un proceso aleatorio discreto . Este proceso tiene una p.d.f. normal con media cero y es no estacionario, debido a que la varianza de las muestras del espectro depende del número del canal de la FFT. Para este tipo de procesos la variable aleatoria tiene una p.d.f. como la que exhibe el máximo absoluto de una secuencia de N muestras de la envolvente de un proceso Gaussiano aleatorio con envolvente compleja .

Los métodos de búsqueda de la p.d.f. de los valores máximos de la envolvente, relacionados con los procesos aleatorios estacionarios, se pueden encontrar en la monografía [17]. Debido a que en el caso general este problema no tiene una solución exacta, entonces en correspondencia con [17] se presentan dos variantes aproximadas:

1. Considerado la correlación entre dos valores, módulos de las componentes espectrales del ruido, en los canales de la FFT (aproximación independiente);
2. Sin considerar la correlación entre los valores de las componentes en 2 canales vecinos de la FFT (aproximación de primer orden).

Para el caso de la aproximación independiente, la c.d.f. de los máximos es igual a

en donde cada una de las multiplicaciones se aproxima con la expresión

En esta expresión, es la varianza de la muestra k del espectro complejo, cuyo valor se encuentra con la ayuda de la expresión (20) para el canal k, con frecuencia f_k, del procesador como

Para la aproximación de primer orden la función del máximo absoluto del proceso aleatorio estacionario, conformado por N muestras distribuidas de la misma forma se encuentra con la expresión

En donde es la c.d.f. bidimensional de los módulos de las componentes espectrales en los canales k y k+1. Gracias a que consideramos que las componentes espectrales tienen una distribución normal, entonces la función puede encontrarse como la función bidimensional de distribución cumulativa de la envolvente de un proceso aleatorio normal

en donde [18]; es la función Gamma incompleta; es el modulo del coeficiente de correlación de las componentes espectrales en los canales k y k+1 de la FFT.

El análisis de la expresión para el coeficiente de correlación cuando el volumen de la FFT N_fft es el doble de la cantidad de muestras de la señal y se cumple que N_fft >> 2p está determinado por:

Por lo tanto, el modulo del coeficiente de correlación no depende del número del canal k y aproximadamente es igual a

El valor del umbral puede ser encontrado si se resuelve la ecuación (38) teniendo en cuenta las expresiones (43) y (44) para una probabilidad de alarma falsa dada. Sin embargo, para utilizar estas formulas en la práctica es necesario calcular la varianza en N canales del procesador de Fourier.

Encontremos ahora la probabilidad de detección correcta D, para ello utilizamos la p.d.f. de la envolvente del proceso aleatorio en la salida del canal de la FFT, encontrada anteriormente (19). La probabilidad D la encontramos entonces como:

en donde

es la función Q de Marcum [16]. La varianza del ruido en la salida del canal del procesador de Fourier s_v² y la amplitud de la señal А_m se encuentran con las formulas dadas (20) y (37) correspondientemente. El valor del umbral V_n se escoge como un compromiso para garantizar una probabilidad de alarma falsa del detector.

En la figura 4 se comparan las curvas de probabilidad de detección correcta en función de la relación S/N para 4 velocidades de transmisión, obtenidas con las expresiones encontradas en este artículo y experimentalmente en un modelo simulado del detector. El tamaño de la muestra para el análisis fue de 1000 realizaciones.

Figura 4. Curvas de probabilidad de detección de la señal en función de la relación S/N (F=0,01)
Figure 4. Detection probability curves in function of the S/N relation (F=0,01)

A partir de las funciones encontradas es posible hacer las siguientes conclusiones:

• Las curvas analítica y experimental prácticamente coinciden.
• Si la relación S/N es de 1,5 dB en la entrada del detector, la probabilidad de detección D, para la señal que tiene una velocidad de 1,2 kbps es cercana a 0,99. La probabilidad D para las otras velocidades de transmisión prácticamente es igual a 1,0 bajo la probabilidad de alarma falsa dada y dentro de las restricciones planteadas en el experimento.
• El hecho de que las curvas para mayores velocidades arrojen mejores resultados, obedece principalmente a las características específicas del receptor. Dicho dispositivo permite la recepción unificada de varios tipos de señales, a varias tasas de transmisión. Para lograr esto se diseñó el filtro DLPF con un ancho de banda fijo para la máxima tasa de transmisión, teniendo en cuenta el máximo desplazamiento Doppler. Es decir, que para las otras tasas el excedente de ancho de banda relativo causa una disminución en la relación S/N.

 

5. CONCLUSIONES

En este trabajo se han encontrado las características estocásticas de un detector no lineal de comunicaciones basado en análisis espectral. Se presentan las expresiones generales para calcular la función de autocorrelación y el espectro energético en la salida del bloque cuadrático. Se obtuvo una expresión cerrada para el espectro energético para el caso cuando la respuesta en frecuencia del filtro es ideal y cuadrada. Se presentaron las curvas de relación S/N de salida en la FFT en función de la relación S/N de entrada, las cuales mostraron que para un tamaño de N_s =2048 muestras de la señal se garantiza la detección confiable de la señal. Se presentaron las expresiones para el cálculo de la probabilidad de alarma falsa y la probabilidad de detección correcta de la señal.

Para esta última, se muestran las curvas de detección. Se determinó que para una relación S/N de 1,5 dB la probabilidad de detección correcta es igual a 0,99 para una probabilidad de alarma falsa de 0,01.

 

REFERENCIAS

[1] NASA/CR-2005-213587. Technology Assessment for the Future Aeronautical Communications System. ITT Industries. Reston, Virginia, May 2005         [ Links ]
[2] WRITERS, S. The Modern Airborne Military Communications Market. Newtown CT, Sep 25, 2008.         [ Links ]
[3] CAMANA, P. INTEGRATED COMMUNICATIONS, navigation, identification avionics(ICNIA)-the next generation. Aerospace and Electronic Systems Magazine, IEEE. Volume 3, Issue 8, Aug 1988.         [ Links ]
[4] SINAY, J. Integrated communication radio navigation and identification system (ICRNI). Proceedings of the National Aerospace and Electronics Conference. NAECON. Dayton, 1990.         [ Links ]
[5] PANASIK, C. M.; NAWMAN, L. R. Programmable filter technology for integrated communication, navigation and identification systems. Proceedings of the National Aerospace and Electronics Conference. NAECON. Dayton, 1982.         [ Links ]
[6] SADIKU, M.N.O.; Akujuobi, C.M. Software-defined radio: a brief overview Potentials, IEEE Volume 23, Issue 4, Oct.-Nov. 2004.         [ Links ]
[7] ZHODZISHSKY, M.I., Mazepa, E.P: Ovsiannikov, M.I. Sistemas Receptores de radio digitales. Spravochnik. Moscú: Radio y Sviaz, 1990. (en ruso)         [ Links ]
[8] VELEZ, J.C. Características de detección de señales BPSK para diferentes tasas de transmisión. Vestnik MPEI. Moscú: MPEI Publisher, 2003. (en ruso)         [ Links ]
[9] ANTONOV-ANTIPOV, YURY N., VÉLEZ, J.C. Algoritmo de pre-procesamiento de detección de señales BPSK realizable en una FPGA. 9ª CI en ingenierías radio electrónica, electrónica y eléctrica. Moscú: MPEI Publisher, 2003. (en ruso)         [ Links ]
[10] TIKHONOV, V.I. Radioelectrónica Estadística. Moscú: Radio y Sviaz, 1982. (en ruso)         [ Links ]
[11] EVSIKOV, Y.A., Chapurskii V.V. Transformación de Procesos Aleatorios. Moscú: Vishaia Schola, 1977. (en ruso)         [ Links ]
[12] SHAKHTARIN, B.I., Detección de señales. Moscú: Bauman (MSTU) Publisher, 2005. (en ruso)         [ Links ]
[13] LINATTI J.H. On the Threshold Setting Principles in Code Acquisition of DS-SS Signals. IEEE Journal On Selected Areas In Communications, Vol. 18, No. 1, January 2000.         [ Links ]
[14] DAVENPORT, W., Root, W. An introduction to the theory of random signals and Noise. New York:,McGraw-Hill, 1958.         [ Links ]
[15] HAYKIN, S. Communication Systems. Third edition. New York: John Wiley & Sons, INC, 1994.         [ Links ]
[16] BORISHOV, V.I., Inmunidad de los Sistemas Radioelectrónicos de banda ancha. Moscú: Radio y Sviaz, 2000. (en ruso)         [ Links ]
[17] FOMIN, Y.A., Teoría de los extremos de procesos aleatorios. Moscú: Radio y Sviaz, 1980. (en ruso)         [ Links ]
[18] ABRAMOVITZ, M., Stigan, I. Manual de Funciones Especiales con Formulas, Gráficas y Tablas. Moscú: Nauka, 1979. (en ruso).         [ Links ]