Fast Fourier Transform (FFT)
Yarko Ni
no C. y Paulo Herrera R.
Departamento de Ingeniera Civil, Universidad de Chile
Semestre Primavera 2011
MODELACION NUMERICA CON APLICACIONES EN INGENIERIA HIDRAULICA Y AMBIENTAL
Transformada de Fourier
Un proceso fsico puede ser descrito en el dominio del tiempo como una
funcion h(t) o en el dominio de frecuencias como H(f ). La relacion
entre las dos funciones est
a dada por,
Z
H(f ) = h(t)e2if t dt
Z
h(t) = H(f )e2if t df
Si h(t) es funcion del tiempo, entonces f tiene unidades de ciclos por
segundo. Tambien puede ser que h es funci on de la posicion, entonces
f es la inversa de la longitud de onda y tiene unidades de ciclos por
metro. Otras unidades tambien son posibles.
2 / 26
Transformada de Fourier (cont.)
En fsica es com
un expresar la frecuencia como frecuencia angular
[radianes/segundo]. En ese caso,
w = 2f H(w) = [H(f )]f =w/2
y el par de transformadas se escribe como,
Z
H(w) = h(t)eiwt dt
Z
1
h(t) = H(w)eiwt dw
2
Es decir, el uso de H(w) rompe la simetra del par de
transformadas. Por eso es mejor usar H(f ) para problemas
computacionales.
3 / 26
Propiedades de transformada de Fourier
h(t) es real H(f ) = [H(f )]
h(t) es imaginaria H(f ) = [H(f )]
h(t) = h(t) (par) H(f ) = H(f )
h(t) = h(t) (impar) H(f ) = H(f )
h(t) es real y par H(f ) es real y par
h(t) es real e impar H(f ) es imaginaria e impar
h(t) es imaginaria y par H(f ) es imaginaria y par
h(t) es imaginaria e impar H(f ) es real e impar
4 / 26
Propiedades de transformada de Fourier
Transformada de Fourier es una transformacion lineal.
1 f
h(at) |a| H a time scaling
1 t
|b| h b H(bf ) frequency scaling
h(t t0 ) H(f) e2if t0 time shifting
h(t)e2if0 t H(f f0 ) frequency shifting
gh G(f )H(f ) Convolution Theorem
Corr(g, h) G(f )H(f )) Correlation Theorem
R
gh g( )h(t ) d
R
Corr(g, h) g( + t)h( ) d
5 / 26
Propiedades de transformada de Fourier (cont.)
Teorema de Wiener-Kinchin (autocorrelaci
on):
Corr(g, g) |G(f )|2
donde usamos G(f ) = [G(f )] si g(t) es real.
Theorema de Parseval (Energa Total):
Z Z
2
E |h(t)| dt = |H(f )|2 df
6 / 26
Power Spectral Density (PSD)
PREGUNTA: Cuanta energa est
a contenida entre f y f + df ?
Si la funcion h(t) es real entonces,
Ph (f ) |H(f )|2 + |H(f )|2 = 2|H(f )|2 0f <
Ph es llamada one-sided power spectral density (PSD).
Normalmente, h(t) es evaluada o medida en un intervalo acotado de
tiempo t. En ese caso, se define la PSD por unidad de tiempo,
dividiendola por t.
7 / 26
Ejemplo: Topografa.
8 / 26
Ejemplo: Topografa.
9 / 26
Ejemplo: Topografa.
10 / 26
Series discretas de valores
En la mayora de los casos h(t) es medida o evaluada a intervalos
regulares. Si el espaciamiento de las muestras es , entonces los valores
evaluados son:
hn = h (n) n = . . . , 2, 1, 0, 2, 3, . . .
y la frecuencia de muestreo es f = 1/.
Teorema de Nyquist:
Sea,
1
fc
2
Si la funcion h(t) contiene un espectro acotado de frecuencias, tal que
H(f ) = 0, |f | fc , entonces la funci
on h(t) es completamente
determinada por los valores medidos hn .
11 / 26
Teorema de Nyquist
Si la funcion h(t) contiene un espectro acotado de frecuencias en el
intervalo fc < f < fc , entonces pueder ser reconstruida a partir de los
valores medidos con intervalo , usando,
X sin[2fc (t n)]
h(t) = hn
n=
(t n)
NOTAS:
I Este teorema establece que la informaci
on contenida en una
funcion continua que tiene un espectro de frecuencia limitado,
puede ser representada por un numero finito de valores.
I Tambien es conocido como Nyquist-Shanon Sampling
Theorem en Teora de la Informaci
on.
12 / 26
Teorema de Nyquist (cont.)
Uno de los corolarios mas importantes del Teorema de Nyquist es que:
Si la funcion h(t) contiene componentes de frecuencia m as alta que
fc , la energa de esos componentes es movida al intervalo de
frecuencias fc < f < fc (aliasing).
En general, aliasing no se puede remover.
Sin embargo podemos chequear si el contenido de la PSD se acerca a
cero a media que las frecuencia se acerca a fc o fc . Si no lo hace, es
probable que aliasing es importante.
13 / 26
Aliasing
14 / 26
Transformada de Fourier Discreta
Considerando N valores consecutivos,
hk h(tk ) = h(k) k = 0, 1, 2, . . . , N 1
Supuestos:
I h(t) es no cero solo en los N puntos muestreados. Si no,
I h(t) es peri
odica
I N es par (por ahora).
Tenemos N valores podemos producir N resultados. Definiendo,
n N N
fn , n= ,...,
N 2 2
I N + 1 valores en total, pero N valores independientes.
I Frecuencias extremas corresponden a las frecuencias de Nyquist.
15 / 26
Transformada de Fourier Discreta (cont.)
Aproximando la transformada de Fourier como una sumatoria,
Z N
X 1 N
X 1
2ifn t 2ifn tk
H(fn ) = h(t)e dt hk e = hk e2ikn/N
k=0 k=0
La transformada de Fourier discreta de los hk valores es definida como,
N
X 1
Hn hk e2ikn/N
k=0
NOTA: Hn no depende de ning un par
ametro dimensional.
La transformada de Fourier continua puede ser calculada a partir de la
transformada discreta usando,
H(fn ) Hn
16 / 26
Transformada de Fourier Discreta (cont.)
Z N
X 1 N
X 1
H(fn ) = h(t)e2ifn t hk e2ifn tk = hk e2ikn/N
k=0 k=0
Hn es periodica en n con perodo N , entonces,
Hn = HN n n = 1, 2, . . .
Por lo tanto, mas facil usar n = 0, 1, . . . , N 1. Con este sistema:
I f = 0 corresponde a n = 0.
I fc > f > 0 corresponde a 1 n N/2 1
I fc < f < 0 corresponde a N/2 + 1 n N 1
I fc y fc corresponde a N/2
17 / 26
Transformada de Fourier Discreta (cont.)
Todas las propiedades de la transformada de Fourier tambien son
alidas para Hn , con los cambios H(f ) Hn y H(f ) HN n .
v
La transformada inversa de Hn para recuperar los hk es,
N 1
1 X
hk = Hn e2ikn/N
N
n=0
NOTA: En general, se usa la misma rutina para calcular Hn y
recuperar hk , por lo que, muchas veces hay que normalizar los
resultados de la transformada inversa.
18 / 26
Transformada de Fourier Discreta (cont.)
La version discreta del Teorema de Parseval es,
N 1 N 1
X
2 1 X
|hk | = |Hn |2
N
k=0 n=0
19 / 26
Transformada de Fourier Discreta (cont.)
Cu
anto cuesta calcular Hn ?
Procedimiento original: definir W e2i/N , entonces
N
X 1
Hn = W nk hk
k=0
Equivalente a multiplicar W hk , donde
Wn,k = W nk
Multiplicacion de la matriz requiere N 2 operaciones complejas + otras
operaciones para generar Wn,k .
CONCLUSION: El calculo ingenuo de la DFT es O(N 2 ). Esto es
demasiado para valores grandes de N !!!.
20 / 26
Fast Fourier Transform (FFT)
Danielson and Lanczos (1942), demostraron que:
N
X 1
Fk = e2ijk/N fj
j=0
N/21 N/21
X X
= e2ik(2j)/N f2j + e2ik(2j+1)/N f2j+1
j=0 j=0
N/21 N/21
X X
2ikj/(N/2) k
= e f2j + W e2ikj/(N/2) f2j+1 = Fke + W k Fko
j=0 j=0
Es decir, el calculo de la DFT se puede dividir en el calculo de dos
DFT: una sobre valores pares (Fko ) y otra sobre valores impares (Fke ).
21 / 26
Fast Fourier Transform (FFT) (cont.)
Si N es par, este procedimiento se puede aplicar recursivamente en
log2 N niveles hasta llegar a calcular la DFT de un solo elemento,
Fkeoeoeo...oee = fn for some n
Por ejemplo, en un segundo nivel de recursi
on cada DFT es calculada
sobre N/4 valores considerando elementos del tipo Fkee (impar-impar) y
Fkeo (impar-par).
PROBLEMA: Que n corresponde a la combinacion eoeoeo . . . oee ?
22 / 26
Fast Fourier Transform (FFT) (cont.)
Si N es par, este procedimiento se puede aplicar recursivamente en
log2 N niveles hasta llegar a calcular la DFT de un solo elemento,
Fkeoeoeo...oee = fn for some n
Por ejemplo, en un segundo nivel de recursi
on cada DFT es calculada
sobre N/4 valores considerando elementos del tipo Fkee (impar-impar) y
Fkeo (impar-par).
PROBLEMA: Que n corresponde a la combinacion eoeoeo . . . oee ?
SOLUCION: Cambiar el sentido de la secuencia de es y os y asignar
e = 0 y o = 1, entonces la representaci
on binaria de la secuencia
corresponde al valor de n.
22 / 26
Bit-reversing order
Valores pueden ser copiado y ordenados en un nuevo arreglo o en el
mismo arreglo (in-place).
23 / 26
Bit-reversing order (cont.)
Depues de ordenar los valores, es f
acil calcular la DFT a diferentes
niveles, porque los coeficientes en cada nivel estan almacenados en
forma contigua.
Algoritmo:
1 Ordenar los valores de usando bit-reversing
2 Calcular DFTs para cada nivel. Esto requiere N operaciones por
nivel, y tenemos log2 N niveles N log2 N operaciones.
24 / 26
O(N 2 ) vs O(N log2 N )
Para N = 106 , N 2 = 1012 , N log2 N 2 x 108
Intel Core i7 980 XE puede alcanzar 109 GFLOPS (FLOPS = floating
point operations per second) 25 / 26
Comentarios
I La FFT es considerada como uno (el) de los algoritmos mas
importantes en computaci
on cientfica.
I Es una buena muestra del uso de ideas ingeniosas para acelerar
un calculo.
I La FFT es un buen ejemplo de algoritmos del tipo Divide &
Conquer.
I Aunque en principio la FFT fue desarrollada para N par, tambien
hay implementanciones para N impar que utilizan soluciones
hbridas (FFT + DFT).
I La FFT se aplica en muchos problemas pr acticos: manipulacion de
imagenes (filtros), analisis de se
nales, metodos espectrales,
estadstica (correlaci
on), evaluaci on de convolucion, etc.
26 / 26
