Capitolo 12 ) Trasformazioni veloci di Fourier

La trasformazione dei valori di un segnale campionato nel tempo in uno spettro di frequenza richiede lo svolgimento di un gran numero di operazioni aritmetiche.

Le sole moltiplicazioni occorrenti per trasformare N campionamenti in N/2 frequenze sono N², e poichè si è visto che per ottenere una buona definizione occorre un N molto elevato, il tempo richiesto per svolgere queste operazioni può costituire un limite alle applicazioni pratiche.

Con un campionamento a 1000 punti si devono svolgere 1 milione di moltiplicazioni, operazioni che richiedono un tempo relativamente elevato rispetto ad altre come le addizioni ( il tempo di esecuzione delle quali diventa quindi trascurabile). Anche se un microprocessore può svolgere una moltiplicazione in un decimillesimo di secondo, sono necessari 100 secondi per svolgere solo le moltiplicazioni, il che rende impossibile l’analisi di una forma d’onda in tempo reale, cioè nello stesso tempo in cui questa viene prodotta.

Un notevole miglioramento può essere ottenuto applicando un procedimento noto come FFT (Fast Fourier Transform) alternativo a quello DFT (Discrete Fourier Transform) visto finora.

Il metodo FFT elimina in un certo modo la ripetizione di moltiplicazioni ridondanti, quindi rende più rapida l’esecuzione.

Si deve evidenziare che l’FFT non è un’approssimazione, ma dà esattamente gli stessi risultati del DFT, eseguendo però solamente (N/2)·log₂(N) moltiplicazioni.

Il vantaggio è più notevole quanto maggiore è N, e precisamente il rapporto fra numero di moltiplicazioni richieste con FFT rispetto a quelle richieste da DFT è (log₂(N)) / (2·N).

Unica limitazione del FFT è che il numero di campionamenti N deve sempre essere un’esatta potenza di 2.

Nell’esempio precedente anzichè 1000 punti dovrebbero essere considerati 1024 punti (2¹⁰), ma si avrebbe una riduzione di circa 200 volte del tempo di esecuzione.

Nell’analisi di Fourier la serie delle armoniche è definita dai valori A_k  dei seni e B_k dei coseni, e passando alla rappresentazione complessa è possibile indicare ciascuna armonica come

Y_k = B_k·cos(2·p·k·f·t) + j·A_k·sen(2·p·k·f·t)

dove f è la fondamentale ( =1/T, con T intervallo di osservazione).

La struttura di calcolo di A_k e B_k per ciascuna armonica k prevede 2 sommatorie di N moltiplicazioni fra gli N campionamenti nel tempo, y_n, per le costanti che rappresentano i valori rispettivamente di seno e coseno delle frequenze k·f in corrispondenza di ciascun istante di campionamento.

Essendo N/2 il numero di armoniche dello spettro base si ha, come già detto, un totale di 2N·N/2 = N² moltiplicazioni.

Il principio su cui si basa il metodo FFT è la reiterata suddivisione per 2 del numero N di punti campionati.

Suddividendo infatti i campionamenti in 2 gruppi di N/2 punti (rispettivamente di ordine pari e di ordine dispari), si possono combinare le corrispondenti uscite del primo e del secondo gruppo in modo da ottenere il calcolo, in termini complessi, delle armoniche. Ma a sua volta ciascun gruppo può essere suddiviso per 2 , e così via fino ad ottenere gruppi di 2 punti.

La struttura di calcolo che permette questa trasformazione è una forma a ‘farfalla’ (butterfly), come indicato nella Fig. 12.1.

Fig. 12.1 - Struttura a farfalla per il calcolo della FFT.

Dati due ingressi (complessi) I_a e I_b ed un coefficente W_k (pure complesso), si ottengono le uscite

U_a = I_a + I_b· W_k e U_b = I_a - I_b· W_k

Applicando questa struttura a ritroso fra le corrispondenti uscite dei gruppi, in log₂(N) passi si completa lo schema di calcolo.

Nell’esempio considerato di N = 8, la struttura completa appare come in Fig. 12.2.

Si deve però osservare che agli ingressi delle struttura la sequenza dei campionamenti risulta alterata. La ‘regola’ di tale sequenza è piuttosto singolare ed è chiamata ‘ordine binario inverso’.

Esprimendo ciascun indice temporale n in numero binario, quindi con una successione di bit, per ricavare il nuovo indice corrispondente si deve invertire l’ordine di tale successione.

Nell’esempio considerato l’indice di y₁ in codice binario è 001, ed invertendone l’ordine si ottiene 100 (cioè 4), dunque al secondo posto della struttura (dopo y₀) dovrà essere applicato y₄, e y₁ andrà al quinto posto.

Così dovranno essere scambiati anche y₃ e y₆, mentre gli indici simmetrici non subiscono variazioni.

Le costanti W_k = cos(2pkf/N) + j·sen(2pkf/N) rappresentano infine i coefficienti per cui sono moltiplicati i valori di campionamento.

Fig. 12.2 - Struttura di calcolo per FFT

La procedura di calcolo richiede di eseguire tutte le ‘farfalle’ del primo passo, poi con i risultati di queste si devono eseguire quelle del secondo passo ed infine quelle del terzo.

La Fig. 12.3 riporta un esempio numerico di tale procedura, con le tabelle dei risultati intermedi z_n,p ad ogni passo.

I risultati dell’ultimo passo devono poi essere normalizzati (moltiplicati per N/2) per essere confrontabili con i valori di A_k e di B_k ottenuti con il normale procedimento DFT.

Vi è da notare che il primo valore, cioè Y₀, dovrebbe corrispondere alla componente continua C. In realtà questo risulta essere il doppio in entrambi i casi.

La spiegazione è che ricavando la componente continua come coefficienti di Fourier A₀ e B₀ , cioè iniziando da k=0 anzichè da k=1, si ottiene B₀=2C.

Infatti A₀ risulta sempre uguale a zero (è la somma dei vari campionamenti per sen(0) ) e B_k =(2/N)·Sy_n·cos(0), cioè il doppio della media dei campionamenti.

Si può inoltre osservare che il metodo calcola uno spettro doppio di quello base (N frequenze anzichè N/2), pur non richiedendo moltiplicazioni in più, ma solo differenze.

In effetti questo è dovuto al fatto che i campionamenti sono numeri reali, mentre il metodo considera numeri complessi. Si può trarre ulteriore vantaggio da questa osservazione sia ‘mescolando’ due segnali reali, ciascuno con N campionamenti, ed elaborandoli nello stesso FFT come se uno fosse reale e l’altro immaginario, oppure dividendo per due il numero dei 2N campionamenti dello stesso segnale (in pari e dispari) e trattando le due serie una come reale e l’altra come immaginaria.

Ottenuto lo spettro totale (di N frequenze) si possono ricostruire gli spettri singoli (di N/2 frequenze) sfruttando delle proprietà di simmetria intrinseche nelle componenti reali ed immaginarie delle armoniche corrispondenti, ‘riflesse’ attorno alla frequenza zero. Un esempio sarà dato alla fine di questo capitolo.

L’illustrazione del metodo FFT che, come abbiamo visto, risulta di difficile implementazione, è stata qui data solo per conoscenza: in pratica non occorre procedere al suo svolgimento in quanto già sono disponibili per l’applicazione numerose versioni standard, sia hardware che software.

In particolare in ambiente Mathcadâ, qui utilizzato per il calcolo e la grafica dei vari esempi riportati nelle figure, esistono particolari funzioni di conversione che permettono l’immediata trasformazione dal tempo alla frequenza e viceversa.

Stabilita una variabile indicizzata (array) y_n che rappresenta i campionamenti nel tempo, la serie di frequenze Y_k è ottenuta semplicemente scrivendo Y=FFT(y).

Fig. 12.3 - Esempio di calcolo FFT

Il risultato Y_k , con k variable da 0 ad N/2-1, fornisce in forma complessa l’ampiezza di ciascuna armonica dello spettro base (le armoniche fuori da questo spettro, come ormai noto, non sono altro che una duplicazione speculare di queste).

In particolare si ha: A_k = 2·Im(Yk), cioè il doppio della parte immaginaria di Y_k , e B_k = 2· Re(Y_k) , cioè il doppio della parte reale di Y_k .

La componente continua è invece C= Re(Y₀).

La Fig. 12.4 dà un esempio di tali trasformazioni, in cui i campioni sono generati dalla somma di sinusoidi e cosinusoidi di cui è fissata l’ampiezza. Come è facile vedere le componenti immaginarie e reali di Y_k ricavate dalla trasformazione Mathcadâ FFT, corrispondono alla metà dei valori assegnati A_k e B_k .

Fig. 12.4 - Esempio di applicazione FFT ad una forma d’onda costruita con sinusoidi e cosinusoidi di ampiezza nota.

Viceversa dato lo spettro di frequenze Y_k si ottiene l’andamento nel tempo con la trasformazione inversa y=IFFT(Y). Naturalmente l’indice k deve essere variabile fra 0 e N/2-1, dove N è sempre una potenza di 2.

Si è già detto che il metodo FFT lavora con numeri complessi e che N campioni complessi forniscono N frequenze. Si è pure detto che partendo invece da N

Fig. 12.5a - Applicazione della FFT complessa al calcolo dello spettro della somma di due funzioni.

campionamenti reali si può semplificare il calcolo ottenendo N/2 frequenze, cioè lo spettro base, come si è appena visto.

Il Mathcadâ consente anche il calcolo con campioni complessi, semplicemente cambiando il nome della funzione. Con Y=CFFT(y) si elaborano gli N valori complessi y_n ottenendo un spettro base di N frequenze.

Ciò permette ad esempio di elaborare due distinte forme d’onda g_n e h_n , campionate negli stessi tempi t_n , con un solo calcolo FFT.

La Fig. 12.5 costituisce un esempio di tale procedura.

Fig 12.5b - Verifica della procedura illustrata nella figura precedente.

Naturalmente risulta più semplice il calcolo diretto, cioè separato delle due funzioni, come appare nella verifica, ma l’esempio è stato riportato per dare l’idea della possibilità di ottimizzazione delle procedure (soprattutto nei casi di realizzazioni in hardware) con conseguenti diminuzioni dei tempi di calcolo globali.