PONTIFICIA UNIVERSIDAD CATÓLICA DE CHILE
Escuela de Ingenierı́a
Departamento de Ingenierı́a Estructural y Geotécnica
ICE 3723 Dinámica de Estructuras
Algoritmo de la Transformada rápida de Fourier (FFT)
La transformada de Fourier se define como
Z ∞
F (ω) = f (t)e−iωt dt (1)
−∞
Sin embargo, si en vez de tener una función continua en el tiempo tenemos un registro discreto de datos,
como ocurre en un movimiento sı́smico, se tiene para un tiempo tn que
T
f (tn ) = xn tn = n∆t ∆t =
N
donde T es el tiempo total del registro, N es la cantidad de muestras, xn es el valor de la función en el tiempo
tn y n = 0, 1, · · · , N − 1. Ası́, la transformada de Fourier en tiempo discreto (DTFT) queda definida como
N −1
T X T
FDT F T (ω) = xn e−iωn N (2)
N n=0
Ahora, si además calculamos la transformada para sólo algunos valores de ω, es decir se discretiza en la
dimensión de la frecuencia, se tiene
2π
ωm = mω = m
T
ası́, la transformada discreta de Fourier (FT) es
N −1
T X 2π
FF T [ωm ] = xn e−imn N (3)
N n=0
Pese a que la transformada se puede calcular para un m en particular, la idea es calcularla para m = n, lo
cual implicarı́a N 2 operaciones. Para eso empleamos la transformada rápida de Fourier (FFT) que consiste
en dividir el conjunto de términos en dos conjuntos, uno de pares y otro de impares. Ası́, si trabajamos sólo
con la sumatoria, tenemos
N N
N −1 −1
2 −1
2
−imn 2π −im2j 2π 2π
X X X
Xm = xn e N = x2j e N + x2j+1 e−im(2j+1) N
n=0 j=0 j=0
N N
2 −1 2 −1
X −imj 2π 2π
X −imj 2π
+ e− N im
N N
= x2j e 2 x2j+1 e 2 (4)
j=0 j=0
Si llamamos x2j = yj y x2j+1 = zj , nos queda
N N
−1
2 −1
2
−imj 2π −imj 2π
− 2π
X N
X N
Xm = yj e 2 +e N im zj e 2
j=0 j=0
2π
Xm = Ym + e− N im Zm (5)
1 En caso de identificar cualquier error, se solicita por favor informarlo.
1
donde Ym y Zm son las transformadas de Fourier discretas de los términos pares e impares respectivamente.
Por otro lado, como buscamos los valores m = n podemos separar la transformada en dos: Xm y X N +m ,
2
donde m = 0, 1, · · · , N2 −1. Claramente la primera expresión no cambia, sin embargo si calculamos la segunda:
N N
−1
2 2−1
−i( N
2 +m)j N
2π
−i( N
2 +m)j N
2π
N i( 2 +m)
− 2π
X N X
X N +m = yj e 2 +e zj e 2
2
j=0 j=0
N N
2 −1 −1
2
−imj 2π −imj 2π
−3πi − 2π
X X
−2πi N
N im
N
= e yj e 2 +e e zj e 2
j=0 j=0
N N
−1
2 2−1
−imj 2π −imj 2π
− 2π
X N
X N
= yj e 2 −e N im zj e 2
j=0 j=0
2π
= Ym − e− N im Zm
Luego, la transformada rápida de Fourier se puede resumir como
2π
Xm = Ym + e− N im Zm
2π (6)
X N +m = Ym − e− N im Zm
2
para m = 0, 1, · · · , N2 − 1. Implementando este algoritmo en MatLab, conocido también como el algoritmo
de Cooley-Tukey, tenemos
1 function X= FFT(x) % Fast Fourier Transform
2 N= length(x);
3 if N==1 % Caso trivial
4 X= x;
5 else % N términos, N=2ˆk, k>0
6
7 X= zeros(N,1); % Transformada
8 xp= x(2*(1:N/2)−1); % Pares del input
9 xnp= x(2*(1:N/2)); % Impares del input
10
11 Y= FFT(xp); % FFT para los términos pares
12 Z= FFT(xnp); % FFT para los términos impares
13
14 % Combinamos FFT
15 X(1:N/2)= Y(1:N/2)+exp(−2*pi*1j*((1:N/2)'−1)/N).*Z(1:N/2);
16 X(N/2+(1:N/2))= Y(1:N/2)−exp(−2*pi*1j*((1:N/2)'−1)/N).*Z(1:N/2);
17
18 end
Como vemos, este algoritmo divide cada conjunto de términos en 2, por lo que se requiere que el número
total de términos sea una potencia de 2 y ası́ sólo realizar N log N operaciones. Por otro lado, podemos notar
que los términos a partir de N2 en realidad corresponden a ((frecuencias negativas)), por lo que, o se grafica
sólo la primera mitad, o se reordenan.
2
Aplicándolo para una señal cualquiera de 16 puntos tenemos que, reordenando el vector resultante, el
gráfico del FAS (módulo de la transformada de Fourier) es
FAS 5
0
−8 −6 −4 −2 0 2 4 6
ω
Figura 1: FAS para un registro random
La señal utilizada fue
0,7898
0,1023
0,7679
0,6621
0,7731
0,9320
0,1769
0,0204
x=
0,8336
0,4732
0,1086
0,0891
0,1151
0,1680
0,9725
0,2860
Vale destacar que todo este análisis se realizó sin considerar el ∆t, el cual se multiplica al final de todos los
cálculos.