• Author / Uploaded
• santiago

• Categories
• Transformada rápida de Fourier
• Transformada discreta de Fourier
• Algoritmos
• Matriz (Matemáticas)
• Áreas de la informática

Citation preview

UNIVERSIDAD NACIONAL DE PIURA ESCUELA PROFESIONAL DE INGENIERIA ELECTRONICA Y TELECOMUNICACIONES

16-7-2015

TRANSFORMADA RAPIDA DE FOURIER

PROFESOR: ING.MIGUEL ANGEL PANDURO CURSO: PROCESAMIENTO DIGITAL DE SEÑALES I ALUMNO: JHONATHAN SANCHEZ RAMIREZ

PROCESAMIENTO DIGITAL DE SEÑALES I

INDICE

1 RESUMEN 2 FUNDAMENTOS 2.1 ANALISIS DE FOURIER EN TIEMPO DISCRETO 2.2 TRANSFORMADA DISCRETA DE FOURIER (DFT) 2.3 EJEMPLO DE LA DFT 3 DEDUCCIÓN DEL ALGORITMO FFT DE COOLEY – TUKEY 3.1 EJEMPLO DEL ALGORITMO FFT 4 FFT EN MATLAB 4.1 4.2 4.3 4.4

SINTAXIS DEFINICIONES DESCRIPCIÓN EJEMPLO

5 BIBLIOGRAFIA

PROCESAMIENTO DIGITAL DE SEÑALES

TRANSFORMADA RÁPIDA DE FOURIER (FFT) 1. RESUMEN La transformada rápida de Fourier (Fast Fourier Transform FFT) es un algoritmo extremadamente rápido para calcular la transformada discreta de Fourier (Discrete Fourier Transform DFT), ambos son métodos que en la práctica ejecuta un computador sobre datos discretos. A menudo cuando se presentan señales en el tiempo de larga duración, se hace inviable ejecutar la DFT, por esto fue necesaria la creación de la FFT, originalmente descubierta por Gauss (1805), fue redescubierta por J. W. Cooley y O. W. Tukey en IBM durante 1960, C.S. Burrus, de la Universidad de Rice University siendo jefe del departamento de Ingeniería, literalmente "escribió el libro" de los algoritmos de la rápida Transformada Discreta de Fourier DFT. Existen distintos métodos para calcular la FFT, en general los podemos dividir en 2 tipos: decimación en el tiempo, y diezmado en la frecuencia. El algoritmo busca reducir el número de multiplicaciones efectuadas en la DFT, reduciendo el número de cálculos para N datos de 2N2 a 2N·log2N, donde N debe ser una potencia de 2.Usualmente la presentación del algoritmo FFT se realiza de forma polinomial pero también puede ser presentado de forma matricial. La FFT explota las simetrías en la matriz W (ec.2.7) para aproximarse a una matriz diagonal. En la actualidad existen algoritmos aun más eficientes de calcular la DFT que inclusive el algoritmo FFT de Cooley y Tukey.

2. FUNDAMENTO: 2.1 ANALISÍS DE FOURIER EN TIEMPO DISCRETO: Una señal discreta x[n] será periódica si se cumple: x[n]  x[n  N ] en donde N, periodo fundamental, es un entero mínimo. La exponencial compleja e ejemplo de función periódica discreta.

 j 2N n

es un

El análisis de Fourier en tiempo discreto es similar a su análogo en tiempo continuo, pero una de las grandes diferencias que presenta en que las series ahora no presentaran infinitos términos sino que estarán determinados por el número del periodo N. Una señal periódica puede representarse en términos de exponenciales complejas de la forma:

x [n ] 

N2



k  N1

ak e

jk 2N n

con N2  ( N1 )  N

(2.1)

Esta es la representación de una serie de Fourier de una señal discreta periódica; para hallar el k-ésimo coeficiente ak multipliquemos por e

e

 jr 2N n

x [n ] 

N2



k  N1

TRANSFORMADA RAPIDA DE FOURIER

ak e

 jr 2N n

ambos miembros en (2.1):

jr 2N n

e

 jk 2N n

-1-

PROCESAMIENTO DIGITAL DE SEÑALES

Puesto que x[n] es periódica da lo mismo que n  [ N1, N2

o n  [0, N . Ahora

tomando sumatoria para 0  n  N : N 1

e n 0

 jr 2N n

N 1 N 1

x[n ]   ak e

jr 2N n  jk 2N n

e

n 0 k 0

N 1

N 1

k 0

n 0

  ak  e

j ( r k ) 2N n

(*)

Pero veamos que:

Entonces en (*): N 1

e

 jr 2N n

n 0

N 1

N 1

k 0

n 0

x[n ]   ak  e

j ( r  k ) 2N n

 ar N

Luego:

ak 

1 N 1  jk 2N n  e x [n ] N n 0

(2.2)

Esta última se llama ecuación de análisis, es aplicable solo a una función periódica para obtener su serie discreta de Fourier (SDF). Veamos ahora que en analogía a la variable continua nuestros resultados se pueden extender para hallar la SDF de señales de duración finita como se ve en la figura:

x[n] (a)

-N1

N2

n

xp[n] (b)

-N1

N2

n

FIGURA 1:

TRANSFORMADA RAPIDA DE FOURIER

-2-

PROCESAMIENTO DIGITAL DE SEÑALES

Ahora, sea x[n] una señal aperiódica de duración N podemos construir una señal periódica xp[n] de periodo N tal que:

 x [n ] ; N1  n  N2 x [n ]   p 0 ; N2  n y n  N1

(**)

Entonces podemos hallar la representación de la SDF de xp[n] sobre N1  n  N2 en donde se debe cumplir que:

ak 

1 N



N 1 n 0

e

 jk 2N n

x p [n ]

Ahora para que xp[n] se acerque mas a x[n] podemos hacer que el periodo sea más grande, es decir que en la figura.1.b los ciclos de xp[n] estarán cada vez más alejados y como x[n]  0 n fuera de  N1  n  N2 , podemos escribir:

1   jk 2N n ak   e x [n ] N n  Definamos la función:

(2.3)

X (e j )   n x[n]e  jn ,

Entonces en (2.3) vemos que:

ak  N1 X (e jok ) , con o  2N De lo cual llegamos a:

X (e jok ) 



 x[n]e

 jko n

(2.4)

n 

Esta expresión se conoce como la transformada de Fourier en tiempo discreto. A partir de (2.1) se observa que x[n] se puede expresar también como:

x [n ]  Como X (e luego:

jo k

1  X (e jo k )e jko n  N k 

(2.5)

) es periódica podemos coger el intervalo de la sumatoria de 0 a N-1,

x [n ] 

1 N 1 X (e jo k )e jko n  N k 0

N 1

X (e jo k )   x[n]e  jko n

TRANSFORMADA RAPIDA DE FOURIER

(2.6)

n 0

-3-

PROCESAMIENTO DIGITAL DE SEÑALES

2.2 TRANSFORMADA DISCRETA DE FOURIER (DFT): Con el resultado de (2.6) podemos intentar calcular la transformada para un conjunto j k de N datos, por simplicidad hagamos X (e o )  X (k )

n x[n] Y e

jo k

e

 jk 2N n

0 1

1 2

2 3

3 4

 w Nkn , desarrollando: X (0)  X (1) 

 

N 1 n 0

e

N 1 n 0

 j (0) 2N n

x [n ]  (1  1 

 j (1) 2N n

x[n ]  (w 1·0  w 1·1  N N

e

X (N  1)   Nn 01e

 j ( N 1) 2N n

 1) N 1)  w 1·( ) N

x[n ]  (w (NN 1)0  w (NN 1)1 

 w (NN 1)( N 1) )

Puede ser expresado de forma matricial como:

1  X (0)  1  (1)(1)   1 w N  X (1)          ( N 1)(1)  X (N  1)  1 w N

  x [0]    wN   x [1]      ( N 1)( N 1)  wN   x [N  1]  1

(1)( N 1)

(2.7)

Equivalente a: X  W  x donde W es denominada matriz de Fourier. Un hecho muy importante y evidente es que W es una matriz simétrica. Ahora veamos algunos casos prácticos en donde se usa la DFT:

2.3 EJEMPLOS DE LA DFT a.- Ejecutemos la DFT para el siguiente caso (N=4) por el método ordinario: Calculando los valores de wN  e kn

n

 ( j )kn :

0

1

2

3

e0·0=1 e1·0=1 e2·0=1 e3·0=1

e0·1=1 -j -1 j

e0·2=1 -1 1 -1

e0·3=1 j -1 -j

k

0 1 2 3

 j ( 24 )kn

Entonces de la formula (2.7):

1 1 1 1   1   X (0)   10  1  j 1 j   2   X (1)   2  2 j         1 1 1 1  3   X (2)   2          1 j 1  j   4   X (3)   2  2 j 

TRANSFORMADA RAPIDA DE FOURIER

(2.8.1)

-4-

PROCESAMIENTO DIGITAL DE SEÑALES

b.- Veamos ahora el caso para n=8: n x[n]

0 1

1 2

2 3

3 4

4 5

5 6

7 8

6 7

Procediendo como en el caso anterior encontramos ahora que: 1 1  2 2 1 2  j 2 1 j   1  22  j 22 W  1 1 1  2  j 2 2 2  j 1  2 2 1 2  j 2

1

1

j



2 2

1

j

1

2 2

j 2 2

j

2 2

1

j

2 2

2 2

j 

j

2 2

2 2

j

1

2 2

j

2 2

j

1

j

2 2

2 2

2 2

1

1

2 2

j

1 

2 2

j

  2 j  j 22  2  1 j   j  22  j 22   1 1  2 2  j  2 j 2  1 j   2 j  j 22  2 1

j

1

j

1

1  1

j

1

1

2 2

1

(2.8.2)

Entonces ejecutando W x  X , obtenemos: k X(k)

0 36

1 -4 -j9.65

2 -4 -j4

3 -4 -j1.65

4 -4

5 -4+j1.65

6 -4 +j4

7 -4+j9.65

(*) OBSERVACIÓN: Veamos que en la ec. (2.8.1) pudimos haber realizado la siguiente factorización: Observemos que: w 4 Entonces: X (k ) 



nk

3 n 0

e

 j 24 nk

 ( j )nk

x[n ] w 4 nk

Desarrollando:

X (k )  x[0]  x[1](  j )k  x[2](  j )2k  x[3](  j )3k  x[0]  x[1]( j )k  x[2]( 1)k  x[3]( j )k Reordenando: X (k )  ( x[0]  x[2]( 1)k )  ( x[1]  x[3]( 1)k )(  j )k Y sea ahora:

(*.1)

 ( x[2·0]  x[2·1]( 1) )  ( x[2·0  1]  x[2·1  1]( 1) )(  j )k xpar [r ]  x[2r ] , xipm [r ]  x[2r 1] ; r  0,1 k

k

Y además como el valor de (-1)k solo depende de si k es par o impar, podemos hacer que:

X par (r )  ( x par [0]  x par [1]( 1)r ) , y X ipm (r )  ( xipm [0]  xipm [1]( 1)r ) .

TRANSFORMADA RAPIDA DE FOURIER

-5-

PROCESAMIENTO DIGITAL DE SEÑALES

Entonces para visualizar esto mejor desarrollemos (*.1) para este caso:

X (0)  ( x[0]  x [2](1) )  ( x [1]  x [3](1)

)(1)

X (1)  ( x [0]  x [2](1))  ( x [1]  x [3](1) )( j )

(*.2)

X (2)  ( x[0]  x [2](1) )  ( x [1]  x [3](1) )(1) X (3)  ( x[0]  x [2]( 1) )  ( x [1]  x [3](1))( j )

En donde x[0]=1, x[1]=2, x[2]=3, x[3]=4, entonces X par (0)  4 ; X par (1)  2

y Xipm (0)  6 ; Xipm (2) , entonces: X (0)  X par (0)  X ipm (0)(1)  10 X (1)  X par (1)  X ipm (1)(  j )  2  2 j

(*.3)

X (2)  X par (0)  X ipm (0)( 1)  2 X (3)  X par (1)  X ipm (1)(1)  2  2 j

Que es el mismo resultado que obtuvimos anteriormente, pero ahora se necesitaron efectuar menos multiplicaciones. La ec. (*.3) nos invita a hacer el siguiente diagrama de desarrollo:

x[0]

Xpar(0)

X(0)

x[2]

Xpar(1)

X(1)

x[1]

Ximp(0)

X(2)

x[3]

Ximp(1) FIGURA 2

X(3)

En la ecuación (2.7) y de los ejemplos a y b notamos que el número de multiplicaciones complejas a ejecutar es N2 (sin considerar el hecho que w N0  1 ). El método matricial que hemos visto en la práctica es ejecutado por un computador; sin embargo cuando N es muy grande como en señales de video y audio digitales N  106 y se debe realizar 1012 multiplicaciones aproximadamente, entonces los cálculos resultaran muy tediosos, ocuparían mucha memoria y además tomaran mucho tiempo. Esto último fue la motivación para desarrollar un método más directo para calcular la DFT, el algoritmo de la transformada rápida de Fourier FFT, el cual acabamos de introducir implícitamente en la observación anterior.

TRANSFORMADA RAPIDA DE FOURIER

-6-

PROCESAMIENTO DIGITAL DE SEÑALES

3. DEDUCCIÓN DEL ALGORITMO FFT DE COOLEY – TUKEY: Partimos de la DFT para una señal dada polinomialmente como: N 1

X (k )   x[n] w N k ·n

(3.1)

n 0

Ahora asumiendo que el número de datos N es par, descomponemos la sumatoria en sus términos pares e impares: N / 2 1

X (k ) 





x[2n] w N k ·2n  w N k

n 0 N / 2 1



n 0

Sea ahora: x[2n]  f [n] y x[2n  1]  g [n] notemos que podemos hacer que: k ·2n

wN

N / 2 1

x[2n] w N k ·2n 

e

 j 2N k ·2n

 n 0

x[2n  1] w N k ·(2n 1)

N / 2 1

 n 0

x[2n  1] w N k ·2n n  [0, N / 2  1] ,

donde el nuevo

e

 j N2/ 2 k ·n

 wN / 2kn

Y definamos a:

F (k ) 

N / 2 1

 n 0

G(k ) 

f [n] w Nk ·/n2

N / 2 1

 g [n ] w n 0

k ·n N/2

(3.2)

Además veamos que F (k) y G (k) son también periódicas:

F (k  N / 2) 

N / 21

 n 0

f [n] w N( k/2N / 2)·n 

N / 21

 f [n ] w n 0

k ·n N/2

e

 j N2 / 2 ( N / 2)·n

 F (k )

Análogamente para G (k); con N/2 como periodo mínimo. Entonces nuestro problema de calcular una DFT para N datos se redujo ahora a calcular 2 DFT para N/2 datos cada una. Nuevamente asumamos que N/2 es par y ahora, para F (k):

F (k ) 

N / 2 1

 f [n ] w n 0



N / 21

 f [n ] w n 0

k ·n N/2





k ·n N/2

N / 4 1

 f [2n] w n 0

N / 4 1

 f [n ] w n 0

f

k ·2 n N/2

k ·2 n N/2

w



N / 4 1

 f [2n  1] w n 0

N / 4 1 k N/2

 g [n ] w n 0

f

k ·2 n N/2

k ·(2 n 1) N/2

n  [0, N / 4  1]

Donde ff , gf son lo mismo para f[n] que para x[n]. De igual forma para g[n]:

G(k ) 

N / 4 1

 g [2n] w n 0

k ·2 n N/2



N / 4 1

 g [2n  1] w n 0

TRANSFORMADA RAPIDA DE FOURIER

k ·(2 n 1) N/2

-7-

PROCESAMIENTO DIGITAL DE SEÑALES



N / 4 1

 n 0

fg [n] w Nk ·2/ 2n  w Nk / 2

N / 4 1

 g [n ] w f

n 0

k ·2 n N/2

n  [0, N / 4  1]

Es decir que ahora sobre f[n] y g[n] se realizan 4 DFT de longitud N/4. Entonces podríamos hacer múltiples divisiones del intervalo [0,N-1] mientras se pueda dividir N entre 2. Ahora estamos listos para generalizar el método, entonces: Sea x un vector de datos, de longitud N =2m. Entonces sobre el intervalo [0,N-1] se pueden realizar m particiones como se mostró anteriormente hasta llegar a una DFT de longitud 2 ,esta es la unidad básica del FFT conocida como mariposa (o butterfly en ingles) en donde solo se necesitara una multiplicación y 2 sumas complejas: como se muestra en la figura:

Figura 3: Los elementos computacionales básicos de la transformada rápida de Fourier es la mariposa. Toma dos números complejos, representados por a y b, y forma las cantidades mostradas. Cada mariposa requiere una multiplicación compleja y dos sumas complejas.

3.1 EJEMPLO DEL ALGORITMO FFT Para especificar la idea realicemos nuevamente el ejemplo b, para N=23 por el algoritmo FFT: Desarrollando:

X (k )   7n 0 x[n] w 8 k ·n

(3.1.1)

 x[0]w80k  x[1]w8k  x[2]w82k  x[3]w83k  x[4]w84k  x[5]w85k  x[6]w86k  x[7]w87k Donde se cumple para w:

w 8 2n·k  e

 j 28 2 n·k

w 8(2n 1)·k  e

e

 j 24 n·k

 j 28 (2 n 1)·k

e

 w 4n·k

 j 24 n·k

e

 j 28 k

 w 4 n·k w 8k

Agrupando términos pares con impares:

X (k )  ( x[0]w 40k  x[2]w 4k  x[4]w 42k  x[6]w 43k ) ( x[1]w 40k  x[3]w 4k  x[5]w 42k  x[7]w 43k )w8k

TRANSFORMADA RAPIDA DE FOURIER

-8-

PROCESAMIENTO DIGITAL DE SEÑALES

Tomamos: x[2r ]  f [r ] y x[2r  1]  g [r ] , r  0,1,2,3 luego:

X (k )  (f [0]w 40k  f [1]w 4k  f [2]w 42k  f [3]w 43k ) (g[0]w 40k  g[1]w 4k  g[2]w 42k  g[3]w 43k )w8k

(3.1.2)

Pero dentro de cada paréntesis de (3.1.2) podemos realizar una nueva factorización, sabiendo que:

w 4 2 n·k  e

 j 24 2 n ·k

w 4(2 n 1)·k  e Entonces:



e

 j 22 n ·k

 j 24 (2 n 1)·k

e

 w 2 n·k

 j 22 n ·k

e

X (k )  (f [0]  f [2]w 4 2k )  (f [1]  f [3]w 4 2k )w 4 k



 j 24 k

 w 2n·k w 4 k





 (g [0]  g [2]w 4 )  (g [1]  g [3]w 4 2k )w 4 k w 8 k 2k

Nuevamente: f [2s]  ff [s] y f [2s  1]  fg [s ] , s  0,1

Luego:



X (k )  (ff [0]  ff [1]w 2k )  (fg [0]  fg [1]w 2k )w 4 k







 (gf [0]  gf [1]w 2 )  (g g [0]  g g [1]w 2k )w 4 k w 8 k k

Remplazando los índices iniciales:



X (k )  ( x [0]  x[4]e  j k )  ( x[2]  x[6]e  j k )e



 j 2 k

 ( x [1]  x [5]e  j k )  ( x[3]  x[7]e  j k )e

=>



 j 2 k

X (k )  ( x[0]  x[4]( 1)k )  ( x[2]  x[6]( 1)k )(  j )k



 

e

(***)

 j 4 k



 ( x[1]  x[5](1)k )  ( x[3]  x[7](1)k )(  j )k e

 j 4 k

(3.1.3)

Esta es la forma general de se obtiene k-ésimo termino para N=8, además vemos que se deberán ejecutar 4 DFT de longitud 2, y los resultados que se obtengan se usaran para ejecutar 2 DFT de longitud 4, para lo cual se realizo 3 particiones. Para interpretar lo que indica la ec. (3.1.3) veamos la figura 4, donde cada flecha en diagonal representa una suma y las acompañan sus factores multiplicativos:

TRANSFORMADA RAPIDA DE FOURIER

-9-

PROCESAMIENTO DIGITAL DE SEÑALES

FIGURA 4

El número total de cálculos que se realizara serán N=8 sumas para cada etapa y log2N=3 etapas, haciendo el número de procesos básicos de (8)(3)=N log2N. En general se cumple que para la FFT el número de cálculos es: N log2N. Comparando con la DFT de N2 cálculos. Para cuantificar la diferencia veamos la siguiente tabla para distintos valores de N: N N2 Nlog2N

2

4

8

16

32

64

128

256

512

1024

4 2

16 8

64 24

256 64

1024 160

4096 384

16384 896

65536 2048

262144 4608

1048576 10240

El número de cálculos es directamente proporcional al tiempo de solución:

FIGURA 5: Esta figura muestra que tan lento crece el tiempo de solución de un proceso de Nlog2N.

TRANSFORMADA RAPIDA DE FOURIER

- 10 -

PROCESAMIENTO DIGITAL DE SEÑALES

Esta diferencia en la rapidez de cálculos de la FFT + computación digital fue completamente responsable de la "explosión" del Procesamiento Digital de Señales DSP en los años 60's.

4. FFT EN MATLAB Para terminar hay que mencionar que el algoritmo FFT se encuentra implementado en muchos programas de computación. Veamos las instrucciones que se usan en MATLAB.

4.1 Sintaxis Y = fft (x) Y = fft(X, n) Y = fft(X, [], dim) Y = fft(X, n, dim)

4.2 Definiciones Las funciones Y = fft(x), y = ifft(X) implementan el par transformar y transformada inversa dado para vectores de longitud N por:

Dónde:

Es un N º raíz de la unidad.

4.3 Descripción Y = fft(x) Y = fft(x) devuelve el discreto Transformada de Fourier (DFT) de vector x, calculado con una transformación (FFT) algoritmo de Fourier rápida. Si la entrada X es una matriz, Y = fft(X) devuelve la transformada de Fourier de cada columna de la matriz. Si la entrada X es una matriz multidimensional, fft opera en la primera dimensión más de un elemento. Y= fft(X, n), Y = fft(X, n) devuelve el n -point DFT. Fft(X) es equivalente a fft(X, n) donde n es el tamaño de X en la primera dimensión más de un elemento. Si la longitud de X es menor que n, X se rellena con ceros a la longitud n. Si la longitud de X es mayor que n, la secuencia X se trunca. Cuando X es una matriz, la longitud de las columnas se ajustan de la misma manera.

TRANSFORMADA RAPIDA DE FOURIER

- 11 -

PROCESAMIENTO DIGITAL DE SEÑALES

Y = fft(X, [], dim) Y = fft(X, [], dim), Y = fft(X, n, dim) se aplica la operación FFT través de la dimensión dim.

4.4 EJEMPLO Ejecutar el ejemplo a y b:

A.-Para n=22, x= (1, 2, 3,4): >> x=1:1:4 x = 1 2 3 >> y=fft(x)' y = 10.0000 -2.0000 - 2.0000i -2.0000 -2.0000 + 2.0000i

4

B.-Para n=23, x= (1, 2, 3, 4, 5, 6, 7,8) >> x=1:1:8 x = 1

2

3

4

5

6

7

8

>> y=fft(x)' y = 36.0000 -4.0000 -4.0000 -4.0000 -4.0000 -4.0000 -4.0000 -4.0000

- 9.6569i - 4.0000i - 1.6569i + 1.6569i + 4.0000i + 9.6569i

TRANSFORMADA RAPIDA DE FOURIER

- 12 -

PROCESAMIENTO DIGITAL DE SEÑALES

5. BIBLIOGRAFIA  OPPENHEIM, Alan V., WILLSKY, Alan S. y NAWAB, S. Hamid. Señales y Sistemas. 2da. Edición, Prentice Hall, 1998.  ROBERTS, Michael J. Señales y Sistemas: Análisis mediante métodos de transformada y MATLAB. 1ra. Edición, Mc Graw Hill, 2004.  PROAKIS, John G. y MANOLAKIS, Dimitris G. Tratamiento Digital de Señales. 4ta Edición. Prentice Hall, 2007.

Referencias en Internet:  http://www.mathworks.com/help/matlab/ref/fft.html  https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm

TRANSFORMADA RAPIDA DE FOURIER

- 13 -