Trasformata di Fourier veloce

In matematica, la trasformata di Fourier veloce, spesso abbreviata con FFT (dall'inglese Fast Fourier Transform), è un algoritmo ottimizzato per calcolare la trasformata discreta di Fourier (DFT) o la sua inversa.

La FFT è utilizzata in una grande varietà di applicazioni, dall'elaborazione di segnali digitali alla soluzione di equazioni differenziali alle derivate parziali agli algoritmi per moltiplicare numeri interi di grandi dimensioni grazie al basso costo computazionale.

Definizione

Sia $x_{0},\dots ,x_{N-1}$ una $N$ -pla di numeri complessi. La DFT è definita dalla formula:

X_{q}={\mathcal {F}}_{d}(x_{n})=\sum _{k=0}^{N-1}x_{k}e^{-i{\frac {2\pi }{N}}kq},\qquad q=0,1,\dots ,N-1.

Calcolare direttamente questa somma richiede una quantità di operazioni aritmetiche $O(N^{2})$ . Un algoritmo FFT ottiene lo stesso risultato con un numero di operazioni $O(N\log(N))$ . In generale questi algoritmi si basano sulla fattorizzazione di $N$ , ma esistono algoritmi FFT per qualunque $N$ , anche per numeri primi.

Poiché l'antitrasformata discreta di Fourier è uguale alla DFT, ma con esponente di segno opposto e $1/N$ a fattore, qualsiasi algoritmo FFT può essere facilmente invertito.

Algoritmo di Cooley-Tukey

L'algoritmo FFT più diffuso è l'algoritmo di Cooley-Tukey. Questo algoritmo si basa sul principio di divide et impera, e spezza ricorsivamente una DFT di qualsiasi dimensione $N$ , con $N$ numero composto tale che $N=N_{1}N_{2}$ in DFT più piccole di dimensioni $N_{1}$ e $N_{2}$ , insieme a $O(n)$ moltiplicazioni per l'unità immaginaria, detti fattori twiddle.

Questo metodo (e in generale l'idea di una trasformata di Fourier veloce) fu reso popolare da una pubblicazione di James Cooley e John Wilder Tukey nel 1965, ma in seguito si scoprì che i due autori avevano indipendentemente reinventato un algoritmo noto a Carl Friedrich Gauss nel 1805^[1] (e successivamente riscoperto in molte altre forme limitate).

L'uso più conosciuto dell'algoritmo di Cooley-Tukey è di dividere e trasformare in due pezzi di $n/2$ ad ogni passo, ed è quindi ottimizzato solo per dimensioni che siano potenze di due, ma in generale può essere utilizzata qualsiasi fattorizzazione (com'era noto sia a Gauss che a Cooley e Tukey). Questi casi sono chiamati rispettivamente casi radice 2 e radice mista (ma esistono anche altri approcci con nomi specifici). Anche se l'idea di base è ricorsiva, la gran parte delle implementazioni tradizionali riorganizzano l'algoritmo per evitare la ricorsione esplicita. Inoltre, poiché l'algoritmo di Cooley-Tukey spezza la DFT in DFT più piccole, può essere arbitrariamente combinato con qualsiasi altro algoritmo per la DFT, come quelli descritti qui sotto.

Un'implementazione iterativa in C++ della FFT basata sull'algoritmo di Cooley-Tukey è la seguente:

#include <iostream>
#include <complex>
#define MAX 200

using namespace std;

int log2(int N)    //funzione per calcolare il logaritmo in base 2 di un intero
{
  int k = N, i = 0;
  while(k) {
    k >>= 1;
    i++;
  }
  return i - 1;
}

int check(int n)    //usato per controllare se il numero di componenti del vettore di input è una potenza di 2
{
  return n > 0 && (n & (n - 1)) == 0;
}

int reverse(int N, int n)    //calcola il reverse number di ogni intero n rispetto al numero massimo N
{
  int j, p = 0;
  for(j = 1; j <= log2(N); j++) {
    if(n & (1 << (log2(N) - j)))
      p |= 1 << (j - 1);
  }
  return p;
}

void ordina(complex<double>* f1, int N)     //dispone gli elementi del vettore ordinandoli per reverse order
{
  complex<double> f2[MAX];
  for(int i = 0; i < N; i++)
    f2[i] = f1[reverse(N, i)];
  for(int j = 0; j < N; j++)
    f1[j] = f2[j];
}

void transform(complex<double>* f, int N)     //calcola il vettore trasformato
{
  ordina(f, N);    //dapprima lo ordina col reverse order
  complex<double> W[N / 2]; //vettore degli zeri dell'unità.
                            //Prima N/2-1 ma genera errore con ciclo for successivo
                           //in quanto prova a copiare in una zona non allocata "W[N/2-1]"
  W[1] = polar(1., -2. * M_PI / N);
  W[0] = 1;
  for(int i = 2; i < N / 2; i++)
    W[i] = pow(W[1], i);
  int n = 1;
  int a = N / 2;
  for(int j = 0; j < log2(N); j++) {
    for(int i = 0; i < N; i++) {
      if(!(i & n)) {
        /*ad ogni step di raddoppiamento di n, vengono utilizzati gli indici */
        /*'i' presi alternativamente a gruppetti di n, una volta si e una no.*/
        complex<double> temp = f[i];
        complex<double> Temp = W[(i * a) % (n * a)] * f[i + n];
        f[i] = temp + Temp;
        f[i + n] = temp - Temp;
      }
    }
    n *= 2;
    a = a / 2;
  }
}

void FFT(complex<double>* f, int N, double d)
{
  transform(f, N);
  for(int i = 0; i < N; i++)
    f[i] *= d; //moltiplica il vettore per il passo in modo da avere il vettore trasformato effettivo
}

int main()
{
  int n;
  do {
    cout << "specifica la dimensione del vettore,che sia potenza di 2" << endl;
    cin >> n;
  } while(!check(n));
  double d;
  cout << "inserisci la dimensione del passo di campionamento" << endl;
  cin >> d;
  cout << "passo di campionamento = " << d << endl;
  complex<double> vec[MAX];
  cout << "inserisci il vettore di campionamento" << endl;
  for(int i = 0; i < n; i++) {
    cout << "inserisci la componente di indice " << i << endl;
    cin >> vec[i];
    cout << "indice " <<i<< "= " << vec[i]  << endl;
  }
  FFT(vec, n, d);
  cout << "vettore trasformato" << endl;
  for(int j = 0; j < n; j++)
    cout << vec[j] << endl;
  return 0;
}

Altri algoritmi per calcolare la FFT

Ci sono altri algoritmi per la FFT oltre al Cooley-Tukey. Per N=N₁N₂ con N₁ e N₂ numeri coprimi può essere utilizzato l'algoritmo di Good-Thomas PFA (Prime-factor FFT Algorithm), basato sul teorema cinese del resto, che fattorizza la DFT in un modo simile al Cooley-Tukey. L'algoritmo FFT di Rader-Brenner è un sistema di fattorizzazione simile al Cooley-Tukey ma con fattori twiddle immaginari puri, riducendo il numero delle moltiplicazioni al costo di un aumento delle addizioni e della instabilità numerica.

Gli algoritmi che suddividono ricorsivamente la DFT in operazioni più piccole diverse dalla DFT includono l'algoritmo di Bruun e il QFT. (Gli algoritmi di Rader-Brenner e il QFT sono stati proposti per dimensioni che siano potenze di 2, ma è possibile che vengano adattati per numeri composti qualunque. L'algoritmo di Bruun si può applicare su dimensioni qualunque anche composte). In particolare l'algoritmo per la FFT di Bruun è basato sull'interpretare la FFT come la fattorizzazione ricorsiva del polinomio zⁿ-1 in polinomi a coefficienti reali nella forma z^m-1 e z^2m+az^m+1.

Un altro punto di vista polinomiale è sfruttato dall'algoritmo FFT di Schmuel Winograd, che fattorizza zⁿ-1 in polinomi ciclotomici, che spesso hanno coefficienti 1, 0 o −1, e quindi richiedono poche (se non nessuna) moltiplicazione, quindi può essere utilizzato per ottenere FFT con un numero minimo di moltiplicazioni ed è spesso usato per trovare algoritmi efficienti per piccoli fattori. In effetti Winograd dimostrò che la DFT può essere calcolata con solo O(n) moltiplicazioni; sfortunatamente questo si ottiene con un numero considerevolmente superiore di addizioni, uno scambio non più favorevole sui moderni processori dotati di chip dedicati per le moltiplicazioni in virgola mobile. In particolare, l'algoritmo di Winograd fa anche uso dell'algoritmo PFA e di quello di Rader per le FFT con dimensioni che siano numeri primi.

Un altro algoritmo FFT per numeri primi è dovuto a L. I. Bluestein ed è a volte chiamata algoritmo chirp-z: riesprime la DFT come una convoluzione, la cui dimensione può essere portata ad una potenza di due e valutata dalla FFT di Cooley-Tukey.

Algoritmi FFT specializzati per dati reali e/o simmetrici

In molte applicazioni i dati di input per la DFT sono reali puri, nel qual caso il risultato soddisfa la simmetria

f_{-j}=f_{j}^{*},

e sono stati conseguentemente sviluppati algoritmi FFT per questa situazione^[2]. Un approccio consiste nel riutilizzare un algoritmo ordinario e rimuovere le parti di calcolo ridondanti, risparmiando approssimativamente un fattore due di occupazione di memoria e di tempo impiegato. In alternativa, è possibile esprimere la trasformata discreta di Fourier di una n-pla con un numero pari di elementi tutti reali come la trasformata discreta complessa della metà della dimensione originaria le cui parti reali ed immaginarie sono gli elementi pari e dispari degli elementi originari, seguiti da O(n) operazioni di post-processo.

Si pensava una volta che le n-ple reali di cui calcolare la DFT potessero essere calcolate in modo più efficiente con la trasformata discreta di Hartley (detta DHT dall'acronimo di Discrete Hartley Transform), ma si scoprì in seguito che in genere possono essere individuati algoritmi FFT specializzati che richiedono meno operazioni del corrispondente algoritmo DHT. Anche l'algoritmo di Bruun fu inizialmente proposto per avvantaggiarsi dei numeri reali come input, ma non è mai diventato popolare.

Accuratezza e approssimazioni

Tutti gli algoritmi FFT presentati finora calcolano esattamente la DFT (in senso aritmetico, trascurando cioè gli errori dovuti ai calcoli in virgola mobile). Tuttavia sono stati anche proposti algoritmi FFT che calcolano la DFT approssimativamente, con un errore che può essere reso arbitrariamente piccolo al costo di un maggiore sforzo computazionale. Questi algoritmi scambiano l'errore di approssimazione a favore di una maggiore velocità od altre caratteristiche. Alcuni esempi sono l'algoritmo FFT di Edelman et al. (1999), la FFT di Guo e Barrus (1996) basata sui wavelet o quella di Shentov et al. (1995).

Anche gli algoritmi FFT "esatti" hanno degli errori quando viene utilizzata l'aritmetica a virgola mobile a precisione finita, ma questi errori sono in genere molto piccoli; la maggior parte degli algoritmi FFT hanno eccellenti proprietà numeriche. Il limite superiore dell'errore relativo per l'algoritmo di Cooley-Tukey è $O(\varepsilon \log n),$ contro $O(\varepsilon n^{3/2})$ per la formula originaria della DFT (Gentleman e Sande, 1966). Questi risultati, comunque, sono molto sensibili verso l'accuratezza dei fattori twiddle usati nella FFT (che sono poi i valori delle funzioni trigonometriche), e non è raro che implementazioni poco accurate della FFT abbiano precisione molto peggiori, ad esempio se utilizzano tabelle dei valori trigonometrici poco accurate. Alcuni algoritmi di FFT, come il Rader-Brenner, sono intrinsecamente meno stabili.

In aritmetica a virgola fissa, la precisione degli errori accumulati dagli algoritmi di FFT sono peggiori, e crescenti come $O({\sqrt {n}})$ per l'algoritmo di Cooley-Tukey^[3]. Inoltre, raggiungere questa precisione richiede un'accurata attenzione e nello riscalare i fattori per minimizzare la perdita di precisione, e gli algoritmi di FFT in virgola fissa richiedono il riscalamento ad ogni stadio intermedio delle decomposizioni come nel Cooley-Tukey.

Per verificare la correttezza di un'implementazione di un algoritmo FFT, garanzie rigorose possono essere ottenute in tempo $O(n\log n)$ con una semplice procedura che controlla le proprietà di linearità, della risposta impulsiva e della tempo-invarianza per dati in input casuali^[4].

Note

^ (LA) Carl Friedrich Gauss, Theoria interpolationis methodo nova tractata, pp. 266-327.
^ H.V. Sorensen, D.L Jones, M.T. Heideman, and C.S. Burrus. (1987, June). Real-valued fast Fourier transform algorithms. IEEE Transactions on Signal Processing, 35(6), 849-863.
^ Digital Signal Processing, A. V. Oppenheim, R. W. Schafer, Prentice Hall, 1975
^ Funda Ergün, 1995, Testing multivariate linear functions: Overcoming the generator bottleneck, Proc. 27th ACM Symposium on the Theory of Computing: 407–416.

Bibliografia

(EN) Gilbert Strang, Wavelets, in American Scientist, vol. 82, n. 3, maggio-giugno 1994, p. 253. URL consultato l'8 ottobre 2013.
(EN) D. F. Elliott, & K. R. Rao, 1982, Fast transforms: Algorithms, analyses, applications. New York: Academic Press.
(EN) Funda Ergün, 1995, Proc. 27th ACM Symposium on the Theory of Computing: 407–416.
(EN) M. Frigo and S. G. Johnson, 2005, "The Design and Implementation of FFTW3," Proceedings of the IEEE 93: 216–231.
(EN) Carl Friedrich Gauss, 1866. "Theoria Interpolationis Methodo Nova Tractata", Werke band 3, 265–327. Göttingen: Königliche Gesellschaft der Wissenschaften.
(EN) W. M. Gentleman and G. Sande, 1966, "Fast Fourier transforms—for fun and profit," Proc. AFIPS 29: 563–578.
(EN) H. Guo and C. S. Burrus, 1996, Proc. SPIE Intl. Soc. Opt. Eng. 2825: 250–259.
(EN) H. Guo, G. A. Sitton, C. S. Burrus, 1994, Proc. IEEE Conf. Acoust. Speech and Sig. Processing (ICASSP) 3: 445–448.
(EN) Steve Haynal and Heidi Haynal, "Generating and Searching Families of FFT Algorithms", Journal on Satisfiability, Boolean Modeling and Computation vol. 7, pp. 145–187 (2011).
(EN) S. G. Johnson and M. Frigo, 2007. "A modified split-radix FFT with fewer arithmetic operations," IEEE Trans. Signal Processing 55 (1): 111–119.
(EN) T. Lundy and J. Van Buskirk, 2007. "A new matrix approach to real FFTs and convolutions of length 2^k," Computing 80 (1): 23-45.
(EN) Kent, Ray D. and Read, Charles (2002). Acoustic Analysis of Speech. ISBN 0-7693-0112-6. Cites Strang, G. (1994)/May–June). Wavelets. American Scientist, 82, 250-255.
(EN) James C. Schatzman, 1996, Accuracy of the discrete Fourier transform and the fast Fourier transform, SIAM J. Sci. Comput. 17: 1150–1166.

Voci correlate

Altri progetti

Wikimedia Commons contiene immagini o altri file sulla trasformata di Fourier veloce

Collegamenti esterni

FFT, su Treccani.it – Enciclopedie on line, Istituto dell'Enciclopedia Italiana.
FFT, in Dizionario delle scienze fisiche, Istituto dell'Enciclopedia Italiana, 1996.
Trasformata di Fourier veloce, in Dizionario delle scienze fisiche, Istituto dell'Enciclopedia Italiana, 1996.
Lorenzo Seno, FFT (Fast Fourier transform), in Enciclopedia della scienza e della tecnica, Istituto dell'Enciclopedia Italiana, 2007-2008.
FFT, su sapere.it, De Agostini.
(EN) fast Fourier transform, su Enciclopedia Britannica, Encyclopædia Britannica, Inc.
(EN) Eric W. Weisstein, Trasformata di Fourier veloce, su MathWorld, Wolfram Research.
(EN) Denis Howe, Trasformata di Fourier veloce, in Free On-line Dictionary of Computing. Disponibile con licenza GFDL
(EN) Fast Fourier Algorithm, su cs.pitt.edu.
(EN) Fast Fourier Transforms, Connexions online book edited by C. Sidney Burrus, with chapters by C. Sidney Burrus, Ivan Selesnick, Markus Pueschel, Matteo Frigo, and Steven G. Johnson (2008).
(EN) Links to FFT code and information online., su fftw.org.
(EN) National Taiwan University – FFT, su cmlab.csie.ntu.edu.tw. URL consultato il 12 gennaio 2014 (archiviato dall'url originale il 30 aprile 2010).
(EN) FFT programming in C++ — Cooley–Tukey algorithm., su librow.com.
(EN) Online documentation, links, book, and code., su jjj.de.
(EN) Using FFT to construct aggregate probability distributions, su vosesoftware.com.
(EN) Sri Welaratna, "Thirty years of FFT analyzers Archiviato il 12 gennaio 2014 in Internet Archive.", Sound and Vibration (January 1997, 30th anniversary issue). A historical review of hardware FFT devices.
(EN) FFT Basics and Case Study Using Multi-Instrument (PDF), su multi-instrument.com.
(EN) FFT Textbook notes, PPTs, Videos at Holistic Numerical Methods Institute.
(EN) ALGLIB FFT Code GPL Licensed multilanguage (VBA, C++, Pascal, etc.) numerical analysis and data processing library.
(EN) MIT's sFFT MIT Sparse FFT algorithm and implementation.
(EN) VB6 FFT VB6 optimized library implementation with source code.
FFTW, libreria software in linguaggio C per la FFT