Trasformata discreta di Fourier, relazione con la trasformata zeta

Prev Sezione 4.4: Trasformata di Fourier di sequenze Su Capitolo 4: Campionamento quantizzazione ed elaborazione numerica Sezione 4.6: Filtraggio numerico via DFT Next

4.5 Trasformata discreta di Fourier

Disponendo di una sequenza di lunghezza finita composta da N valori x_n, n = 0, 1, ..., N − 1, corrispondenti a campioni di un segnale x(t) prelevati con ritmo f_c = 1 T_c , si indica come Discrete Fourier Transform (DFT) la nuova sequenza [174] [174] La (10.75) può essere fatta discendere dalla (10.73) vincolando f ad assumere i soli valori discreti f = mN 1 T_c , e limitando l’indice della sommatoria ad un insieme finito di campioni, vedi fig. 4.24.

(10.75) X_m = ^N− 1⎲⎳_n = 0x_n e^{−j2π m N n}

univocamente definita per m = 0, 1, ..., N − 1, e che costituisce una approssimazione [175] [175] Una prima fonte di approssimazione deriva dall’operazione di finestratura legata all’uso di un numero finito di campioni, operando quindi su x_w(t) = x(t)w(t_c) anziché su x(t). Per analizzare le altre fonti di approssimazione, iniziamo a scrivere l’espressione di X_w(f) = F {x_w(t)} per f = m N f_c:

X_w⎛⎝f = m N f_c⎞⎠ = ∫^{(N − 1)T_c}₀x(t)e^{−j2πmN f_ct}dt ≃ ∑^N− 1_n = 0x_n ⋅ ∫^{(N − 1)T_c}₀ sinc(f_c(t − nT_c))e^{−j2πnN f_ct}dt

in cui la seconda eguaglianza utilizza l’interpolazione cardinale x(t) = ∑^∞_{n = −∞}x_n ⋅ sinc(f_c(t − nT_c)) fornita dalla (10.66), 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 f_c) di sinc(f_c(t − nT_c)), pari quest’ultima a T_c rect_{f_c}(f)e^−j2πfnT_c, che quando valutata per f = m N f_c, fornisce il risultato

X_w⎛⎝f = m N f_c⎞⎠ ≃ T_c∑^N− 1_n = 0x_ne^−j2πmNn

per valori |m| ≤ N 2 , a causa della estensione limitata (in frequenza) di rect_{f_c}(f). E’ però facile verificare che X_w⎛⎝m N f_c⎞⎠ è 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) = F {x(t)}, calcolata per f = m N f_c, e moltiplicata per f_c:

(10.76) X_m ≃ f_c X⎛⎝f = m N f_c⎞⎠

Notiamo subito che la (10.75) è valida per qualsiasi m, ed ha un andamento periodico con periodo N, a cui corrisponde la frequenza f_c = 1 T_c (vedi fig. 4.24), in accordo con la periodicità in frequenza evidenziata per X^•(f) (vedi (10.73) e (10.67)); per questo motivo, qualora il segnale originario x(t) contenga componenti a frequenze maggiori di f_c 2 , gli X_m con indici prossimi ad N2 presenteranno errore di aliasing [176] [176] Come osservato al § 4.1.1, lo spettro X^•(f) di un segnale campionato a frequenza f_c è costituito dalle repliche del segnale originario, distanziate di multipli di f_c: X^•(f) = ∑^∞_{n = −∞}X(f − nf_c), e coincide con X(f) per − f_c ⁄ 2 < f < f_c ⁄ 2, se X(f) è limitata in banda tra ± W ed f_c ≥ 2W. Al contrario, se f_c < 2W, allora le repliche X(f − nf_c) si sovrappongono, e la (10.76) si riscrive come X_m ≃ f_cX^•⎛⎝f = m N f_c⎞⎠.. Notiamo inoltre come i valori X_m siano tutti relativi a frequenze ≥ 0, ma nel caso di una sequenza x_n reale la proprietà di simmetria coniugata, associata alla periodicità in frequenza, rende il risultato interessante solo per indici m ≤ N2, dato che successivamente si trovano valori coniugati a quelli della prima metà.

Figure 4.24 Valori della DFT e campioni di X^•(f) associati per f = n N f_c

DFT come prodotto matriciale

Notiamo infine che la (10.75) può essere espressa in forma matriciale: ad esempio, per N = 4 si ottiene

(10.77)
⎡⎢⎢⎢⎢⎢⎣ X₀ X₁ X₂ X₃ ⎤⎥⎥⎥⎥⎥⎦ = ⎡⎢⎢⎢⎢⎢⎢⎢⎣ 1 1 1 1 1 e^{− j π 2} e^−jπ e^−j3π2 1 e^−jπ e^−j2π e^−j3π 1 e^{−j 3π 2} e^−j3π e^−j9π2 ⎤⎥⎥⎥⎥⎥⎥⎥⎦⎡⎢⎢⎢⎢⎢⎣ x₀ x₁ x₂ x₃ ⎤⎥⎥⎥⎥⎥⎦

da cui notiamo la proprietà di simmetria per la matrice dei coefficienti.

Esempio Allo scopo di concretizzare le differenze tra la trasformata di Fourier ed i valori forniti dalla dft, in fig. 4.25-a sono riportati i valori |X_m|, normalizzati in ampiezza, per la dft di una sinusoide a 10 Hz, adottando due diverse finestre di analisi (vedi § 3.8.4), 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 [177] [177] 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.25-b) con quello ottenibile per via analitica calcolando la F -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 F {} mediante la dft migliora, anche se persiste un ridotto potere di risoluzione spettrale.

a) modulo della DFT di una sinusoide a 10 Hz, f_c = 100 Hz, con finestre di diversa lunghezza

b - modulo della F -trasformata di una sinusoide a 10 Hz, con finestre di diversa lunghezza

Figure 4.25 Confronto tra DFT ed F -trasformata con uguale estensione temporale

Osservazione Qualche lettore potrebbe stupirsi di non trovare due linee spettrali, come da attendersi per una sinusoide. Ma questo 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 = ^T⁄_{T_c} che ricadono entro uno stesso periodo T della sinusoide [178] [178] Con la ovvia condizione che sia M > 2 per rispettare il vincolo f_c > ²⁄_T, con k che esprime quanti periodi entrano in N campioni, ed M quanti campioni/periodo. In tal caso l’applicazione della idft (10.78) produce una sequenza ancora periodica. Nel nostro esempio, essendo T = ¹⁄₁₀ = 100 msec e scegliendo N = 64 punti e k = 6 periodi, il vincolo N = kM = kTf_c permette di ottenere f_c = N kT = 64 6 ⋅ 10^− 1 = 106, 6 Hz, ovvero M = Tf_c = 10^− 1 ⋅ 106, 6 = 10, 6 campioni/periodo. La fig. 4.26 mostra questo risultato, evidenziando come la riga spettrale si manifesti per m = 6, ossia alla frequenza f = mN f_c = 6 64 106.6 = 10 Hz, mentre la riga presente in m = 58 è in realtà il ripiegamento periodico di quella a − 10 Hz.

Figure 4.26 Sinusoide campionata su di un numero intero di periodi e relativo modulo di dft

Il passaggio dai campioni nel tempo x_n a quelli in frequenza X_m è invertibile [179] [179] Sostituendo infatti la (10.75) nella (10.78), otteniamo

x_n = 1 N ∑^N− 1_m = 0(∑^N− 1_k = 0x_k e^−j2πmNk)e^j2πmNn = 1 N ∑^N− 1_k = 0x_k∑^N− 1_m = 0 e^{j2πmN(n − k)}

ma, dato che ∑^N− 1_m = 0 e^{j2πmN(n − k)} = ⎧⎨⎩ N se k = n 0 altrimenti , nella sommatoria esterna sopravvive solo il termine x_n, dimostrando l’uguaglianza.

, ricorrendo alla Inverse Discrete Fourier Transform (IDFT)

(10.78) x_n = 1 N ^N− 1⎲⎳_m = 0X_m e^{j2π 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 X_m rappresentano una approssimazione [180] [180] La relazione (10.79) si dimostra combinando le relazioni (10.34) e (10.76):

X_n ≃ f_cX⎛⎝nN f_c⎞⎠ = f_cX⎛⎝nNT_c ⎞⎠ = f_cX⎛⎝nT⎞⎠ = f_cX(nF) = f_cTX^SF_n = 1 T_c TX^SF_n = NX^SF_n

dei rispettivi coefficienti della serie di Fourier X^SF_m = 1 T ∫^T2_−T2x_T(t)e^{−j2πm T t}dt, calcolati a partire da un segmento x_T(t) estratto da x(t), e moltiplicati per N:

(10.79) X_m ≃ N ⋅ X^SF_m

Per approfondire i risvolti di questo risultato, affrontiamo la sezione successiva.

4.5.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 X(z) = Z{x_n}, funzione complessa della variabile complessa zeta, e valutata come

(10.80) X(z) =^∞⎲⎳_{n = −∞}x_n z^− n

che, nel caso in cui per |z| = 1 la (10.80) converga, può essere fatta corrispondere alla dtft (10.73) della stessa sequenza x_n semplicemente ponendo z = e^jω, ovvero calcolando X(z) sul cerchio unitario:

(10.81) X( e^jω) = ^∞⎲⎳_{n = −∞}x_n e^−jωn = X(z)|_{z = e^jω} = X^•(f)|_{f = ω 2πT_c}

Infatti, nelle consuete condizioni in cui gli x_n sono i campioni di un segnale x(t) limitato in banda tra ± W e prelevati con ritmo f_c ≥ 2W, la (10.81) effettivamente coincide (per − π ≤ ω < π) con la X^•(f) (eq. (10.73), per − ^f_c⁄₂ ≤ f ≤ ^f_c⁄₂) in cui si ponga f = ω2πT_c , mettendo cioè in corrispondenza le frequenze ± ^f_c⁄₂ di X^•(f) con le pulsazioni ± π di X(e^jω). Al di fuori dell’intervallo − π ≤ ω < π la X(e^jω) è periodica in ω con periodo 2π, analogamente a ciò che risulta (con periodo f_c) per la trasformata di Fourier X^•(f) di sequenze; se invece x_n è sottocampionata, ossia f_c < 2W, anche X(e^jω) è affetta da aliasing, così come avviene per X^•(f).

Figure 4.28 Sequenza x_n = aⁿ e modulo della relativa trasformata di Fourier a tempo discreto

Esempio Consideriamo la sequenza x_n = ⎧⎨⎩ aⁿ se n ≥ 0 0 altrim. il cui andamento per a = 0.7 è mostrato in fig. 4.28, la cui trasformata zeta X(z) = ∑^∞_{n = −∞}aⁿ z^− n risulta pari a [181] [181] Il risultato si ottiene ricordando che ∑^∞_n = 0αⁿ = 1 1 − α qualora |α| < 1. X(z) = 1 1 − az^− 1 , ed il cui modulo, dopo aver scritto la variabile complessa z come z = x + jy, è espresso come |X(z)| = ¹⁄_{√ ⎛⎝x² − ax + y² x² + y² ⎞⎠² + ⎛⎝ayx² + y² ⎞⎠²}.

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 la forma di

(10.82) |X(e^jω)| = 1√1 + a² − 2a cosω

ossia della dtft della sequenza aⁿ, che a sua volta è mostrata in fig. 4.28 per − π < ω < π, e nella figura a lato come una linea rossa.

Se la X(e^jω) ottenuta per una sequenza x_n aperiodica nel tempo è campionata in N punti equispaziati e disposti sul cerchio unitario, ossia ponendo ω = 2π m N con m = 0, 1, …, N − 1, allora si ottiene una sequenza periodica in frequenza [182] [182] Infatti e^{−j2πm + N N n} = e^{−j2π m N n}e^−j2πn, ed il secondo termine vale 1 per qualsiasi n. Indichiamo qui ed al prossimo §, una sequenza periodica mediante la tilde .̃.

(10.83) X̃_m = ^∞⎲⎳_{n = −∞}x_n e^{−j2π m N n} = X( e^jω)|_{ω = 2π m N} = X^•(f)|_{f = f_c m N}

che può coincidere con la sequenza X_m ottenuta calcolando la dft (10.75) di una sequenza x_n, qualora questa abbia una durata limitata ≤ N. D’altro canto, è possibile applicare la idft (10.78) ad un periodo della sequenza X̃_m, ed ottenere quindi una nuova sequenza di valori nel tempo, anch’essa periodica di periodo N, espressa come

(10.84) x̃_n = 1 N ^N− 1⎲⎳_m = 0X̃_m e^{j2π m N n}

Infatti, i valori x̃_n dipendono da quelli x_n = x(t)|_{t = nT_c} del segnale originario x(t), campionato agli istanti t = nT_c, mediante la relazione [183] [183] Infatti, sostituendo la (10.83) in (10.84), otteniamo x̃_n = 1 N ∑^N− 1_k = 0 ∑^∞_{h = −∞}x_h e^{−j2π k N h}e^{j2πk N n}. Scambiando ora l’ordine delle sommatorie risulta

x̃_n = ∑^∞_{h = −∞}x_h⎛⎝1 N ∑^N− 1_k = 0 e^{−j2πkN(h − n)}⎞⎠

Dato che 1 N ∑^N− 1_k = 0 e^{−j2πkN(h − n)} = ⎧⎨⎩ 1 se h = n + rN 0 altrimenti , con r intero, si ottiene il risultato (10.85).

(10.85) x̃_n = ^∞⎲⎳_{r = −∞}x_n + rN

e quindi i primi N valori di x̃_n coincidono con i campioni di x(t) solo se quest’ultimo ha durata limitata, con estensione minore di NT_c, ossia se N è sufficientemente elevato in modo che NT_c copra tutta la durata di x(t), e la (10.80) si riconduca alla somma di un numero finito di termini. D’altra parte, se x(t) ha durata maggiore di NT_c, ovvero X(z) è stata campionata su di un numero di campioni troppo ristretto, allora l’applicazione della IDFT (10.84) ad X̃(k) provoca il fenomeno di aliasing temporale.

Esempio Nella parte sinistra di fig. 4.29 viene mostrato il modulo della sequenza X_m ottenuta campionando la X(e^jω) dell’esempio precedente, utilizzando 16, 8 o 4 campioni/periodo. Nella parte destra della stessa figura sono quindi rappresentate le corrispondenti sequenze x_n ottenute mediante idft. Si può notare che, mentre con 16 campioni/periodo la ricostruzione della sequenza x_n = aⁿ è piuttosto fedele, con 8 campioni si inizia a verificare il fenomeno di aliasing temporale, che diviene ancor più evidente per 4 campioni/periodo.

Figure 4.29 Aliasing temporale al diminuire della risoluzione del campionamento in frequenza

Riepilogo

La figura 4.30 riassume le relazioni che legano la trasformata di Fourier per segnali limitati in banda ai suoi campioni ed alla relativa dtft, così come viene illustrata la relazione tra dtft, dft, e trasformata zeta.

Figure 4.30 Grafo di corrispondenza tra trasformate continue e discrete

4.5.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 X_m della (10.75), considerando i valori e^{−j2πm N n} precalcolati (vedi (10.77)), richieda N moltiplicazioni complesse (equivalenti ognuna a 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) ≃ 8N²

operazioni.

Al contrario, gli algoritmi fft più efficienti riducono il numero di operazioni ad 8Nlog₂N: ad esempio, ponendo N = 1024, si ottiene un miglioramento di ^{2³(2¹⁰)²}⁄_{2³2¹⁰ ⋅ 10} = ^2¹⁰⁄₁₀ ≃ 100 volte! Queste prestazioni sono legate all’adozione di un valore di N che sia una potenza di due (ossia N = 2^Mcon M intero), ma successivamente sono stati individuati metodi [184] [184] 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.5.3 Relazione tra DFT e DCT

Anche per la dft risulta valida la proprietà di simmetria coniugata (§ 84) e quindi, se i valori della sequenza x_n di lunghezza N che compare nella (10.75) sono reali anziché complessi, allora i coefficienti di dft X_m presentano parte reale pari e parte immaginaria dispari. In particolare, se immaginiamo di estendere la lunghezza della sequenza a 2N punti, ottenuti ribaltando sugli indici negativi la sequenza di partenza come x_− n = x_n (vedi prima riga di fig. 4.31), allora siamo nelle condizioni di sequenza reale pari, che

Figure 4.31 Estensione pari di sequenza reale

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 al passo precedente in modo da renderla effettivamente pari (seconda riga di fig. 4.31). Considerando che per segnali reali pari le componenti sinusoidali della base della dft non danno contributi al risultato [185] [185] Scriviamo la (10.75) come

X_m = ^{N − ¹⁄₂}⎲⎳_{n’ = − N + ¹⁄₂}x_{n’ − ¹⁄₂} e^{−j2πn’2Nm} = ^{N − ¹⁄₂}⎲⎳_{n’ = − N + ¹⁄₂}x_{n’ − ¹⁄₂}cos⎛⎝2πn’2Nm⎞⎠ − j^{N − ¹⁄₂}⎲⎳_{m’ = − N + ¹⁄₂}x_{n’ − ¹⁄₂}sin⎛⎝2πn’2Nm⎞⎠ = = 2^{N − ¹⁄₂}⎲⎳_{n’ = ¹⁄₂}x_{n’ − ¹⁄₂}cos⎛⎝2πn’2Nm⎞⎠ = 2^N− 1⎲⎳_n = 0x_ncos⎛⎝2πn + ¹⁄₂ 2N m⎞⎠ = 2^N− 1⎲⎳_n = 0x_ncos⎡⎣πN ⎛⎝n + 12⎞⎠m⎤⎦

in cui x_n’ è quella disegnata per seconda in fig. 4.31. 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 x_{n’ − ¹⁄₂} pari e sin⎛⎝2π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’ − ¹⁄₂, mentre l’ultima è (a parte il fattore 2) la definizione della DCT data in (10.86). , e adottando un nuovo cambio di variabile, si ottiene in definitiva la formula di calcolo della dct come

(10.86) X_m = ⎲⎳^N− 1_n = 0x_ncos⎡⎣π N ⎛⎝n + 1 2 ⎞⎠m⎤⎦

a cui è associata la trasformazione inversa idct

(10.87)
x_n = 12 X ₀ + ⎲⎳^N− 1_m = 1 X_m cos⎡⎣π N ⎛⎝m + 1 2 ⎞⎠n⎤⎦

La dct verrà usata nell’ambito della compressione di immagini (§ 304). Infatti, i valori di luminanza dei pixel in cui si scompone una immagine sono tutti valori reali.

4.5.4 DFT come un banco di filtri

Un banco di filtri è qualcosa che esegue simultaneamente diverse operazioni di filtraggio sullo stesso segnale, emettendo tutti i risultati contemporaneamente. I valori

X_m = ⎲⎳^N− 1_n = 0 x_n e^{−j2π m N n}

di una dft ad N punti per una sequenza x_n equivalgono all’output decimato di N filtri, ciascuno con risposta impulsiva h_m(n) = e^{j2π N m ⋅ n} con m, n = 0, 1, ⋯, N − 1, rappresentate nella figura a lato:

viene mostrata solo la parte reale di h_m(n);
la decimazione avviene con periodo pari alla separazione delle finestre di segnale (sovrapposte o meno);
ogni filtro h_m ha il modulo della risposta in frequenza approssimativamente pari ad un sinc (dovuto ad una finestra rettangolare) centrato alla frequenza digitale 2πmN (vedi figura seguente);
un progetto più specifico del filtro può ottenere una risposta in frequenza più selettiva, vedi cap. 5.

Un’analisi più approfondita di ciò che accade può essere fatta partendo dal filtro che produce X₀ e che è un filtro a media mobile (pagina 1), la cui uscita è esattamente uguale a X₀ = ∑^N− 1_n = 0x_n e^{−j2π 0 N n} = ∑^N− 1_n = 0x_n, che è la media dei valori dei campioni: pertanto il filtro risultante ha una risposta in frequenza passa-basso della forma mostrata a pagina 1, con N − 1 zeri disposti sul cerchio unitario. Per gli altri elementi della dft con m ≠ 0 la risposta impulsiva a media mobile viene modulata moltiplicandola per e^−j2π^{m N n}, in modo che la sua risposta in frequenza trasli verso l’indice m, corrispondente alla frequenza digitale ω = 2π m N . In definitiva, la risposta in frequenza composita viene mostrata sotto [186] [186] Immagine presa da https://wirelesspi.com/discrete-fourier-transform-dft-as-a-filter-bank/ con una scala verticale in dB, per il caso di una dft con 12 punti.

Prev Sezione 4.4: Trasformata di Fourier di sequenze Su Capitolo 4: Campionamento quantizzazione ed elaborazione numerica Sezione 4.6: Filtraggio numerico via DFT Next