Sezione 4.3: Trasformata di Fourier di sequenze Su  Capitolo 4: Campionamento, quantizzazione ed elaborazione numerica Sezione 4.5: Filtraggio numerico via DFT 

4.4  Trasformata discreta di Fourier

Disponendo di una sequenza di lunghezza finita composta da N valori xn, n = 0, 1, ..., N − 1, corrispondenti a campioni di un segnale x(t) prelevati con ritmo fc  = (1)/(Tc), si indica come Discrete Fourier Transform (DFT) la nuova sequenza[105] [105] La (8.28↓) può essere fatta discendere dalla (8.26↑) vincolando f ad assumere i soli valori discreti f = (m)/(N)(1)/(Tc), e limitando l’indice della sommatoria ad un insieme finito di campioni.
(8.28) Xm  = N − 1n = 0xne − j2π(m)/(N)n
univocamente definita per m = 0, 1, ..., N  − 1, e che costituisce una approssimazione[106]  [106]  Una prima fonte di approssimazione deriva dall’operazione di finestratura legata all’uso di un numero finito di campioni, operando quindi su xw(t) =  x(t)w(tc) anziché su x(t). Per analizzare le altre fonti di approssimazione, iniziamo a scrivere l’espressione di Xw(f) =  ℱ{xw(t)} per f = (m)/(N)fc:
Xwf  = (m)/(N)fc  =  (N  − 1)Tc0x(t)e − j2π(m)/(N)fctdt N − 1n = 0xn(N  − 1)Tc0 sinc(fc(t  − nTc))e − j2π(n)/(N)fctdt
in cui la seconda eguaglianza utilizza l’interpolazione cardinale x(t) = n  =  − ∞xnsinc(fc(t  − nTc)) fornita dalla (8.21↑), ed introduce una seconda fonte di approssimazione legata all’intervallo finito di variazione per n: infatti, benché l’integrale abbia estensione limitata, i valori di x(t) che cadono entro tale estensione, dovrebbero dipendere da tutti i suoi campioni. L’ultimo integrale è a sua volta una approssimazione (a causa degli estremi di integrazione limitati, e peggiore per i sinc centrati in prossimità dei confini della finestra) della trasformata (calcolata in f = (m)/(N)fc) di sinc(fc(t − nTc)), pari quest’ultima a Tcrectfc(f)e  − j2πfnTc, che quando valutata per f = (m)/(N)fc, fornisce il risultato
Xwf  =  (m)/(N)fcTcN − 1n = 0xne  − j2π(m)/(N)n
per valori |m| ≤ (N)/(2), a causa della estensione limitata (in frequenza) di rectfc(f). E’ però facile verificare che Xw(m)/(N)fc è periodica in m con periodo N, cosicché i valori assunti per m = (N)/(2)  + 1, (N)/(2)  + 2, … sono uguali a quelli per m  =  − (N)/(2)  + 1,  − (N)/(2)  + 2, .
del campionamento in frequenza della trasformata X(f) = ℱ{x(t)}, calcolata per f = (m)/(N)fc, e moltiplicata per fc:
(8.29) XmfcXf  = (m)/(N)fc
Notiamo subito che la (8.28↑) è valida per qualsiasi m, ed ha un andamento periodico con periodo N, a cui corrisponde la frequenza fc  = (1)/(Tc), in accordo con la periodicità in frequenza evidenziata per X(f) (vedi (8.26↑) e fig. 4.2↑); per questo motivo, qualora il segnale originario x(t) contenga componenti a frequenze maggiori di (fc)/(2), gli Xm con indici prossimi ad (N)/(2) presenteranno errore di aliasing[107] [107] Come osservato al § 4.1.1↑, lo spettro X(f) di un segnale campionato a frequenza fc è costituito dalle repliche del segnale originario, distanziate di multipli di fc: X(f) = n  =  − ∞X(f  − nfc), e coincide con X(f) per  − fc  ⁄ 2 < f < fc ⁄ 2, se X(f) è limitata in banda tra ±W ed fc ≥ 2W. Al contrario, se fc  < 2W, allora le repliche X(f − nfc) si sovrappongono, e la (8.29↑) si riscrive come XmfcXf  = (m)/(N)fc.. Notiamo inoltre come i valori Xm siano tutti relativi a frequenze  ≥ 0, ma nel caso di una sequenza xn reale la proprietà di simmetria coniugata, associata alla periodicità in frequenza, rende il risultato interessante solo per indici m ≤ (N)/(2), dato che successivamente si trovano valori coniugati a quelli della prima metà. Notiamo infine che la (8.28↑) può essere espressa in forma matriciale: ad esempio, per N  = 4 si ottiene
(8.30) X0 X1 X2 X3   =  1 1 1 1 1 e − j(π)/(2) e  − jπ e − j(3π)/(2) 1  e  − jπ e − j2π e  − j3π 1  e − j(3π)/(2)  e − j3π e  − j(9π)/(2) x0 x1 x2 x3
da cui notiamo la proprietà di simmetria per la matrice dei coefficienti.
EsempioAllo scopo di concretizzare le differenze tra la trasformata di Fourier ed i valori forniti dalla dft, in fig. 4.18↑-a sono riportati i valori |Xm|, normalizzati in ampiezza, per la dft di una sinusoide a 10 Hz, adottando due diverse finestre di analisi (vedi § 3.9.3↑), prelevando alla medesima frequenza di campionamento (100 Hz) un numero variabile di campioni (mostrato in figura), e ponendo i rimanenti a zero, per calcolare in tutti i casi la medesima DFT a 256 punti[108]  [108] Il metodo esposto di porre a zero i campioni fino al raggiungimento di una potenza di due è detto zero padding. Il calcolo della DFT su di un numero di punti pari ai campioni di segnale disponibili, non avrebbe dato luogo all’effetto finestra, ma avrebbe fornito in tutti i casi andamenti simili a quello osservabile per 256 punti. Infine, notiamo che nelle figure sono mostrati solo i primi 128 valori, essendo i rimanenti speculari.. Il risultato è quindi confrontato (fig. 4.18↑-b) con quello ottenibile per via analitica calcolando la -trasformata dello stesso segnale, adottando le medesime finestre temporali, di durata uguale al primo caso. Le curve ottenute nel caso di 80 msec (e 8 campioni!) dipendono da meno di un periodo di segnale, e perciò presentano una componente continua apprezzabile. Aumentando la durata della finestra, l’approssimazione di calcolare una {} mediante la dft migliora, anche se persiste un ridotto potere di risoluzione spettrale.


a - modulo della DFT di una sinusoide a 10 Hz, fc  =  100 Hz, con finestre di diversa lunghezza
Trasformata discreta di Fourier Trasformata discreta di Fourier
b - modulo della -trasformata di una sinusoide a 10 Hz, con finestre di diversa lunghezza
Trasformata discreta di Fourier Trasformata discreta di Fourier

Figura 4.18 Confronto tra DFT ed -trasformata con uguale estensione temporale
OsservazioneProbabilmente a questo punto qualche lettore può trovarsi stupito di non avere incontrato due linee spettrali, come sarebbe stato lecito aspettarsi per una sinusoide. In realtà ciò può accadere, a patto che il numero di campioni N su cui si effettua la dft sia un multiplo intero k del numero di campioni M  = TTc che ricadono entro uno stesso periodo T della sinusoide[109] [109] Con la ovvia condizione che sia M > 2 per rispettare il vincolo fc  > 2T, ovvero k esprime quanti periodi entrano in N campioni, ed M quanti campioni/periodo: infatti in questo caso, l’applicazione della idft (8.31↓) produce una sequenza ancora periodica. Nel nostro esempio, essendo T  = 110  = 100 msec e scegliendo N = 64 punti e k = 6 periodi, il vincolo N  = kM = kTfc permette di ottenere fc  = (N)/(kT) = (64)/(6⋅10  − 1) = 106, 6 Hz, ovvero M  = Tfc = 10 − 1⋅106, 6  = 10, 6 campioni/periodo. La fig. 4.19↓ mostra questo risultato, evidenziando come la riga spettrale si manifesti per m = 6, ossia alla frequenza f  = (m)/(N)fc  = (6)/(64)106.6 = 10 Hz, mentre la riga presente in m = 58 è in realtà il ripiegamento periodico di quella a   − 10 Hz.
Discrete fourier trasform di sinusoide Discrete fourier trasform di sinusoide
Figura 4.19 Sinusoide campionata su di un numero intero di periodi e relativo modulo di dft
Il passaggio dai campioni nel tempo xn a quelli in frequenza Xm è invertibile[110] [110] Sostituendo infatti la (8.28↑) nella (8.31↓), otteniamo
xn = (1)/(N)N  − 1m = 0(N  − 1k = 0xke − j2π(m)/(N)k)ej2π(m)/(N)n = (1)/(N)N  − 1k = 0xkN − 1m = 0 ej2π(m)/(N)(n  − k)
ma, dato che N − 1m  = 0 ej2π(m)/(N)(n  − k) =  N sek = n 0  altrimenti , allora nella sommatoria esterna sopravvive solo il termine xn, dimostrando l’uguaglianza.
, ricorrendo alla Inverse Discrete Fourier Transform (IDFT)
(8.31) xn  = (1)/(N)N  − 1m = 0Xmej2π(m)/(N)n
che per n esterno a [0, N − 1] continua a valere, ed assume valori periodici, coerentemente a quanto accade per lo sviluppo in serie di Fourier. Infatti il legame tra dft e serie di Fourier è molto stretto, in quanto i valori Xm rappresentano una approssimazione[111] [111] La relazione (8.32↓) si dimostra combinando le relazioni (8.2↑) e (8.29↑): XnfcX(n)/(N)fc  = fcX(n)/(NTc) = fcX(n)/(T) = fcX(nF)  = fcTXSFn  = (1)/(Tc)TXSFn  = NXSFn dei rispettivi coefficienti della serie di Fourier XSFm  = (1)/(T)(T)/(2) − (T)/(2)xT(t)e  − j2π(m)/(T)tdt, calcolati a partire da un segmento xT(t) estratto da x(t), e moltiplicati per N:
(8.32) XmNXSFm
Per approfondire i risvolti di questo risultato, affrontiamo la sezione successiva.

4.4.1  Relazione tra DTFT, DFT e trasformata zeta

Così come per i segnali analogici sussiste una relazione (vedi pag. 1↓) tra la trasformata di Fourier e quella di Laplace, così nel contesto delle sequenze, esistono legami tra dtft e trasformata zeta, definita quest’ultima come
(8.33) X(z) = n  =  − ∞xnz  − n
che, nel caso in cui converga per |z| = 1, può essere fatta corrispondere alla dtft (8.26↑) semplicemente ponendo z  = ejω, ovvero calcolando X(z) sul cerchio unitario:
(8.34) X(ejω)  = n =  − ∞xne − jωn = X(z)|z  = ejω = X(f)|f = (ω)/(2πTc)
Infatti, nelle consuete condizioni in cui gli xn sono i campioni di un segnale x(t) limitato in banda tra ±W e prelevati con ritmo fc  ≥ 2W, la (8.34↑) effettivamente coincide (per  − π  ≤ ω < π) con la X(f) (eq. (8.26↑), per  − (fc)/(2)  ≤ f ≤ (fc)/(2)) in cui si ponga f = (ω)/(2πTc), mettendo cioè in corrispondenza le frequenze ±(fc)/(2) di X(f) con le pulsazioni ±π di X(ejω). Al di fuori di tale intervallo, X(ejω) è periodica in ω con periodo 2π, analogamente a ciò che risulta (con periodo fc) per la trasformata di Fourier X(f) di sequenze; se invece xn è sottocampionata, ossia fc  < 2W, anche X(ejω) è affetta da aliasing, così come avviene per X(f).
trasformata zeta di esponenziale trasformata zeta di esponenziale

                
Figura 4.20 Sequenza xn  = an e modulo della relativa trasformata di Fourier a tempo discreto
EsempioConsideriamo la sequenza xn  =  anse  ≥  0    0     altrim. il cui andamento per a =  0.7 è mostrato in fig. 4.20↓, la cui trasformata zeta X(z)  = n  =  − ∞anz  − n risulta pari a[112] X(z)  = 11 − az − 1, ed il cui modulo, dopo aver scritto la variabile complessa z come z  = x + jy,è espresso come |X(z)| = 1x2  − ax + y2x2  + y22  + ayx2  + y22.
figure f3.36.png
Figura 4.22 Estensione pari di sequenza reale
Facendo ora variare z nell’intervallo [  − 2 − j2,  2 + j2] si ottiene per il modulo di X(z) l’andamento mostrato nella figura a lato, in cui è anche raffigurato un cilindro di raggio unitario, la cui intersezione con |X(z)| individua l’andamento di |X(ejω)| = 11 + a2  − 2acosω, ossia la dtft della sequenza an, che a sua volta è mostrata in fig. 
4.20↓ per  − π < ω  < π.

Se la X(ejω) ottenuta per una sequenza xn aperiodica nel tempo è campionata in N punti equispaziati e disposti sul cerchio unitario, ossia ponendo ω = 2π(m)/(N) con m = 0, 1, …, N − 1, si ottiene allora una sequenza periodica in frequenza[113] [113] Infatti, e − j2π(m +  N)/(N)n = e  − j2π(m)/(N)ne − j2πn, ed il secondo termine vale 1 per qualsiasi n. Indichiamo qui ed al prossimo §, una sequenza periodica mediante la tilde .
(8.35) m  = n  =  − ∞xne  − j2π(m)/(N)n  = X(ejω)|ω  = 2π(m)/(N)  = X(f)|f  = fc(m)/(N)
che può coincidere con la sequenza Xm ottenuta calcolando la dft (8.28↑) di una sequenza xn, qualora questa abbia una durata limitata  ≤ N. D’altro canto, è possibile applicare la idft (8.31↑) ad un periodo della sequenza m, ed ottenere quindi una nuova sequenza di valori nel tempo, anch’essa periodica di periodo N, espressa come
(8.36) n  = (1)/(N)N  − 1m = 0mej2π(m)/(N)n
Infatti, i valori n dipendono da quelli xn  = x(t)|t  = nTc del segnale originario x(t), campionato agli istanti t = nTc, mediante la relazione[114] [114] Infatti, sostituendo la (8.35↑) in (8.36↑), otteniamo n = (1)/(N)N − 1k  = 0h  =  − ∞xhe  − j2π(k)/(N)hej2π(k)/(N)n. Scambiando ora l’ordine delle sommatorie risulta
n = h  =  − ∞xh(1)/(N)N  − 1k = 0  e − j2π(k)/(N)(h − n)
Dato che (1)/(N)N − 1k  = 0 e − j2π(k)/(N)(h  − n) =  1   seh = n + rN 0  altrimenti , con r intero, si ottiene il risultato (8.37↓).
(8.37) n  = r  =  − ∞xn  + rN
e quindi i primi N valori di n coincidono con i campioni di x(t) solo se quest’ultimo ha durata limitata, con estensione minore di NTc, ossia se N è sufficientemente elevato in modo che NTc copra tutta la durata di x(t), e la (8.33↑) si riconduca alla somma di un numero finito di termini. D’altra parte, se x(t) ha durata maggiore di NTc, ovvero X(z) è stata campionata su di un numero di campioni troppo ristretto, allora l’applicazione della IDFT (8.36↑) ad (k) provoca il fenomeno di aliasing temporale.
EsempioNella parte sinistra di fig. 4.21↓ viene mostrato il modulo della sequenza Xm ottenuta campionando la X(ejω) dell’esempio precedente, utilizzando 16, 8 o 4 campioni/periodo. Nella parte destra della stessa figura sono quindi rappresentate le corrispondenti sequenze xn ottenute mediante idft. Si può notare che, mentre con 16 campioni/periodo la ricostruzione della sequenza xn = an è piuttosto fedele, con 8 campioni si inizia a verificare il fenomeno di aliasing temporale, che diviene ancor più evidente per 4 campioni/periodo.
Aliasing temporale Aliasing temporale
Aliasing temporale
Aliasing temporale
Aliasing temporale
Aliasing temporale
Figura 4.21 Aliasing temporale al diminuire della risoluzione del campionamento in frequenza

4.4.2  Fast Fourier Transform

La sigla fft descrive una classe di algoritmi di calcolo della dft e della sua inversa, caratterizzati dall’uso di un numero molto ridotto di operazioni, rendendo così computazionalmente praticabile l’elaborazione numerica dei segnali. Analizziamo innanzitutto come il calcolo di ognuno degli N termini Xm della (8.28↑), considerando i valori e  − j2π(m)/(N)n precalcolati (vedi (8.30↑)), richieda N moltiplicazioni complesse (equivalenti ognuna ad 4 moltiplicazioni e 2 somme reali) ed N − 1 somme complesse (ognuna pari a 2 somme reali): pertanto una dft richiede N(N(4 + 2) + 2(N − 1)) = N(8N  − 2)≃8N2 operazioni. Al contrario, gli algoritmi fft più efficienti riducono il numero di operazioni ad 8Nlog2N: ad esempio, ponendo N = 1024, si ottiene un miglioramento di 23(210)223210⋅10  = 21010≃100 volte! Queste prestazioni sono legate all’adozione di un valore di N che sia una potenza di due (ossia N = 2Mcon M intero), ma successivamente sono stati individuati metodi[115] [115] Vedi ad es. http://it.wikipedia.org/wiki/Trasformata_di_Fourier_veloce che permettono una efficienza di calcolo comparabile anche per finestre di analisi di lunghezza qualsiasi.

4.4.3  Relazione tra DFT e DCT

Anche per la dft risulta valida la proprietà di simmetria coniugata (§ 46↑) e quindi, se i valori della sequenza xn di lunghezza N che compare nella (8.28↑) sono reali anziché complessi, allora i coefficienti di dft Xm presentano parte reale pari e parte immaginaria dispari. In particolare, se
discrete cosine transform
Figura 4.22 Estensione pari di sequenza reale
immaginiamo di estendere la lunghezza della sequenza a 2N punti, ottenuti ribaltando sugli indici negativi la sequenza di partenza come x − n = xn (vedi prima riga di fig. 4.22↑), allora siamo nelle condizioni di sequenza reale pari, che determina una trasformata solo reale (e pari), con parte immaginaria nulla.
Per arrivare a definire la Discrete Cosine Transform (dct) si calcola una dft bilatera sulla sequenza lunga 2N ottenuta traslando quella descritta in modo da renderla effettivamente pari (seconda riga di fig. 4.22↑). Considerando che per segnali reali pari le componenti sinusoidali della base della dft non danno contributi al risultato[116] [116] Scriviamo la (8.28↑) come
Xm  =  N  − 12n’ =   − N + 12xn’  − 12  e − j2π(n)/(2N)m =   =  N − 12n’ =   − N + 12xn’  − 12cos2π(n)/(2N)m  − jN  − 12m’ =   − N + 12xn’  − 12sin2π(n)/(2N)m  =   =  2N − 12n’ = 12xn’  − 12cos2π(n)/(2N)m  = 2N − 1n = 0xncos2π(n  + 12)/(2N)m  =   =  2N − 1n = 0xncos(π)/(N)n  + (1)/(2)m
in cui xn è quella disegnata per seconda in fig. 4.22↑. La quarta eguaglianza tiene conto del fatto che il termine immaginario si annulla, in quanto sommatoria bilatera di una funzione dispari (ottenuta come prodotto di xn’  − 12 pari e sin2π(n)/(2N)m dispari), e del fatto che essendo i termini coseno pari, la sommatoria può essere ristretta ai soli indici positivi, raddoppiati. La penultima eguaglianza rappresenta il semplice cambio di variabile n = n’ − 12, mentre l’ultima è (a parte il fattore 2) la definizione della DCT data in (8.38↓).
, e adottando un nuovo cambio di variabile, si ottiene in definitiva la formula di calcolo della dct come
(8.38) Xm  = N − 1n = 0xncos(π)/(N)n + (1)/(2)m
a cui è associata la trasformazione inversa idct
(8.39) xn  = (1)/(2)X0 + N − 1m = 1Xmcos(π)/(N)m + (1)/(2)n
La dct verrà usata in questo testo nell’ambito della compressione di immagini (§ 228↓): infatti i valori di luminanza dei pixel in cui si scompone una immagine, sono tutti valori reali.
  Sezione 4.3: Trasformata di Fourier di sequenze Su  Capitolo 4: Campionamento, quantizzazione ed elaborazione numerica Sezione 4.5: Filtraggio numerico via DFT 
x Logo

Trasmissione dei Segnali e Sistemi di Telecomunicazione

http://infocom.uniroma1.it/alef/libro/

Un esperimento divenuto nel tempo un riferimento culturale. Scopri come effettuare il download, ricevere gli aggiornamenti, e contribuire!