Fast Fourier Trasform: matlab

La trasformata di Fourier veloce (spesso indicata come FFT, dall'inglese Fast Fourier Transform) è un algoritmo ottimizzato per calcolare la trasformata discreta di Fourier (detta DFT) e la sua inversa. La FFT è di grande importanza per 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.
Sia x₀, ..., x_N-1 una n-pla di numeri complessi. La DFT è definita dalla formula

per q = 0, 1, ..., N-1

Calcolare direttamente questa somma richiede una quantità di operazioni aritmetiche O(N²). 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 (questi ci piacciono di più).
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.

L'utilizzo della FFT in programmazione, richiede la conoscenza di alcune basi (veramente semplici) dei segnali in frequenza. Ci concentriamo qua sulla descrizione di una semplice funzione che calcoli la FFT di una serie di dati.

Dopo questa semplice intro sulla FFT, ecco una funzione in matlab abile di calcore e plottare i dati trasformati in frequenza. Gli input da fornire sono il vettore dati e la frequenza di campionamento. L'utilizzo della funzione è semplicissimo, sia V il vettore dati e Fs la frequenza di campionamento, si digiti:

z=calcFFT(V,Fs,type)

type è un indice da 1 a 4 che permette di scegliere il tipo di output:

1) plot amplitude  

2) plot phase       

3) plot aplitude and phase (using subplot)       

4) power spectral density

Il codice utilizzato, è molto semplcie, basta salvarlo con estension .m ed aggiungere la directory dove si trova il codice al path di matlab. Ed ecco fatto un semplice programmino che calcoli e plotti la FFT per una qualsiasi serie di dati, e ci dia in automatico il plot di questi.

CODICE:

%% computation and direct plotting of fft
% input parameters:
%
%   signal= input signal to trasform
%   Fs= sample frequency
%   type= 1) plot amplitude
%         2) plot phase
%         3) plot aplitude and phase (using subplot)
%         4) power spectral density
%


function z=calcFFT(signal,Fs,type);




z=fft(signal);

F=(0:(numel(signal))-1)*Fs/numel(signal);

if nargin<3
 type=1;  
end

if type==1
plot(F,abs(z));
xlim([0 Fs/2])
xlabel('Frequency (hz)')
ylabel('Amplitude')
end

if type==2
  plot(F,imag(z));
xlim([0 Fs/2])
xlabel('Frequency (hz)')
ylabel('Phase')
end

if type==3
  subplot(121)
    plot(F,imag(z));
xlim([0 Fs/2])
xlabel('Frequency (hz)')
ylabel('Phase')
  subplot(122)
plot(F,abs(z));
xlim([0 Fs/2])
xlabel('Frequency (hz)')
ylabel('Amplitude')

end  

if type==4
plot(F,abs(z).^2);
xlim([0 Fs/2])
xlabel('Frequency (hz)')
ylabel('Energy')
title('Power spectral density')
end   

MORE: 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₁N₂ in DFT più piccole di dimensioni N₁ e N₂, 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 William 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 (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 di 2 e radice mista (e le altre varianti hanno i propri nomi. 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.

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 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

e sono stati conseguentemente sviluppati algoritmi FFT per questa situazione (cfr. Sorensen, 1987). 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 FFT (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(ε log n), contro O(ε 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(√n) per l'algoritmo di Cooley-Tukey (Oppenheim & Schafer, 1975). 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 (Ergün, 1995).