REVISTA DE LA ESCUELA DE FÍSICA, UNAH • Vol. 5, No.

1 • 6-10 6

Trasformada rápida de Fourier utilizando Python

Michael Spilsbury1 y Armando Euceda2
1
Escuela de Fı́sica - UNAH , mail: [Link]@[Link]
2
Escuela de Fı́sica - UNAH, mail: aeunah@[Link]

Recibido: 28 de febrero del 2017 / Aceptado: 30 de abril del 2017

Resumen

A continuación se presenta un programa de computadora para calcular la transformada discreta de Fourier utilizando
el algoritmo de la transformada rápida de Fourier (FFT por sus siglas en inglés). Desde 1965[2], cuando James
W. Cooley y John W. Tukey publicaron dicho algoritmo, su uso se ha expandido rápidamente y las computadoras
personales han impulsado una explosión de aplicaciones adicionales de la FFT. Como lenguaje de programación
se usará Python, que es un lenguaje de programación multiparadigma, esto significa que más que forzar a los
programadores a adoptar un estilo particular de programación, permite varios estilos: programación orientada
a objetos, imperativa y funcional. Los usuarios de Python se refieren a menudo a la filosofı́a Python que es
bastante análoga a la filosofı́a de Unix. Al mismo tiempo se integrará un módulo de Fortran para mejorar el desempeño.

Palabras clave: Transforma rápida de Fourier, algoritmo, transformada discreta de Fourier, Cooley, Tu-
key, programació, Python, Fortran.

Next a computer program is presented to calculate the discrete Fourier transform using the fast Fourier transform
algorithm. Since 1965[2], when James W. Cooley and John W. Tukey published this algorithm, its use has expanded
quickly and personal computers have generated an explosion of further FFT applications. The programming
language to use is Python, which is a multi-paradigm programming language, this means that rather than forcing
programmers to adopt a particular style of programming allows several styles: object-oriented, functional and
imperative programming. Python users often refer to the Python philosophy is quite analogous to the Unix philosophy.
At the same time Fortran module will be integrated to improve performance.

Keywords: Fast fourier transform, algorithm, discrete Fourier transform, Cooley, Tukey, programming,
Python, Fortran.

I. Programa en Python para cálculos Como se indica, Numpy contiene dentro de sus venta-
usando FFT jas varias funciones para transformadas de Fourier, que
nos servirán para comparar los tiempos de cálculo y los
omo se apuntó inicialmente, el lenguaje base pa- datos obtenidos. Para conocer sobre Python, ver [3]

C ra la programación de la FFT en este trabajo es

Python, que es un lenguaje interpretado y debido
a esto, cuando se tiene un gran número de muestras la
II. Programa FFT
velocidad de procesamiento se ve reducida; por lo que para Se desarrollará a continuación el programa en códi-
dar velocidad y poder de procesamiento a los cálculos se go Python, siguiendo el diagrama de flujo expuesto por
utilizará un módulo de Fortran (lenguaje compilado) que Brigham[1] (Fig. 1). Cabe mencionar que en el desarrollo
será utilizado por Python para realizar los cálculos. se considera que el número de muestras, N , es una po-
Python tiene desarrolladas muchas extensiones que tencia de 2 (N = 2γ ). Además recordar que el diagrama
le dan mayor versatilidad, una de ellas es Numpy que surge de la aplicación de la transformada discreta de Fou-
es el paquete fundamental para el cálculo cientı́fico, que rier y sus propiedades, y en sı́ el algoritmo no requiere
contiene entre otras cosas: una análisis profundo ya que solamente ayuda a calcular
más rápida y eficientemente la transformada discreta de
Un poderoso arreglo de objetos N-dimensionales.
Fourier, por lo tanto no se dará mayor detalle. Se puede
Funciones (transmitidas) sofisticadas. consultar la referencia antes citada.
Herramientas para la integración de código C/C++ y
Fortran.
Útiles funciones de álgebra lineal, transformada de
Fourier y números aleatorios.

REF-UNAH / Vol. 5 - No. 1 / 6-10

 REVISTA DE LA ESCUELA DE FÍSICA, UNAH • Vol. 5, No. 1 • 6-10 7

Comienzo

Datos Entrada
x(k), k = 0, 1, . . . , N − 1;
N = 2γ , γ es un entero;
NU = γ

Inicializadores
l = 1;
N 2 = N/2;
N U 1 = γ - 1;
k =0

SI ¿Es
i = I BR(k)
l > γ?

k = k + 1 NO
¿Es SI
i ≤ k?

NO I = 1

T 3 = x(k); ¿Es
x(k) = x(i); k = N − 1? NO
x(i) = T 3
l = l + 1;
SI M = Entero de k/(2N U 1 ); 2 N 2 = N 2/2;
P = I BR(M ) N U 1 = N U 1 − 1;
Final k =0

T 1 = W P x(k + N 2);
x(k + N 2) = x(k) − T 1;
x(k) = x(k) + T 1

I = I + 1 k = k + 1

SI

NO ¿Es SI ¿Es NO
Subrutina de Inversión de Bit k = k + N2
I = N 2? k < N − 1?

Entrada
M, N U

Inicio
¿Es NO J 2 = M/2;
I 1 = 1; I1 = I1 + 1
I1 > N U ? I BR = 2 ∗ I BR + (M − 2 ∗ J 2);
I BR = 0 M = J2

SI

Retornar

Figura 1: Diagrama de flujo del programa FFT. Fuente: Brigham[1]

Para ejemplo se utilizarán muestras de la fun- 22 T = nT ∗ T 0

ción A sen(2πf0 t), tomando A = 1, f0 = 250 Hz y 23 delta t = 1 / n ##f a c t o r e n t r e TFC y TFD
d e l t a t 1 = T / n ##D e l t a t d e l m u e s t r e o
N = 210 = 1024 muestras:
24
25 delta f = 1 / T ##D e l t a f d e l m u e s t r e o
26 t = np . a r r a y ( d e l t a t 1 ∗ k m u e s t r a s )
27 x r = np . s i n ( 2 ∗ np . p i ∗ f 0 ∗ t )
1 ##Transformada r á p i d a de F o u r i e r 28 xi = [ 0 ] ∗ n
2 from future im por t p r i n t f u n c t i o n , d i v i s i o n 29 xz = np . a r r a y ( x r ) + 1 j ∗ np . a r r a y ( x i )
3 from p y l a b im por t ∗ 30
4 31 ##I n i c i a l i z a r v a r i a b l e s a u x i l i a r e s
5 ##V a r i a b l e s de e n t r a d a 32 l = 1
6 ##x r = p a r t e r e a l de d a t o s 33 n2 = np . i n t ( n / 2 )
7 ##x i = p a r t e i m a g i n a r i a de d a t o s 34 nu1 = nu − 1
8 ## W = exp (−1 j ∗ 2 ∗ np . p i / N) 35 k = 0
9 ##T1 = v a l o r t e m p o r a l (Wˆp ) ∗x ( k + N2) 36
10 x r = [ ] 37 ##I m p l e m e n t a c i ó n Python−I n i c i o
11 x i = [ ] 38 xza = np . a r r a y ( xz )
12 nu = 10 39 w h i l e l <= nu :
13 n = np . power ( 2 , nu ) 40 w h i l e k <= n − 1 :
14 41 q = 1
15 ##Datos de p r u e b a s 42 w h i l e q <= n2 :
16 k m u e s t r a s = np . a r a n g e ( n ) 43 M = np . i n t ( k / np . power ( 2 , nu1 ) )
17 44 P =
18 ##F r e c u e n c i a y p e r i o d o f u n c i ó n s e n o np . i n t ( ’ { : 0 { width }b} ’ . f o r m a t (M, width = nu )
19 f 0 = 250 ##F r e c u e n c i a 45 [:: −1] ,2)
20 T 0 = 1 / f 0 ##P e r i o d o 46 W = np . exp(− 1 j ∗ 2 ∗ np . p i / n )
21 nT = 20 ##Veces e l p e r i o d o a g r a f i c a r

REF-UNAH / Vol. 5 - No. 1 / 6-10

 REVISTA DE LA ESCUELA DE FÍSICA, UNAH • Vol. 5, No. 1 • 6-10 8

47 T1 = xza [ k + n2 ] ∗ W∗∗P Gráfico función del tiempo

48 xza [ k + n2 ] = xza [ k ] − T1 1.0
49 xza [ k ] = xza [ k ] + T1 0.5

Amplitud
50 k += 1 0.0
51 q +=1 0.5
52 k += n2 1.0
53 l += 1 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
54 n2 = np . i n t ( n2 / 2 ) tiempo (s)
55 nu1 −= 1 0.6 Gráfico función de la frecuencia
56 k = 0 0.4
57 w h i l e k < n − 1 : 0.2

Amplitud
58 i = np . i n t ( ’ { : 0 { width }b} ’ . f o r m a t ( k , width 0.0
= nu ) 0.2
59 [:: −1] ,2) 0.4
60 if i > k: 0.6

-4000

-2000

-250
0
250

2000

4000
61 T3 = xza [ k ]
62 xza [ k ] = xza [ i ] Frecuencia (Hz)
63 xza [ i ] = T3
64 k +=1 Figura 2: Gráfico de la función A sen(2πf0 t) y su transfor-
65 xza ∗= d e l t a t
mada. A = 1 y f0 = 250 Hz
66 ##I m p l e m e n t a c i ó n Python−F i n a l
67
68 ##A j u s t e de d a t o s de f r e c u e n c i a Gráfico función del tiempo
69 f = [] 1.0
70 f o r s i n r a n g e ( l e n ( k m u e s t r a s ) ) :
Amplitud 0.5
71 i f s <= ( n / 2 + 1 ) : 0.0
72 f . append ( s ∗ d e l t a f ) 0.5
73 else :
1.0
74 f . append ( ( s − n ) ∗ d e l t a f ) 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
tiempo (s)
Se puede apreciar que en el programa en Python no se 0.6 Gráfico función de la frecuencia
consideró la subrutina de inversión de bit ya que se utilizó
0.4
una de las muchas ventajas de Python para realizar esto
Amplitud

0.2
en una sola lı́nea de comando, siendo estas lı́neas la 44 y
0.0
la 58 (estas lı́neas son la caja 1 y la caja 2 del diagrama
de flujo, que utilizan la subrutina de inversión de bit). -0.2
-4000

-2000

-250
0
250

2000

4000
Además la lı́nea 12 corresponde a la potencia de la base 2 Frecuencia (Hz)
que genera el número de muestras tomadas en la función
antes mencionada, definiendo parámetros y muestras de Figura 3: Gráfico de la función A cos(2πf0 t) y su transfor-
la misma entre las lı́neas 19 a la 29. Las lı́neas entre la 31 mada. A = 1 y f0 = 250 Hz
a la 66 corresponden al programa indicado en el diagrama
de flujo de la Fig. 1. La lı́nea 65 corresponde al factor de Precisamente, ya que la transforma de Fourier en ge-
escala resultante del análisis de la transformada discreta neral es una función compleja
de Fourier (Brigham[1]). Las lı́neas entre la 69 a la 74
corresponden a los ajustes de los datos de frecuencia con Z∞
los que se representa la transformada contı́nua de Fourier, H (f ) = h(t)e−j2πf t dt
siendo los lı́mites de −f0 a f0 . −∞
En la Fig. 2 se muestran los datos obtenidos para el = R(f ) + jI (f ) = |H (f )|ejθ (f ) (3)
caso mencionado [A sen(2πf0 t)], y al comparar con el par
transformado para esta función en el programa la parte real [R(f )] deberá graficarse con
lı́nea de color diferente o tipo de lı́nea a la de la parte
A A
A sen(2πf0 t) ⇔ −j δ (f − f0 ) + j δ (f + f0 ) (1) imaginaria [I (f )] [en estos gráficos R(f ) e I (f ) aparecen
2 2
con color de lı́nea azul y verde respectivamente]. Para
se observa coincidencia con los resultados arrojados por ejemplo tómese la función:
el programa.
βe−αt t > 0

Ademas en la Fig. 3 se observa el caso para h(t) = (4)
A cos(2πf0 t), tomando A = 1, f0 = 250 Hz y N = 0 t<0
210 = 1024 muestras y cuyo par transformado de Fourier cuya transformada de Fourier es:
es
A A
A cos(2πf0 t) ⇔ δ (f − f0 ) + δ (f + f0 ) (2) H (f ) = R(f ) + jI (f )
2 2
βα j2πf β
Nótese el color de lı́nea diferente en el gráfico de la trans- = 2 − (5)
α + (2πf )2 α2 + (2πf )2
formada de Fourier de la Fig. 3, en contraste con el co-
rrespondiente de la Fig. 2. En la Fig. 4 se muestra el gráfico correspondiente.

REF-UNAH / Vol. 5 - No. 1 / 6-10

 REVISTA DE LA ESCUELA DE FÍSICA, UNAH • Vol. 5, No. 1 • 6-10 9

1.2 Gráfico función del tiempo

1.0
0.8
Amplitud

0.6
0.4
0.2
0.0
0.0 0.2 0.4 0.6 0.8 1.0
tiempo (s)
Gráfico función de la frecuencia
0.10 R(f)
I(f)
0.05
Amplitud

0.00

0.05
100 50 0 50 100
Frecuencia (Hz)
Figura 4: Gráfico de e−αt para t > 0 y su transformada.
α = 10 Hz

III. Análisis de resultados del suma de las diferencias entre el número de muestras mul-
tiplicado por 100):
programa
Revisando el tiempo de ejecución del programa utili- Diferencia mayor en los datos = 0.000607593968998
zando código Python únicamente, le tomaba alrededor de Error relativo = 0.0178151019639
200 veces el tiempo que le toma a la extensión Numpy
de Python. Esta diferencia se debe a que Numpy realiza A continuación se muestran las figuras obtenidas los
el cálculo basado en Fortran, por lo que se elaborará un tres métodos
modulo de Fortran que es llamado por Python y de esta
Para que Python llame el módulo de Fortran, en pri-
forma reducir el tiempo de ejecución de las rutinas de
mer lugar se escribio en lenguaje Fortran el programa
cálculo. Con esto se obtendrá una mejora notable, casi
(siguiendo el diagrama de flujo) y luego con el paquete de
igual al tiempo de la extensión Numpy. Por ejemplo, para
Pyhton F2PY se genera dicho módulo. Los pasos seguidos
la función pulso cuadrado, que tiene el par transformado.
fueron:

 A |t| < T0
Generar un archivo de firma (signature file), usando
h(t) = A/2 |t| = T0 ⇔
el comando:
0 |t| > T0

f2py mod-fpy.f -m mod-fpy -h [Link]
sen (2πT0 f ) El archivo con extensión .pyf es el generado.
H (f ) = 2AT0 (6)
2πT0 f
Generar el módulo con el comando:
con un número de muestras igual a 1024, se realizó el
calculo por los tres métodos obteniéndose los siguientes f2py -c [Link] mod-fpy.f95
tiempos de ejecución:
lo anterior se hace desde la lı́nea de comandos del sistema
Tiempo implementación Python: 0.064370020345 s operativo. Con la generación del módulo en el programa
Tiempo implementación Python_Fortran: 0.000663668916201 s se debe agregar en la lı́nea 4:
Tiempo implementación FFT_Numpy: 0.000339165116232 s
4 import mod-fpy
En cuando a los errores obtenidos respecto de la fun-
ción continua, arrojó (el error relativo se calculó de la y se deben reemplazar las lı́neas de la 31 a la 66 por:

REF-UNAH / Vol. 5 - No. 1 / 6-10

 REVISTA DE LA ESCUELA DE FÍSICA, UNAH • Vol. 5, No. 1 • 6-10 10

32 ###################################### Referencias
33 ##Implementacion Python_Fortran-Inicio
34 ###################################### [1] Brigham, E.O., Ed. ( 1988).
c The fast Fourier
35 mod-fpy.fft_fortran.nu = nu
36 mod-fpy.fft_fortran.n = n
transform and its applications. Prentice Hall, En-
37 mod-fpy.fft_fortran.xzf = list(xz) glewood Cliffs, N.J. ISBN 0-13-307505-2. URL
38 mod-fpy.fft_fortran.fft_it() [Link]
39 spxz = mod-fpy.fft_fortran.xzf
40 spxz *= delta_t [2] Cooley, J.W. y Tukey, J.W., An algorithm
41 ##Implementación Python_Fortran-Final
for the machine calculation of complex Fourier se-
donde fft fortran y fft it es el nombre del módulo y la ries. En Mathematics of Computation, Volumen 19
subrutina dentro del programa en Fortran. A continuación (1965) (90): 297. ISSN 0025-5718. doi:10.1090/
se presenta el archivo de firma: S0025-5718-1965-0178586-1.
1 ! −∗− f 9 0 −∗−
2 ! Note : t h e c o n t e x t o f t h i s file i s case [3] Langtangen, H.P. ( 2012).
c A primer on scienti-
sensitive . fic programming with Python, Volumen 6 desde Texts
3 in computational science and engineering. Springer,
4 python module mod−f p y ! i n
Berlin and New York, 3rd ed enlugar. ISBN 978-3-
5 interface ! i n : mod−f p y
6 module f f t f o r t r a n ! i n 642-30292-3.
: mod−f p y : mod−f p y . f 9 5
7 complex ( k i n d =8) ,
a l l o c a t a b l e , dimension ( : ) : : x z f
8 i n t e g e r : : nu
9 i n t e g e r , parameter , o p t i o n a l : :
dp=s e l e c t e d r e a l k i n d ( 1 5 , 3 0 7 )
10 integer : : n
11 subroutine f f t i t ! in
: mod−f p y : mod−f p y . f 9 5 : f f t f o r t r a n
12 end s u b r o u t i n e f f t i t
13 f u n c t i o n i b r ( val m , v a l n u ) ! i n
: mod−f p y : mod−f p y . f 9 5 : f f t f o r t r a n
14 i n t e g e r : : val m
15 integer : : val nu
16 integer : : ibr
17 end f u n c t i o n i b r
18 end module f f t f o r t r a n
19 end i n t e r f a c e
20 end python module mod−f p y
21
22 ! This f i l e was auto−g e n e r a t e d with f 2 p y
( version :2) .
23 ! See h t t p : / / c e n s . i o c . e e / p r o j e c t s / f 2 p y 2 e /

IV. Conclusiones
1. El desarrollo del programa por medio del algoritmo,
resulta bastante intuitivo, indiferentemente del len-
guaje de programación utilizado. Al mismo tiempo
no se requiere mayor detalle, como se apunto, de la
teorı́a de la transformada de Fourier.
2. El tiempo del cálculo se vio reducido considerable-
mente con el uso del lenguaje compilado (Fortran),
mejorando el desempeño del programa. Sin embar-
go, aún se puede reducir el tiempo utilizando lo
expuesto por Brigham, en el caso que la función del
tiempo sea real [1, Sección 9.3].

REF-UNAH / Vol. 5 - No. 1 / 6-10