Fast Fourier Transform
(FFT)
Corso di Metodi per il Trattamento Numerico di
Dati Multimediali
Laura Patricolo (50/174)
Valeria Mele (50/18)
DFT
• f = [ f0, ...., fN-1]T
• per h = 0,…, N-1
C
N 1
F(h)  Fh   f k e
k 0
C
• F = [F0, ...., FN-1]T
2i

hk
N
DFT
Forma matriciale:
F  WN f
Con
WN  ( whk ) h ,k 0,... N 1
whk  ( w)
we

2i
N
h k
2
2
 cos( )  i  sen( )
N
N
( N / 4)
1,1
...
( N / 4)
1, 2
...
( N / 4)
2,1
...
F
F
F (N)
F
( N / 4)
2,1
F
( N / 2)
1
F
Esempio:
è il vettore di dimensione
N/4, ottenuto come prima
parte del vettore
F2( N / 2)
( N / 2)
2
F
F2(,N2 / 4) ...
a sua volta vettore di
dimensione N/2, seconda
parte del vettore
(N)
F
L’operazione di scomposizione intesa con le due frecce indica che nella
prima parte andranno gli elementi di posto pari, e nella seconda quelli di
posto dispari.
Scomposizione formale
N 1
Fh   f k w 
( N 2 ) 1
hk
k 0
f w
N 2 1
k 0
hk
k
f w
N 2 1
k 0
k
hk

k 0
fk w 
N 1
hk
f w
kN 2
hk
k
 f w
 f k ( N 2) wh ( k  ( N 2)) 
per h = 0,…, N-1
N 2 1
k 0
 f
 f N 2 k whk wh ( N 2) 

N 2 1
k 0
k
hk

 f N 2 k whk  h ( N 2) 

h ( N 2)
hk

f
w

w
k
N 2 k
N 1
Fh   f k w 
hk
k 0
f w
N 2 1
k 0
hk
k
f w
N 2 1
k 0
k
hk
( N 2 ) 1

k 0
fk w 
N 1

hk
kN 2
f k whk 
 f w
 f k  ( N 2) wh ( k  ( N 2)) 
N 2 1
k 0
 f
 f N 2 k whk wh ( N 2) 
k
N 2 1
k 0
k
hk

 f N 2 k whk  h ( N 2) 

 f N 2 k wh ( N 2)  whk
Osserviamo:
h ( N 2)
w
e

2 i N
h
N 2
e
ih
 cos( )  i  sin(  )  (1)
h
h
Componenti pari:
f (1)k
(N)
2 h1
F
1  N 21
    f k  f k ( N
2  k 0
2)
w 
2 m1k
1



N.B.: l’apice
accanto a f
indica il passo di
scomposizione
corrente !!
Componenti dispari:
f(1)
k+(n/2)
N 2 1

1
F2(hN1 )1     f k  f k ( N
2  k 0
h1 = 0,…, (N/2)-1
 w  w 
k
2)
2 m1k
1



w1  w2  ei ( 2 /( N 2) )
due trasformate di ordine N/2
f
(1)
k
Trasformando
Vettore delle
componenti di posto
pari di F
f
(1)
k  ( N / 2)
Trasformando
Vettore delle
componenti di
posto dispari di
F
Esempio N=23=8
7
Fh   f k w hk
k 0
F  W8 f
f   f 0 ,..., f 7 
T

F F
(8)
0

(8) T
7
,..., F






W8  






1
1
1
1
1
1
1
1
w
w2
w3
w4
w5
w6
1
w2
w4
w6
w8
w10
w12
1
w3
w6
w9
w12
w15
w18
1
w4
w8
w12
w16
w20
w24
1
w5
w10
w15
w20
w25
w30
1
w6
w12
w18
w24
w30
w36
1
w7
w14
w21
w28
w35
w42
1 

7
w 
w14 
w21 

w28 
w35 

42
w 
w49 
Passo I)
 f 0(1)   f 0  f 4 
 (1)  

 f1   f1  f 5 
 f (1)   f  f 
 2   2 6 
 f 3(1)   f 3  f 7 
 (1)   

 f4   f0  f4 
 f 5(1)   ( f1  f 5 ) w 
 (1)  

2
 f 6   ( f 2  f 6 )w 
 f (1)   ( f  f ) w 3 
 7   3 7

RICORDA: nella notazione di F l’apice esprime la
dimensione del problema, il pedice porta traccia della
scomposizione del vettore, indicando la relazione del
vettore al passo corrente con quello al passo
precedente.
Nella notazione di f invece l’apice esprime il passo
corrente.
NOTA: In questo momento il pedice di F contiene
anche l’indice dell’elemento generico (h1).
F
(8)
F
F
F
F
(8)
2 h1
(8)
2 h1 1
h1= 0, 1, 2, 3
( 4)
1,h1
( 4)
2,h1
Passo II)
 f 0( 2)   f 0(1)  f 2(1) 
 ( 2)   (1)

(1)
 f1   f1  f 3 
 f ( 2)   f (1)  f (1) 
2
 2   0

 f 3( 2)   ( f1(1)  f 3(1) ) w1 
 ( 2)    (1)

(1)
 f4   f4  f6 
 f 5( 2)   f 5(1)  f 7(1) 
 ( 2)   (1)

(1)
 f6   f4  f6 
 f ( 2)   ( f (1)  f (1) ) w 
7
1
 7   5
( 4)
1, h1
F
F2(,4h1)
F
(8)
2 ( 2 h2 )
(8)
2 ( 2 h2 ) 1
(8)
2 ( 2 h2 1)
(8)
2 ( 2 h2 1) 1
F
F
F
F
h2= 0, 1
( 2)
1,1,h2
( 2)
1, 2 ,h2
( 2)
2 ,1,h2
( 2)
2 , 2 ,h2
F
F
F
Passo III)
h3= 0













f 0(3)  
 
( 3)
f1  
f 2(3)  
f 3(3)  

( 3) 
f4  
f 5(3)  
 
( 3)
f6  
f 7(3)  
f 0( 2 )
f 0( 2 )
f 2( 2 )
f 2( 2 )
f 4( 2 )
f 4( 2 )
f 6( 2 )
f 6( 2 )








f1( 2 ) 

( 2)
f1 
f 3( 2 ) 
f 3( 2 ) 
( 2) 
f5 
f 5( 2 ) 

( 2)
f7 
f 7( 2 ) 
F1(,12,)h2
F1(, 22,)h2
F2(,12,)h2
F2(,22),h2
F4((82)h3 )  F0  F1(,11,)1,h3  f 0(3)
F4((82)h3 ) 1  F1  F1(,11,)2,h3  f 4(3)
F2((84)h3 1)  F2  F1(,12),1,h3  f 2( 3)
F2((84)h3 1) 1  F3  F1(,12), 2,h3  f 6(3)
F2((82)( 2 h3 1))  F4  F2(,11),1,h3  f1(3)
F2((82)( 2 h3 1)) 1  F5  F2(,11), 2,h3  f 5(3)
F2((82)( 2 h3 1)1)  F6  F2(,12),1,h3  f 3( 3)
F2((82)( 2 h3 1)1) 1  F7  F2(,12), 2,h3  f 7(3)
schema butterfly
x
y
a
z
b
butterfly
f0
f 0(1)
f1
w
f2 2
f 1(1) w1
f 2(1)
w
f3 3
w
f4
f5 w
f 6 w2
f 7 w3
f
f
f
f
f
f
f
f
(1)
w1
3
(1)
4
(1)
5 w1
(1)
6
(1)
w1
7
f
f
f
f
f
( 2)
0
( 2)
1
( 2)
2
( 2)
3
( 2)
4
( 2)
5
( 2)
6
( 2)
7
f 0( 3)
f 1( 3)
f 2( 3)
f 3( 3)
( 3)
f4
( 3)
f5
f 6( 3)
f 7( 3)
Passo k : N = 8 addizioni/sottrazioni tra elementi del vettore f distanti N/2k
log2N=3 passi
Abbiamo ridotto il numero delle
addizioni/sottrazioni
necessarie al calcolo della nostra DFT
a
N log2N = 8*3 = 24
invece che
N2=64
come richieste dall’algoritmo classico
Notiamo però che il risultato dell’algoritmo è stato:
 F0 
 
 F4 
F 
 2
 F6 
F 
 1
 F5 
 
 F3 
F 
 7













f 0( 3) 

( 3)
f1 
f 2( 3) 
f 3( 3) 
( 3) 
f4 
f 5( 3) 

( 3)
f6 
f 7( 3) 
i coefficienti Fh che abbiamo ottenuto trasformando i
successivi vettori fi non sono posti nell’ordine naturale
Integer Bit Reversal (IBR)
- considera le rappresentazioni binarie degli interi 0, 1, … , N-1;
- inverte le successioni di bit delle singole rappresentazioni;
- determina i corrispondenti numeri decimali.
h
0
1
2
3
4
5
6
7
inversione
000
001
010
011
100
101
110
111
000
100
010
110
001
101
011
111
0
4
2
6
1
5
3
7
ALGORITMO FFT
for m=1, z
m2=2m-1
nm=N/2m
for i=0, m2 step 2
for k=0, nm-1
g= fi*nm+k+ fnm+k
fnm+k=( fi*nm+k - f(i+1)*nm+k)wk
fi*nm+k=g
endfor
endfor
endfor
CONCLUSIONI
La trasformata di Fourier è una tecnica che proprio nelle applicazioni,
dalla fisica dei plasmi alla sismografia, dalla Tac alla oceanografia, alla
ricostruzione di immagini, ha trovato la sua legittimazione come
strumento principale per la risoluzione effettiva di problemi concreti.
Per dare un'idea del guadagno effettivo ricordiamo che la Nasa,
nell'ambito del Progetto Ciclope per la ricerca delle intelligenze
extraterrestri, ha operato la trasformazione di Fourier su un insieme di un
miliardo di dati; ciò ha richiesto circa 9 ore di tempo con la Fft, contro gli
oltre 36.000 anni necessari con l'algoritmo normale.
Ovvero con 1024 punti il rapporto tra i tempi è pari a circa 150.
Bibliografia
[1] – Cooley J. W. e Tukey J. W. , An algorithm for machine calculation of complex Fourier
series, Maths comput., 19, 297 – 301 (1965).
[2] - D. Bini, M. Capovani, G. Lotti, F. Romani "Complessità numerica", ed. Boringhieri
(1981).
[3] - G. Monegato "Fondamenti di Calcolo Numerico", ed. Levrotto&Bella (1990).
[4] - “www.nr.com”, by Numerical Recipes Software, pubblicato da Cambridge University
Press.
[5] - “www.disi.unige.it/person/BoccacciP”, dott. P. Boccacci, Università di Genova.
[6] - “http://digilander.libero.it/nfragale/metodi/capitolo6.html”, Libreria di Software
Matematico, Manuale per l'utente.
[7] - T. Scapolla dell’Università di Pavia, La trasformata di Fourier, Una formula dell'800
usata anche dalla Nasa, Tuttoscienze, inserto del quotidiano La Stampa, (1999) (reperibile
all’indirizzo http://digilander.libero.it/arti2000/ts99/960904.htm)
[8] – “www.wikipedia.org”, Wikipedia, the free encyclopedia.