Trasformazioni  veloci  di  Fourier (12)

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.

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.

Analisi circuitali con FFT (13)

Dopo aver visto la semplicità formale per l’uso delle  FFT,  è opportuno vederne subito alcune possibilità di applicazione allo studio del comportamento di circuiti elettronici.

Per tale studio sarebbe normalmente indispensabile un prototipo del cicuito ed un’adatta attrezzatura di laboratorio (multimetri, oscilloscopi, generatori di segnali, strumenti digitali, ecc.) ma, se si può esprimere matematicamente sia la forma dei segnali d’ingresso (o di ‘stimolo’, come vengono chiamati negli esperimenti), sia il comportamento dei singoli circuiti in studio, è più semplice ed immediato ricorrere al calcolo diretto dei segnali risultanti.

Questa forma di simulazione dell’esperimento ha indubbiamente il vantaggio di una grande semplicità sia formale che pratica.

Qualsiasi intervento correttivo non solo nei parametri in gioco, ma anche nelle forme dei segnali di stimolo e nella stessa struttura dei circuiti può essere ottenuto con una semplice (almeno in senso formale) intervento sui dati di calcolo.

Ovviamente ciò è merito anche dell’ambiente di sviluppo in cui si opera e sulla conoscenza dell’operatore nell’utilizzo di tali mezzi.

Il  Mathcadâ è in tal senso fra i più completi e flessibili ambienti di sviluppo matematico, e le strutture  di impostazione dei calcoli sono certamente fra le più intuitive. Non occorre però sottovalutare che per un suo uso corretto, anche questa forma di ‘programmazione’ deve essere studiata ed assimilata.

Come primo semplice esempio di simulazione di una risposta di un circuito, del tipo  RC  già visto nel  capitolo 7,  ad uno stimolo costituito da  un segnale a semionda raddrizzata a 50 Hz (con o senza parzializzazione), vengono dati nella  Fig. 13.1 tre diversi casi.

In ciascuno si può notare:

- la generazione di  N  campioni del segnale d’ingresso (x_n  con  n = 0...N-1 )

- la trasformazione in  spettro  X(f)  di  K  frequenze (con  k = 0...N/2-1)

- la moltiplicazione per la funzione di trasferimento  H(f)  (nell’attuale caso di una costante di tempo, quindi   = 1/(1+j·w·T), con  T  posto uguale a 20 ms) per ottenere  lo spettro  Y_k   del segnale d’uscita .

- la  riconversione di quest’ultimo in segnale nel tempo  y_n.

Tutte queste operazioni sono svolte formalmente da una sola riga  di impostazioni, immediatamente sopra il  grafico che riporta l’andamento dei segnali d’ingresso e d’uscita per un tempo di 1/10 di secondo, come potrebbe essere rilevato da un oscilloscopio nel caso di reale esperimento in laboratorio.

Il calcolo è effettuato su un periodo di osservazione P di  0.1 secondi (100 ms) con  1024 campionamenti ( quindi 1 campionamento ogni circa 0.1 ms) e comprende perciò 10 semionde

Fig. 13.1   -   Esempi di analisi con FFT della risposta di circuiti a vari segnali d’ingresso (semionde a varia parzializzazione). [commento audio]

La precauzione di osservare il fenomeno su un periodo di tempo relativamente lungo non ha però fondamento, in quanto il metodo FFT fornisce direttamente risultati  a regime stazionario (cioè esenti da transitori iniziali).

Il calcolo ha come scopo la determinazione del valore medio  del segnale in uscita, cioè la media dei valori  di  y_n , e della sua ondulazione residua  (ripple) in corrispondenza ai vari gradi di parzializzazione (0,  25  e  50 % rispettivamente).

Tali calcoli sono svolti immediatamente sotto i rispettivi grafici.

Come si vede, una volta impostato il calcolo ed i singoli parametri, i risultati sia grafici che numerici vengono immediatamente visualizzati, con evidenti vantaggi rispetto ad un’analoga procedura sperimentale.

Un secondo esempio di applicazione diretta di FFT all’analisi circuitale è dato dalla  Fig. 13.2, in cui viene studiato il comportamento di un circuito differenziatore di secondo ordine all’applicazione di una forma d’onda trapezoidale. 

Fig. 13.2   -    Altro esempio di analisi circuitale con FFT:  forma d’onda trapezia attraverso un circuito derivatore di secondo ordine.

Il circuito derivatore è rappresentato dalla sua funzione di trasferimento H(s), dove s  è  espressione delle frequenze multiple della fondamentale (il cui periodo coincide col tempo di osservazione P), e  T  è la costante di derivazione.

La semplice derivazione sarebbe  s ·T , ma è opportuna una limitazione alle alte frequenze con una costante di tempo distanziata di una decade. L’elevazione al quadrato rende il derivatore di secondo ordine.

Quindi  x  costituisce il segnale d’ingresso ed  y il segnale d’uscita, ricavato calcolando lo spettro X,  moltiplicandolo per la funzione di trasferimento  H  ed antitrasformando il prodotto  Y  così ottenuto.

I grafici di  x  e  y  in funzione del tempo  t  (in ms) sostituiscono  brillantemente gli oscillogrammi che si sarebbero ricavati da una prova di laboratorio con circuiti reali

La  Fig. 13.3  è infine un esempio di configurazione circuitale a caratteristiche  non-lineari, facilmente simulabile mediante procedure di calcolo. 

Fig.13.3   -   Configurazione circuitale con caratteristiche non-lineari.

La procedura di calcolo ed i relativi risultati sono riportati in  Fig. 13.4, nell’ipotesi di una forma d’onda d’ingresso costituita da un segnale sinusoidale  a 10 Hz per 0.5 sec, seguito da una pausa per altri 0.5 sec.

Il periodo d’osservazione è quindi  P=1 sec, mentre il tempo  t_n  è espresso in ms.

Come chiaramente indicato nella  Fig.13.3, il circuito è composto da un blocco raddrizzatore (uscita  a_n)  e da una funzione di trasferimento  H(s)  costituita da una semplice costante di tempo con  T=0.2 sec  e da un fattore di amplificazione di 1.5, da cui esce il segnale  b_n.

La trasformazione in frequenza di a_n consente di ottenere  b_n semplicemente moltiplicando lo spettro A del segnale per la funzione di trasferimento H ed antitrasformando poi da frequenza in tempo.

Un terzo blocco limita il segnale d’ingresso x_n  fra  +0.5  e -0.2, generando il segnale  c_n.

Il segnale d’uscita  y_n  è infine il prodotto fra i segnali  b_n  e   c_n .

Anche nei casi di circuiti non-lineari (salvo ovviamente il blocco di trasferimento), il metodo di simulazione con l’uso di FFT dà quindi  possibilità di facili soluzioni, senza dover ricorrere a sperimentazioni in laboratorio

Fig. 13.4   -   Esempio di analisi in una configurazione circuitale non-lineare