Fasor y Fourier

CAPITULO 4

La Transformada de Fourier

4.1. INTRODUCCIÓN El estudio de las señales cotidianas en el dominio de la frecuencia nos proporciona un conocimiento de las características frecuenciales de éstas. Por ejemplo nos es muy útil el conocer la respuesta en frecuencias de un canal de telecomunicaciones, para poder determinar la máxima frecuencia que puede transmitir sin provocar distorsiones de la señal de modo que ésta sea recibida y reconstruida con total garantía. Otro ejemplo puede ser el análisis de la voz. Transformando la señal de voz en sus componentes frecuenciales podemos distinguir las diferencias entre las voces de distintas personas y determinar las palabras que está diciendo. Esto es muy útil para reconocimiento e identificación de voz, dos aplicaciones que están haciéndose muy familiares, al tiempo que aumenta la velocidad de los DSP’s. También podemos poner como ejemplo los Analizadores de Espectro que realizan transformación de señales en el dominio del tiempo al dominio de la frecuencia. Este tipo de instrumentos se han hecho muy populares en los laboratorios electrónicos. Se utilizan para examinar la señal de salida de sensores conectados a estructuras que soportan esfuerzos mecánicos, como puentes o rascacielos en los cuales un cambio significativo de la respuesta en frecuencia puede suponer un exceso de carga en alguna parte de la estructura que puede provocar un derrumbamiento en el futuro. En este capítulo intentamos presentar cómo se transforman las señales de un dominio al otro. Suponemos que no se tienen conocimientos previos en la utilización de este tipo de transformadas, de modo que el estudio va a ser introductorio, ya que las matemáticas implicadas en este campo pueden llegar a ser muy complejas, y hay cantidad de libros de texto al respecto. Vamos a ver en este capítulo cómo el advenimiento de los sofisticados DSP’s que hay en la actualidad ha hecho posible el uso de estas, relativamente viejas, transformadas. Digo esto porque Jean Baptiste Joseph Fourier derivó su fórmula en el año 1822. En la parte final de este capítulo presentaremos también otro tipo de transformada en frecuencia que se está haciendo muy popular en aplicaciones de procesado de imágenes. Desafortunadamente, es imposible realizar una explicación práctica de las transformadas en frecuencia sin hacer uso de unas pocas matemáticas.

Página

72

CAPITULO 4


4.2. EL FASOR Como punto de inicio necesitamos un método simple de describir una señal, que usualmente llamaremos fasor (phasor model). En la figura 4.1 presentamos un fasor simple.

Figura 4.1 - Modelo de un fasor simple

El fasor es un vector rotando en el plano complejo, con una magnitud A y una velocidad de rotación ω expresada en radianes/segundo. Si nosotros cogemos un valor instantáneo de la señal en cualquier momento, x(t), viene dado por : x(t)=(coordenada real) + j(coordenada imaginaria) = a + jb donde : A = a 2 + b 2

θ

= ω ⋅ t = tan −1

b a

Éste es un método para escribir un numero complejo (recuerda : j2=-1), llamado la forma rectangular. El otro método es la forma polar donde :

Página

73

CAPITULO 4

La Transformada de Fourier x (t ) = Ae j (ω ⋅t ) e j (ω ⋅t ) = cos(ω ⋅ t ) + j sen(ω ⋅ t )

Estas ecuaciones básicas dan lugar a muchas de las funciones de la trigonometría que la mayoría de nosotros hemos tenido que aprender en la escuela. Recuerda también que ω se relaciona con la frecuencia de la siguiente forma : ω

= 2π f

Finalmente π radianes equivale a 180o . Nuestra descripción del fasor se puede extender fácilmente al tiempo discreto o sistemas digitales, en los cuales la señal sólo tiene valor en los intervalos específicos de tiempo definidos por el periodo de muestreo Ts: x (n) = Ae j ( n ⋅ω ⋅T s )

A diferencia de la variable continua de tiempo t, nosotros tenemos ahora una variable discreta n, de modo que el fasor avanza en saltos de Ts. Cogiendo cualquiera de los dos casos (continuo o discreto), si nosotros tenemos un valor inicial de x de : x ( 0) = Ae j (α )

nosotros podemos derivar la forma general de ambas ecuaciones como sigue : x (t ) = Ae j (ω ⋅t +α )

ó x (n) = Ae j ( n ⋅ω ⋅T s +α ) Estas simples ecuaciones son la base de todo el análisis que a continuación haremos y podremos ver cómo usamos este modelo para describir señales que nos son muy familiares : senos y cosenos.

Página

74

CAPITULO 4


4.3. MODELANDO SENOIDES Volviendo atrás en la descripción de e j ⋅

ω ( t )

, nosotros la podemos reescribir así :

e j ⋅θ = cosθ + j sen θ

También : e − j ⋅θ = cosθ − j sen θ

Donde : θ

= (ω ⋅ t + α ) ó (n ⋅ ω ⋅ T s + α )

De estas dos ecuaciones nosotros podemos derivar las siguientes relaciones:

cos θ =

e j ⋅θ + e − j ⋅θ

2 ⋅θ

sen θ =

e j − e

− j ⋅θ

2 j

Esto significa que una señal seno o coseno general, x(t), se puede describir con la suma de dos fasores. Por ejemplo : x ( t ) = R cos(ω ⋅ t + α )

=

R

2

( e j (

ω ⋅ t +α )

+ e− j (

ω ⋅t +α )

)

Por tanto nuestro coseno puede representarse con dos fasores que forman un número complejo conjugado. Esto significa que tiene el mismo valor real (a) e igual y el valor de b igual y de signo contrario, ver figura 4.2. Debemos recordar que nuestros fasores están rotando tanto en dirección positiva (sentido horario), como negativa (sentido antihorario), podemos deducir una propiedad interesante que es que las señales reales están hechas de un par de fasores que forman un número complejo conjugado de manera que la suma de los vectores siempre está sobre el eje real.

Página

75

CAPITULO 4


Figura 4.2 - Representación de una señal coseno mediante un fasor

Página

76

CAPITULO 4


4.4. SERIES DE FOURIER El desarrollo anterior lo hemos realizado sobre señales simples del tipo seno o coseno, de una sola frecuencia, pero existe la forma de dividir señales mucho más complejas en distintos senos o cosenos. Por ejemplo un tren de pulso rectangulares consiste en la suma de un número infinito de senos de distinta amplitud. De este modo podemos describir cualquier señal periódica compleja en forma de una suma de muchos fasores (senos o cosenos). Un método de describir una señal de esta forma es las llamadas series de Fourier, en las cuales asumimos que el conjunto de fasores tienen frecuencias múltiplos de alguna frecuencia fundamental, f 0 (o frecuencia angular ω0) : N

x ( t ) =

∑C e (

k = − N

j k ⋅ω 0 ⋅t )

k

Cualquier señal periódica se puede representar como una serie de Fourier donde N tiende a infinito. Los componentes frecuenciales individuales son conocidos como armónicos. Podemos hacer el modelo de Fourier más general usando fasores cuyas frecuencias no son múltiplos del armónico fundamental, lo que ocurre generalmente cuando la señal no es periódica (la mayor parte de los casos en aplicaciones reales) : N

x ( t ) =

∑C e ( j

k = − N

ω k ⋅ t )

k

Cualquier señal arbitraria se puede representar como una serie de Fourier de esta forma general.

4.4.1 SERIES DE FOURIER DISCRETAS Necesitamos traducir estas ecuaciones de tiempo continuo al dominio discreto o digital para poder derivar algunas fórmulas útiles para los DSP’s. El análisis anterior se puede extender a sistemas de tiempo discreto. Lo único que necesitamos es reemplazar la función continua, t, con otra que varíe en saltos de ω0 ×Ts, así para el caso en que la señal es periódica : N

x ( n) =

∑C e (

j kω 0 T s n )

k = − N

k

Es interesante advertir que cuando el incremento de la fase para el k-ésimo armónico está dado por : kω0T s = 2π ⋅ m

Página

77

CAPITULO 4


donde m es entero, la fase no es distinguible de la fase cuando k=0. Esto ocurre porque 2π=360o y ocurre cuando :

k

ω0 ω s

= 1, 2,3,4,...

Si razonamos un poco lo anterior, podemos ver que esto significa que la respuesta en frecuencia de una señal discreta es periódica, con un periodo de 1/Ts. Usando el modelo del fasor simple hemos conseguido describir una señal discreta. Haciendo uso de esta descripción vamos a explicar cómo podemos realizar la transformación entre los dominios del tiempo y la frecuencia.

4.4.2 LA TRANSFORMADA DE FOURIER. En aplicaciones reales la mayor parte de las señales no son periódicas y debemos transformar nuestras series de Fourier para poder analizarlas. Consideremos la serie de Fourier general en la cual todas las frecuencias son múltiplos de la fundamental, es decir : ωk

= k ω0

El hecho de que la señal no es periódica se puede representar por : ω0

→0

Esta simple ecuación expresa que no hay mínimo común denominador entre las frecuencias de todos lo fasores. Cuando el número de fasores tiende a infinito nuestro sumatorio se convierte en una integral :

1 ∞ j ( x ( t ) = X (ω ) ⋅ e ∫ 2π − ∞

ω ⋅ t )

⋅ d ω

En la anterior ecuación asumimos que la amplitud de la señal se puede definir como función de la frecuencia (ω), es decir, X(ω). La ecuación inversa que define X(ω), viene dada por : ∞

X (ω ) =

∫ x (t ) ⋅ e (

j − ω ⋅ t )

⋅ dt

−∞

Por lo tanto, ahora tenemos una ecuación que nos permite calcular la respuesta en amplitud de una señal continua en el dominio de la frecuencia usando su respuesta en el dominio del tiempo. Estas dos ecuaciones se denominan el par de Transformadas de Fourier, las cuales son muy útiles para los matemáticos pero, Página

78

CAPITULO 4


desafortunadamente no es posible implementarlas directamente en un DSP. Por ello debemos deducir la forma discreta de estas ecuaciones. Obviamente, necesitamos un número infinito de medidas para determinar la función X(ω). En la práctica las transformadas de Fourier no se calculan cada vez, nosotros simplemente usamos tablas de transformadas de Fourier que aparecen en muchos libros de matemáticas o libros de DSP’s.

4.4.3 LA TRANSFORMADA DISCRETA DE FOURIER (DFT). Para poder encontrar el equivalente discreto de la transformada de Fourier debemos transformar la variable continua t por la variable discreta nTs. Fuera de los límites ±π/Ts, el espectro se repite, por lo tanto podríamos cambiar los límites de este valor. Por conveniencia, sin embargo, haremos que la variable de integración sea ωTs así la integral se convierte en : π

1 ( x ( n) = X (ω )e j ∫ 2π −

ωT s n )

⋅ d (ωT s )

π

La transformada inversa es : ∞

X (ω ) =

∑ x ( n )e (

j − ωT s n )

n = −∞

Hay que advertir que la segunda ecuación aún usa un sumatorio en lugar de una integral. Podríamos haber esperado esto porque sabemos que x(n) sólo es válido en los instantes de tiempo nTs . Estas dos ecuaciones forman la DFT (Discrete Fourier Transform), la cual vamos a ver con un poco más de detalle. Estas dos ecuaciones nos permiten transformar señales digitales entre los dominios del tiempo y la frecuencia. El espectro que obtenemos al usar la DFT tiene algunas propiedades interesantes, por ejemplo, es periódico con frecuencia ωs. También para una señal real, como los fasores aparecen en forma de complejos conjugados, el espectro siempre tiene simetría par a lo largo del eje real, y simetría impar en el eje imaginario. Esto simplemente significa que si nosotros sabemos que trabajamos con una señal real, la cantidad de información que necesitamos recordar sobre el espectro en frecuencias es menor, ya que es repetitivo. Haciendo uso de estas ecuaciones podemos calcular la respuesta en frecuencia de una señal de voz digitalizada, o determinar los coeficientes de un filtro digital.

Página

79

CAPITULO 4


4.4.4 CONSIDERACIONES PRÁCTICAS Hasta aquí hemos introducido un método para describir una señal que varía en el tiempo, al cual hemos denominado el modelo del fasor. Hemos hecho la explicación en función de señales senoidales y cosenoidales, y hemos introducido las series de Fourier que nos indican cómo una señal periódica se puede representar por medio de un número determinado de senos (armónicos) equiespaciados en frecuencia. Después de hacer las series de Fourier más generales eliminando las relaciones entre los armónicos, hemos aplicado la forma discreta de nuestra ecuación, la cual finalmente nos ha proporcionado el par de fórmulas de la DFT. De este modo tenemos ahora un par de ecuaciones que nos permiten transformar cualquier señal discreta entre los dominios del tiempo y la frecuencia, pero ahora nos hacemos la siguiente pregunta, ¿cómo podemos programar un DSP para que realice la DFT ? ¿qué problemas podemos encontrar cuando intentemos diseñar un algoritmo que realice la DFT de una señal ?. Volviendo atrás a la ecuación de la DFT para X(ω), observamos dos problemas obvios : primero, es que en el mundo real no es posible realizar un sumatorio infinito. En segundo lugar, nosotros siempre vamos a estar limitados en el tiempo requerido para computar la salida. De este modo el número de frecuencias al que nosotros podemos aplicar la DFT se ve restringido por estas dos causas. El primer problema es fácil de superar, debemos coger sólo una parte de los valores de entrada x(n). Esto usualmente se denomina ‘enventanado’, y también se usa para muchas otras aplicaciones como por ejemplo codificación de voz. Si entonces realizamos la DFT sobre los valores de la ventana, el espectro resultante viene dado por : X N(ω)=X(ω) * W(ω) donde X N(ω) indica el espectro enventanado, N es el número de muestras usadas, W(ω) denota el espectro de la ventana y el * es la convolución entre W(ω) y X(ω). La convolución es una operación que ya conocemos y que podemos expresar de la siguiente manera : N − 1

X N (ω ) =

∑0 X (r ) *W ( N − r ) r =

Si recordamos que el espectro en frecuencia será periódico, con frecuencia ωs, podemos ver que un espectro rectangular en el dominio de la frecuencia implica que no causaremos ninguna interferencia entre lóbulos adyacentes de la respuesta. Desafortunadamente una respuesta en frecuencia rectangular es prácticamente imposible de realizar, por tanto siempre deberemos llegar a un compromiso.

Página

80

CAPITULO 4


Asumamos que hemos escogido una función de enventanado, lo que necesitamos conocer ahora es el efecto de usar sólo un número limitado de frecuencias, y cuántas frecuencias necesitamos para mantener una precisión aceptable. Básicamente no hay una respuesta simple a esta cuestión. Por regla general el número óptimo de fasores es igual al número de puntos iniciales de x(n), esto es, N. La forma más simple de imaginar esto es asumir que nuestra porción enventanada de x(n) consiste en un periodo de una larga secuencia (ver figura 4.3), con periodo NTs y frecuencia ωs/N. Si hacemos esto podemos tratar la secuencia como una serie de Fourier y veremos cómo el espectro consiste de N fasores.

Figura 4.3 - Porción enventanada de la señal como un periodo de la secuencia

En muchos libros detallados sobre DSP’s probablemente no encontraremos el par de ecuaciones de la DFT exactamente igual como aquí lo hemos explicado. Generalmente usan otra variable en la ecuación a la que se denomina factor W (en inglés twiddle factor). Retrocedamos un poco y veamos cómo aparece esta nueva función. Dijimos antes que el espectro de una DFT es repetitivo con periodo ωs, así si decimos que los fasores están espaciados por un valor δ (figura 4.4), así : Nδ = ωs

Página

81

CAPITULO 4


Figura 4.4 - Fasores espaciados por δ

Esto nos permite digitalizar nuestra escala de frecuencias de modo que el espectro puede escribirse ahora en términos de k en lugar de ω : N − 1

X (kδ ) = K (δ ) =

∑0 x (n ) ⋅ e (

j − kδT s n )

n=

Recordando que :

ω s

=

δ

=

2π T s

y que : ωs

N

Podemos reescribir la ecuación de este modo : N − 1

X N (k ) =

∑0

x (n ) ⋅ e

− j (

2πkn N

)

n=

Esta es la forma de la DFT práctica que normalmente se usa y el “factor twiddle”, denominado W N, se define como : W N = e

− j (

2π N

)

Página

82

CAPITULO 4


De modo que finalmente podemos escribir el par de ecuaciones de la DFT en su forma más común : N − 1

X N ( k ) =

∑0 x (n) ⋅ W

kn N

n=

y:

x ( n) =

1 N − 1

∑X N 0 k =

N

( k ) ⋅ W N − kn

4.4.5 LA TRANSFORMADA RÁPIDA DE FOURIER (FFT) La DFT es un método muy eficiente para determinar el espectro en frecuencia de una señal, el único inconveniente de esta técnica es la cantidad de tiempo requerida para computar la salida. Esto se debe a que los índices k y n, deben variar de 0 a N-1, para conseguir el rango de salida completo, y por tanto se debe realizar N2 operaciones. Si miramos las anteriores ecuaciones podemos observar que cada parte de la ecuación realiza una función de multiplicación-y-suma, así que una DFT con N=1000 (DFT de 1000 puntos), usa 106 ciclos máquina en cualquier DSP. Si usamos un DSP con un tiempo de ciclo de 50 ns le costará 0.05 segundos, y por lo tanto la mayor velocidad de muestreo que podrá alcanzar será de 20 Hz, lo cual evidentemente es muy poco teniendo en cuenta que cualquier señal audible puede alcanzar frecuencias del orden de 20 KHz o algo más. Si examinamos de nuevo el ‘factor twiddle’, W N : W N = e

− j (

2π N

)

vemos que calculamos muchas veces los mismos valores de W N durante la DFT debido a que W N es una función periódica con un número limitado de valores distintos. El objetivo de la transformada rápida de Fourier (FFT) y su inversa, la IFFT, es usar esta redundancia para reducir el número de cálculos necesarios. El término genérico “transformada rápida de Fourier” abarca distintos algoritmos con distintas características, ventajas y desventajas. Por ejemplo, una FFT diseñada y optimizada usando un lenguaje de alto nivel probablemente no funcionará bien en un DSP de coma fija. Sin embargo, todas las FFT’s usan la misma aproximación para reducir el algoritmo en un número reducido de DFT’s sucesivas, cortas y simples. Vamos a ver en detalle un método para reducir la DFT en una FFT. Aunque esto requiere una pequeña cantidad de

Página

83

CAPITULO 4


matemáticas, nos permitirá encontrar la explicación de por qué casi todos los DSP’s están diseñados para implementar las FFT’s tan rápido como sea posible. El primer paso en esta reducción es dividir la señal de entrada x(n) en muchas secuencias cortas intercaladas. Esto, usualmente, se llama ‘diezmo en el tiempo’. Si cogemos nuestra señal original con N valores, en primer lugar nos permite dividirla en dos secuencias, una con los números impares y otra con los número pares : N − 1

X N ( k ) =

∑0 x (n) ⋅ W

kn N

n=

se transforma en : N

2

N

−1

∑0 x (2r ) ⋅ W 2

X N ( k ) =

2

rk

N

r =

−1

+ ∑ x (2r + 1) ⋅ W N ( 2 r +1) k r = 0

o: N

X N (k ) =

2

N

−1

∑0 x (2r ) ⋅ (W 2 ) N

r =

2

rk

−1

+ W N k ∑ x (2r + 1) ⋅ (W N 2 )rk r = 0

Volviendo atrás en la definición de W N : W N = e

j ( −

2π N

)

Así : 2

W N = e

j ( −

2π N

2)

=e

j (−

2π N / 2

)

= W N / 2

Y nuestra ecuación se convierte en : N

2

X N ( k )

N

−1

∑0 x (2r ) ⋅ W r =

rk N / 2

2

−1

+ W N k ∑ x (2r + 1) ⋅ W N rk / 2 r = 0

Lo cual escribimos : Página

84

CAPITULO 4

La Transformada de Fourier X N ( k ) = G ( k ) + WN k H (k )

donde G(k) es la DFT de los puntos pares y H(k) es la DFT de los números impares. Desafortunadamente debemos multiplicar la DFT impar por W Nk antes de poderle sumar G(k). De este modo tenemos expresada la DFT original en forma de dos DFT’s más pequeñas de longitud N/2. Si N=1000, el cálculo de esto nos llevará 5002+5002+500=500500 en lugar de 106 de la DFT original. Si asumimos que la longitud de la transformada es una potencia de 2, podemos seguir realizando este diezmado y dividir las dos secuencias anteriores en dos partes cada una. El límite de este proceso nos lleva a simples DFT’s de dos puntos, una de las formas de la FFT más usada, la de base 2. Como ejemplo cogeremos una secuencia pequeña donde N=8, es decir, x(n) viene definida por : n={0,1,2,3,4,5,6,7} Si diezmamos una vez obtenemos dos secuencias : n={0,2,4,6,} y {1,3,5,7} Si diezmamos una vez más obtenemos cuatro secuencias : n={0,4} {2,6} {1,5} {3,7} Entonces la FFT se calcula inicialmente realizando cuatro DFT’s de dos puntos sobre los valores de los pares anteriores. Aunque esto simplifica tremendamente la operación, debemos recordar los factores ‘twiddle’ que aparecen cada vez que diezmamos la secuencia. Estos factores son también diferentes en cada etapa del diezmado. El objetivo de un buen algoritmo FFT es ser capaz de incorporar estos factores extra sin que estorbemos demasiado al DSP. Veamos un simple ejemplo donde N=4 para observar las operaciones a realizar : 3

X 4 ( k ) =

∑0 x (n) ⋅ W 4

kn

n=

1

= ∑ x ( 2r ) ⋅ W2 + W4 rk

k

r = 0

1

∑0 x (2r + 1) ⋅W 2

rk

r =

= [ x ( 0) + x (2 ) ⋅ W2 ] + W4k [x (1) + x (3) ⋅ W 2k ] k

Podemos reescribir W N, de modo que :

Página

85

CAPITULO 4


W2 = e k

j (−

2π k) 2

=e

j( −

2π 2k ) 4

= W42 k

Entonces nuestra FFT se transforma en : X 4 ( k ) = [ x (0) + x (2 ) ⋅ W42 k ] + W4k [x (1) + x (3) ⋅ W 42 k ]

Así desarrollando para todos los valores de k, obtenemos : X 4 ( 0) = [ x (0) + x (2 ) ⋅ W40 ] + W40 [ x (1) + x (3) ⋅ W 40 ] X 4 (1) = [ x (0) + x (2 ) ⋅ W42 ] + W41[x (1) + x (3) ⋅ W 42 ] X 4 ( 2) = [ x (0) + x (2 ) ⋅ W40 ] + W42 [ x (1) + x (3) ⋅ W 40 ] X 4 ( 0) = [ x (0) + x (2 ) ⋅ W42 ] + W43 [ x (1) + x (3) ⋅ W 42 ]

Observar que : W44

=e

j ( −

2π )4 4

= 1 = W 40

y que : W46

=e

j ( −

2π )6 4

= −1 = W 42

Es más sencillo visualizar esta DFT de cuatro puntos si hacemos uso de un gráfico como el de la figura 4.5 donde los números de los círculos representan la potencia de W4 que se necesita en cada etapa. Si no hay símbolo, entonces no se necesita multiplicar. El gráfico ilustra una de las características principales de cualquier tipo de FFT, la mariposa :

Página

86

CAPITULO 4


Figura 4.5 - Mariposa característica de una FFT

cada sección está compuesta de varias mariposas como la mostrada en la figura 4.6.

Figura 4.6 - Diagrama de una mariposa básica

La forma generalmente aceptada para la mariposa es en la que los factores de multiplicación son x=0 e y=-1. Esto se toma como convención de modo que muchos libros asumen esto y dibujan la mariposa tal y como se muestra en la figura 4.7. En el caso de que tengamos factores ‘twiddle’ extra dibujaremos la mariposa tal y como aparece en la figura 4.8. Ahora la mariposa consiste solamente en simples sumas y restas, ya que todos los factores

Página

87

CAPITULO 4


‘twiddle’ son +1 ó -1. Cuando hemos conseguido reducir las ecuaciones a este simple formato se vuelve extremadamente sencillo el diseñar un programa que realice la FFT.

Figura 4.7

- Diagrama generalmente usado para representar una mariposa

Figura

4.8 - Mariposa básica con factores “twiddle” extra

Un DSP es ideal para realizar FFT’s ya que puede ejecutar las matemáticas implicadas en éstas muy rápidamente. Los DSP’s tienen también características importantes que vamos a ver a continuación. En cualquier FFT, con el propósito de obtener los valores de salida en orden ascendiente, debemos desordenar la secuencia de entrada de una forma específica. Veremos que la consecuencia de no desordenar la secuencia de entrada provocará que la secuencia de salida salga desordenada. Si miramos cómo se puede implementar una FFT en un DSP observaremos que son necesarias dos tablas de datos, una con los valores de entrada y otra con los valores de salida. Con una FFT de cuatro puntos, para obtener la salida en orden ascendiente los valores de entrada tienen que cargarse en la tabla en el siguiente orden : n={0,2,1,3} Esto es trivial para una FFT de cuatro puntos, pero si nosotros tenemos que realizarlo para una FFT de 512 puntos o 1024, o incluso más, se convierte en una tarea compleja. La mayor parte de los DSP’s de propósito

Página

88

CAPITULO 4


general actualmente ofrecen una solución a este problema que es lo que se denomina direccionamiento con inversión de bit (bit-reversed adressing). Si vamos a hacer uso de este modo de direccionamiento, los valores iniciales se almacenan en orden ascendente, y con este modo de direccionamiento se van cogiendo de la forma desordenada que necesitamos. A continuación vamos a ver cómo trabaja este modo de direccionamiento : Consideremos una dirección de 3 bits, la cual representa 23 u 8 posibles posiciones de memoria. Por lo tanto estamos delante de una FFT de ocho puntos. Después de empezar en 0, lo que vamos a hacer es sumar la mitad de la longitud de la FFT a cada dirección de memoria, pero sumando el acarreo de izquierda a derecha, en lugar de la forma habitual de derecha a izquierda :

inicio en 000

= x(0)

000 + 100

= 100 = x(4)

100 + 100

= 010 = x(2)

010 + 100

= 110 = x(6)

110 + 100

= 001 = x(1)

001 + 100

= 101 = x(5)

101 + 100

= 011 = x(3)

011 + 100

= 111 = x(7)

En los dispositivos DSP más modernos esta manipulación de las direcciones se hace como una tarea en segundo plano de modo que nos podemos beneficiar del uso de las tablas sin ninguna penalización temporal al utilizar direccionamiento con inversión de bit. Esto es posible en dispositivos que tienen unidades aritméticas asociadas a los generadores de direcciones. Son estas ALU’s especializadas las que realizan las operaciones mostradas anteriormente. Usando estos modernos DSP’s se ha hecho posible implementar FFT’s de 1024 puntos con tiempos por debajo de los 5 ms, y esto es debido a que el tiempo para realizar los cálculos se ha reducido a un valor tan pequeño que actualmente podremos realizar tareas en tiempo real. Hay incluso procesadores cuya única finalidad es realizar FFT´s, obviamente estos dispositivos son los más rápidos en el cálculo de las mismas.

Página

89

CAPITULO 4


Volviendo atrás a nuestro ejemplo de la FFT de 4 puntos, veamos los valores de W que obtenemos en las mariposas : W41

=e

j ( −

= e − j 90

2π ) 4

o

= cos( −90 ) + j sen( −90 ) o

o

=− j y: W42 = (W41 ) 2 = j 2 = −1

finalmente : W43 = W41 W42 = j

Por tanto, podemos simplificar el gráfico tal y como lo representamos en la figura 4.9. Como hemos apuntado anteriormente, multiplicar por -1 es simplemente realizar una resta, y multiplicar por ± j simplemente es una conversión entre números reales e imaginarios.

Página

90

CAPITULO 4


Figura 4.9 - Gráfico simplificado de una DFT de 4 puntos

Por lo tanto la DFT es muy simple y es posible convertir el diagrama anterior a su forma estándar que representamos en la figura 4.10.

Figura 4.10 - Gráfico estándar de una DFT de 4 puntos

Página

91

CAPITULO 4


Hasta aquí hemos reducido las matemáticas tanto como ha sido posible, hasta el punto de tener sólo un factor W con valor -j. Ésta es sólo una forma de describir la FFT, desafortunadamente hay muchos libros en los cuales se explica lo mismo de diferente forma, lo cual puede crear confusiones en un principio. Sin embargo, es muy posible que con la explicación realizada, no haya problemas para comprender todos los métodos que se explican en otros textos. Tal y como dijimos al principio de la explicación, la ventaja de descomponer la DFT original en DFT’s más pequeñas es la reducción del número de operaciones necesarias para conseguir la función final. Si nosotros descomponemos una DFT de N puntos a DFT’s de dos puntos conseguimos reducir el número de operaciones que realiza el procesador desde N2 hasta Nlog2 N. Si cogemos como ejemplo una DFT de 512 puntos para ilustrar lo que acabamos de decir, podemos ver que se reduce el número de operaciones desde 260000 hasta 4600 consiguiendo una reducción por encima de 50 veces, aproximadamente. Recordemos que el término transformada rápida de Fourier es general y lo usamos para muchos métodos que reducen el número de operaciones requerido para realizar DFT’s. Nosotros hemos visto aquí las FFT’s de base 2 que son las más populares pero no las únicas. Fácilmente podemos deducir del ejemplo anterior cómo la descomposición en DFT’s de dos puntos puede fácilmente alcanzarse con FFT’s de cualquier tamaño mientras sea potencia entera de dos. Obviamente si éste no es el caso deberemos buscar un método distinto.

Página

92

CAPITULO 4


4.5. IMPLEMENTACIÓN HARDWARE DE LA FFT En términos generales, la memoria necesaria para ejecutar una FFT de N puntos consiste en 2N posiciones para datos (N partes reales y N partes imaginarias) y N posiciones para los valores senoidales (coeficientes W o twiddle factors). Una vez satisfechos los requerimientos de memoria, el DSP debe ejecutar los cálculos necesarios en el tiempo requerido. Muchos fabricantes de DSP’s ofrecen bien un benchmark para una FFT de un tamaño concreto, o bien el tiempo que precisa el cálculo de una mariposa. Al comparar especificaciones, es importante asegurarse de que se ha utilizado el mismo tipo de FFT en todos los casos. Por ejemplo, el benchmark de una FFT de 1024 puntos pudiera haber sido derivado de una FFT de base 2 o de base 4; ambos casos no son comparables porque la cantidad de cálculos necesarios es diferente. Una vez el hardware cumple los requerimientos, el software ha de hacer que el sistema funcione. Con el mismo hardware, diferentes rutinas harán posible ejecutar algoritmos distintos base 2, base 4, diezmado en el tiempo, o diezmado en la frecuencia, los cuales manejan los mismos datos de diferentes maneras.

Página

93

CAPITULO 4


4.6. CONSIDERACIONES PARA ALGORITMOS DE FFT El primer paso al diseñar una FFT consiste en determinar el número de puntos requeridos, N, o longitud del registro. Hay varias formas de abordar este problema. La frecuencia de muestreo, f s, debe ser al menos dos veces la máxima frecuencia de interés en la señal de entrada (por el teorema del muestreo de Nyquist). Una vez fijada la velocidad de muestreo, la resolución espectral viene dada por f s/N. Cuantos más puntos tenga la FFT, mejor será dicha resolución. Esta es una consideración de principio en aplicaciones de análisis espectral. En análisis del habla en tiempo real, por ejemplo, el ancho de banda de la señal es de unos 4kHz, implicando una frecuencia de muestreo mínimo de 8kHz. El espectro del habla no es estacionario (en realidad, está continuamente cambiando). Por lo tanto, para estudiar la señal hay que dividirla en ventanas de tiempo, Tw, suficientemente cortas como para que los matices individuales no se vean promediados en la FFT ; todos los significados se pierden en FFT’s para tratamiento del habla de largo término. Pero Tw debe ser suficientemente largo como para dar la resolución espectral adecuada. Está comprobado, que para el habla humana, 20ms es un tiempo adecuado. Por tanto, haremos Tw igual a 20ms. Veamos un ejemplo de aplicación de una FFT para el análisis de voz en tiempo real : - BW = 4kHz. ⇒ Frecuencia de muestreo = 8kHz. - Ventana = 20ms. - N > 8kHz × 20ms = 160 ⇒ N = 256 puntos. - ¿ Puede cumplir el DSP estas condiciones ? En general, un DSP “antiguo” a 12MHz de velocidad de reloj es capaz de realizar esta FFT en 0,59ms de modo que le quedan 19,41ms para realizar otras tareas. Ahora, ¿qué determina si un DSP sirve para esa tarea ?. El número de puntos de nuestra ventana Tw es igual a Tw× f s, es decir, 20ms × 8KHz = 160 puntos, lo que, redondeado a la potencia de 2 inmediata superior, nos da 256 puntos. Esto significa que el DSP deberá completar la FFT de 256 puntos en menos del tiempo de adquisición de datos por ventana, Tw. De lo contrario no será posible el proceso en tiempo real. Los tiempos de prueba de cálculo de FFT’s suelen ser suministrados por los fabricantes de DSP’s. Al evaluar diferentes procesadores hay que asegurarse de que todos los datos se consiguieron bajo las mismas condiciones. Por ejemplo, una FFT de base 4 es algo más rápida que una de base 2.

Página

94

CAPITULO 4


4.7. FUGADO ESPECTRAL. El fugado espectral durante el proceso de una FFT puede ser entendido mediante la consideración del caso en que la entrada es una onda senoidal pura. Consideraremos dos situaciones. En la primera, reflejada en la figura 4.11 en su apartado a , la razón entre la velocidad de muestreo y la frecuencia de la señal senoidal de entrada es tal que dentro de la ventana de datos se contiene un número exacto de ciclos de la señal de entrada. Esto provoca que la respuesta espectral nos de un único armónico correspondiente a la señal senoidal de entrada. La figura 4.11 apartado b , muestra la segunda situación, en la que la ventana de datos no contiene un número entero exacto de ciclos de la señal de entrada. Las discontinuidades en los extremos son equivalentes a multiplicar la onda por un pulso rectangular “de ventana”, cuyo espectro tiene la forma sen(x)/x. Las discontinuidades en el dominio del tiempo dan así lugar al fenómeno llamado de fugado espectral (spectral leakage) en el dominio de la frecuencia, debido a que se precisa gran cantidad de términos espectrales para describir la discontinuidad. La respuesta espectral obtenida mediante FFT’s se aleja de la teórica consistente en un solo armónico (figura a), para parecerse más a la de la figura b, con el lóbulo central horizontalmente dilatado y un gran número de lóbulos laterales cuyas características básicas proceden del pulso rectangular en el tiempo.

Figura

4.11 - Ejemplo de aparición del efecto de fugado espectral

Puesto que en las aplicaciones práctica de análisis espectral por medio de la FFT las frecuencias exactas son desconocidas, debe hacerse algo para minimizar los lóbulos laterales. Este algo consiste en elegir una función diferente de la rectangular para la ventana : las muestras de entrada en el tiempo son multiplicadas por una función tal que reduce la señal a cero en los extremos de la ventana. La selección de la señal de ventana apropiada es habitualmente un equilibrio entre el ensanchamiento del lóbulo principal y la minimización de los laterales. El fugado puede ser reducido también ampliando los datos con ceros y ejecutando una FFT por lo tanto más larga. Página

95

CAPITULO 4


4.8. TRANSFORMADA DISCRETA DE COSENO (DCT). En este apartado vamos a ver por encima otro tipo de transformada en frecuencia a la que denominamos la Transformada Discreta de Coseno (DCT). El interés por esta transformada ha crecido en los últimos años gracias a que se ha hecho posible implementarla en los modernos DSP’s para aplicaciones en tiempo real. La DCT es una de un conjunto de técnicas heurísticas de codificación de imágenes, que comúnmente sustituyen algoritmos óptimos pero que consumen grandes cantidades de tiempo en la realización de los cálculos. Estas técnicas no encuentran la solución exacta, pero sí una muy aproximada que en la mayoría de los casos es válida y consumen mucho menos tiempo de computación. Una DCT puede proporcionar tasas de compresión del orden de 20:1, lo cual significa reducir una imagen de 8 bits con 1024×1024 puntos a un fichero de tan sólo 50k bytes. Esto se debe a que la DCT transforma la imagen en sus componentes frecuenciales usando una simple operación, mientras que en comparación, por ejemplo, en codificación de banda baja, es necesario un bloque de filtros pasa banda en paralelo. Invariablemente, para realizar la DCT, necesitamos partir la imagen en bloques de 8×8 puntos antes de poder aplicarla. Una razón de esto es que necesitamos la DCT para encontrar la redundancia existente en un conjunto de puntos. Básicamente cualquier punto en una imagen es muy probable que esté estrechamente relacionado con los cuatro puntos que le rodean y similarmente cada uno de estos cuatro puntos es muy probable que tenga alguna relación con los cuatro puntos que tenga más cercanos, pero el punto original es poco probable que esté relacionado con otro punto que esté muy lejano de él. Por tanto dividiendo la imagen podemos tener la esperanza de formar grupos de puntos que estén estadísticamente relacionados consiguiendo así un alto grado de redundancia. La segunda razón de dividir la imagen en bloques más pequeños es reducir la cantidad de cálculos necesarios. La DCT es similar a la DFT, además no es práctico aplicar la DCT sobre toda la imagen sin dividir, ya que los cálculos llevarían mucho tiempo, y no se podría realizar procesamiento de imágenes en tiempo real. Después de aplicar la DCT en un bloque de puntos obtenemos un conjunto de datos que están relacionados frecuencialmente. Si cogemos una imagen de televisión o una fotografía, generalmente encontramos que las componentes de baja frecuencia (incluyendo la componente de continua) son más grandes que las componentes de alta frecuencia. Entonces podemos ignorar los coeficientes de alta frecuencia , o podemos usar un sistema de codificación de longitud variable. Este tipo de sistemas de codificación de longitud variable lo que hacen es asignar más bits a aquellos coeficientes que contienen la información crucial, esto es, los de baja frecuencia en nuestro caso, mientras que para el resto se usa un número reducido de bits. De esta forma podemos

Página

96

CAPITULO 4


reducir el número de bits necesario para transmitir la imagen y por tanto podemos reducir el ancho de banda de la señal a transmitir. Reducir el ancho de banda de la señal es el fin fundamental que se busca usando técnicas de compresión de imágenes. Si volvemos atrás en el ejemplo anterior, la imagen de 8 bits y 1024×1024 puntos requiere alrededor de 50k bytes después de que se haya aplicado la DFT sobre ella. Teniendo en cuenta esto y que el número de imágenes por segundo que por ejemplo utiliza un vídeo-teléfono está alrededor de las 25, necesitamos un canal capaz de transmitir la señal resultante a una velocidad de 1.25 Mbytes/seg, ¡ aún asumiendo que no vamos a consumir ningún tiempo en decodificar la imagen en su destino !. Si dividimos la imagen en bloques más pequeños, en primer lugar conseguimos que la cantidad de información sea menor, después de que la DCT haya sido aplicada, de modo que el tiempo necesario para aplicar la DCT inversa (IDCT) será también menor. En segundo lugar se ve incrementada la relación entre los puntos que queremos reducir con lo cual se reduce el número de coeficientes a transmitir y conseguimos así reducir la velocidad necesaria para transmitir la señal resultante. Vamos a explicar brevemente las matemáticas envueltas en la DCT. La DCT de una serie de puntos discretos x(n)=0,1,2,3,4,...(N-1), se representa por la siguiente ecuación :

X c ( k ) =

2

2

N

N − 1

⋅ C ( k ) ⋅ ∑ x (n ) ⋅ cos n=0

( 2n + 1) k π 2 N

Para k=0,1,2,3,4,..,(N-1). Donde : C ( k ) =

2

1 2

C ( k ) = 1

y,

para k=0

para k=1,2,3,...,(N-1)

De la misma forma, la IDCT viene dada por :

x (n) =

2

N − 1

 (2n + 1)k π  ⋅ ∑ C (k ) ⋅ X c (k ) ⋅ cos  N k = 0  2 N 

2

En muchas aplicaciones donde no podemos tolerar ningún error de los que se introducen por el hecho de dividir la imagen en bloques más pequeños necesitamos usar un método basado en las DFT’s. Un ejemplo de este caso son las imágenes médicas, en las que la señal recibida es muy débil siendo la relación señal-ruido muy pobre de modo que algunas características son muy difíciles de distinguir y no nos podemos permitir el perder detalles.

Página

97

CAPITULO 4


Sin embargo en vídeo-teléfonos los efectos de dividir la imagen en bloques más pequeños no son distinguibles por el ojo humano y las DCT’s están siendo ampliamente utilizadas en este tipo de aplicaciones, aunque no se pueden encontrar directamente ya que ellas forman parte de algoritmos de compresión de imágenes más extensos y complejos. En los últimos años han aparecido transformadas rápidas de coseno (FCT), pero actualmente aún están en fase de desarrollo e investigación.

Página

98

Fasor y Fourier

Recommend Documents