Elaborazione numerica dei segnali/DFT e FFT

Trasformata di Fourier discreta (DFT)

Definizione di DFT

Data una sequenza $x\left(n\right)$ costituita da un numero $N$ finito di campioni, la sua trasformata di Fourier discreta (DFT) $X\left(k\right)$, all'interno del suo periodo $N$, è definita:

$X\left(k\right)=\sum _{n=0}^{N-1}x\left(n\right)e^{-j2\pi n{\frac {k}{N}}},\quad k=0,1,2,\ldots ,N-1$

e può essere interpretata come la discretizzazione in frequenza della DTFT $X\left(e^{j2\pi f}\right)$, di cui vengono prese $N$ frequenze $f_{k}$ equispaziate:

$f_{k}={k \over N},\quad k=0,1,2,\ldots ,N-1$

Dal punto di vista algoritmico, la DFT è di complessità inferiore rispetto alla DTFT, e inoltre è definita in modo discreto → è rappresentabile su un calcolatore.

Inversione della DFT (IDFT)

L'antitrasformata della DFT, all'interno del suo periodo $N$, è definita:

$x\left(n\right)={\frac {1}{N}}\sum _{k=0}^{N-1}X\left(k\right)e^{j2\pi n{\frac {k}{N}}},\quad n=0,1,2,\ldots ,N-1$

Estensioni periodiche

La DFT e la IDFT sono periodiche di periodo $N$:

${\begin{cases}X\left(k\right)=X\left(k+N\right)\quad k=0,\ldots ,N-1\\x\left(n\right)=x\left(n+N\right)\quad n=0,\ldots ,N-1\end{cases}}$

Si definiscono le estensioni periodiche ${\overline {X}}\left(k\right)$ e ${\overline {x}}\left(n\right)$, rispettivamente di $X\left(k\right)$ e $x\left(n\right)$:

${\begin{cases}{\overline {X}}\left(k\right)=\sum _{n=0}^{N-1}x\left(n\right)e^{-j2\pi n{\frac {k}{N}}}\quad \forall k\\{\overline {x}}\left(n\right)={\frac {1}{N}}\sum _{k=0}^{N-1}X\left(k\right)e^{j2\pi n{\frac {k}{N}}}\quad \forall n\end{cases}}$

Tecnica di aggiunta degli zeri

È possibile teoricamente ricavare per interpolazione la DTFT partendo dagli $N$ campioni della DFT:

$X\left(e^{j2\pi f}\right)={\frac {1}{N}}\sum _{k=0}^{N-1}X\left(k\right){\frac {\sin {\left(\pi Nf-k\pi \right)}}{\sin {\left(\pi f-\pi {\frac {k}{N}}\right)}}}e^{-j\pi \left(N-1\right)\left(f-{\frac {k}{N}}\right)}$

Dimostrazione

$X\left(e^{j2\pi f}\right)=\sum _{n=-\infty }^{+\infty }x\left(n\right)e^{-j2\pi fn}=$

Siccome $x\left(n\right)$ ha per ipotesi durata finita $N$:

$=\sum _{n=0}^{N-1}x\left(n\right)e^{-j2\pi fn}=\sum _{n=0}^{N-1}\left({\frac {1}{N}}\sum _{k=0}^{N-1}X\left(k\right)\sum _{k=0}^{N-1}X\left(k\right)e^{j2\pi n{\frac {k}{N}}}\right)e^{-j2\pi fn}={\frac {1}{N}}\sum _{k=0}^{N-1}X\left(k\right)\left(\sum _{n=0}^{N-1}e^{-j2\pi fn}e^{j2\pi n{\frac {k}{N}}}\right)$

$\sum _{n=0}^{N-1}e^{-j2\pi fn}e^{j2\pi n{\frac {k}{N}}}=\sum _{n=0}^{N-1}e^{-j2\pi n\left(f-{\frac {k}{N}}\right)}={\frac {1-e^{-j2\pi N\left(f-{\frac {k}{N}}\right)}}{1-e^{-j2\pi \left(f-{\frac {k}{N}}\right)}}}={\frac {e^{-j\pi N\left(f-{\frac {k}{N}}\right)}\left(e^{j\pi N\left(f-{\frac {k}{N}}\right)}-e^{-j\pi N\left(f-{\frac {k}{N}}\right)}\right)}{e^{-j\pi \left(f-{\frac {k}{N}}\right)}\left(e^{j\pi \left(f-{\frac {k}{N}}\right)}-e^{-j\pi \left(f-{\frac {k}{N}}\right)}\right)}}=e^{-j\pi \left(N-1\right)\left(f-{\frac {k}{N}}\right)}{\frac {\sin {\left(\pi Nf-k\pi \right)}}{\sin {\left(\pi f-\pi {\frac {k}{N}}\right)}}}$

In pratica la DTFT viene approssimata a una DFT definita su un nuovo periodo $N_{1}\gg N$, applicata a una sequenza $x_{z}\left(n\right)$ che non è altro che la sequenza di partenza $x\left(n\right)$ a cui sono stati accodati $N_{1}-N$ campioni nulli:

$x_{z}\left(n\right)={\begin{cases}x\left(n\right)&n=0,\ldots ,N-1\\0&n=N,\ldots ,N_{1}-1\end{cases}}$

Aggiungere degli zeri in fondo alla sequenza $x\left(n\right)$ permette di aumentare artificialmente il periodo della sua DFT senza alterare la DTFT a cui si cerca di approssimare:

$X\left(e^{j2\pi f_{k}}\right)=\sum _{n=0}^{N-1}x\left(n\right)e^{-j2\pi f_{k}n}=\sum _{n=0}^{N_{1}-1}x_{z}\left(n\right)e^{-j2\pi {\frac {k}{N_{1}}}n}\quad k=0,\ldots ,N_{1}-1$

e siccome il numero di campioni presi dalla DTFT è maggiore, la risoluzione della DTFT risulta molto più alta.

Esempio - Sequenza porta

$x\left(n\right)={\begin{cases}1&0\leq n\leq N-1\\0&n\geq N\end{cases}}$

Considerando solamente l'intervallo $\left[0,N-1\right]$ si ottiene come DFT $\left|X\left(k\right)\right|$ una funzione $\mathrm {sinc}$ campionata a una frequenza molto bassa:

$X\left(k\right)=\sum _{n=0}^{N-1}x\left(n\right)e^{-j2\pi n{\frac {k}{N}}}=\sum _{n=0}^{N-1}e^{-j2\pi n{\frac {k}{N}}}={\frac {1-e^{-j2\pi k}}{1-e^{-j2\pi {\frac {k}{N}}}}}={\begin{cases}0&k\neq 0\\N&k=0\end{cases}}$

$N_{1}=N$

Basta però aggiungere degli zeri a destra della porta per migliorare la risoluzione della DTFT:

$N_{1}=2N$

---

$N_{1}=4N$

---

$N_{1}=100N$

Aliasing nel tempo

Partendo dalla sequenza $x\left(n\right)$ di durata non finita:

$x\left(n\right)=a^{n}u\left(n\right)\quad 0<a<1$

la sua DTFT campionata nel primo periodo $N$:

$X\left(k\right)=X\left(e^{j2\pi {\frac {k}{N}}}\right)={\frac {1}{1-a^{-j2\pi {\frac {k}{N}}}}}\quad k=0,\ldots ,N-1$

non coincide esattamente con la DFT calcolata sulla sequenza troncata:

$X_{N}\left(k\right)={\frac {1-a^{N}}{1-a^{-j2\pi {\frac {k}{N}}}}}\quad k=0,\ldots ,N-1$

perché vi è un effetto di aliasing nel tempo dovuto al fatto che la sequenza $x\left(n\right)$ di partenza non ha durata finita.

Al crescere dell'ampiezza del periodo $N$, la DFT della sequenza troncata tende sempre di più alla DTFT campionata nel primo periodo.

Proprietà della DFT

Accorgimenti

Le proprietà della DFT sono analoghe a quelle della DTFT, ma ci vogliono alcuni accorgimenti:

• la DFT $X\left(k\right)$ è periodica di periodo $N$ → le operazioni sulla DFT $X\left(k\right)$ devono condurre a una sequenza periodica di periodo $N$;
• la IDFT $x\left(n\right)$ è periodica di periodo $N$ → le operazioni sulla IDFT $x\left(n\right)$ devono condurre a una sequenza periodica di periodo $N$.

Operatore di modulo

L'operatore di modulo restituisce sempre un numero compreso in $\left[0,N-1\right]$:

${\left|k\right|}_{N}=k{\text{ mod }}N$

Se $k$ è un numero negativo, si somma $N$ a $k$ tante volte fino a ottenere un numero compreso in $\left[0,N-1\right]$:

${\left|-13\right|}_{16}=3$

Operatore di ritardo circolare

L'operatore di ritardo circolare restituisce una sequenza ritardata che è ancora compresa in $\left[0,N-1\right]$:

$x\left({\left|n-N_{0}\right|}_{k}\right)$

in quanto i campioni che vanno al di là del primo periodo ricompaiono all'inizio del primo periodo stesso.

Modo operativo

• si genera la sequenza periodicizzata ${\overline {x}}\left(n\right)$;
• si applica il ritardo di $N_{0}$ campioni a ${\overline {x}}\left(n\right)$: ${\overline {x}}\left(n-N_{0}\right)$;
• la sequenza ritardata $x\left(n-N_{0}\right)$ è pari alla sequenza periodicizzata ${\overline {x}}\left(n-N_{0}\right)$ troncata nel primo periodo $N$ ($n\in \left[0,N-1\right]$).

Convoluzione circolare

La convoluzione circolare tra due sequenze $x\left(n\right)$ e $y\left(n\right)$ è definita:

$x\left(n\right)\otimes y\left(n\right)=\sum _{k=0}^{N-1}x\left(k\right)y\left({\left|n-k\right|}_{N}\right)$

dove $N$ è la più grande tra le durate delle singole sequenze.

Proprietà commutativa

$x\left(n\right)\otimes y\left(n\right)=y\left(n\right)\otimes x\left(n\right)$

La convoluzione circolare si può rappresentare con un prodotto matrice per vettore. Ad esempio, se la durata $N$ è pari a 4:

$x\left(n\right)\otimes y\left(n\right)={\begin{bmatrix}y\left(0\right)&y\left(3\right)&y\left(2\right)&y\left(1\right)\\y\left(1\right)&y\left(0\right)&y\left(3\right)&y\left(2\right)\\y\left(2\right)&y\left(1\right)&y\left(0\right)&y\left(3\right)\\y\left(3\right)&y\left(2\right)&y\left(1\right)&y\left(0\right)\end{bmatrix}}{\begin{bmatrix}x\left(0\right)\\x\left(1\right)\\x\left(2\right)\\x\left(3\right)\end{bmatrix}}$

Linearità

$z\left(n\right)=a_{1}x\left(n\right)+a_{2}y\left(n\right)\Rightarrow Z\left(k\right)=a_{1}X\left(k\right)+a_{2}Y\left(k\right)$

Ritardo

$z\left(n\right)=x\left({\left|n-N_{0}\right|}_{N}\right)\Rightarrow Z\left(k\right)=X\left(k\right)e^{-j2\pi {\frac {k}{N}}N_{0}}$

Modulazione

$z\left(n\right)=e^{j2\pi {\frac {k_{0}}{N}}n}x\left(n\right)\Rightarrow Z\left(k\right)=X\left({\left|k-k_{0}\right|}_{N}\right)$

Convoluzione circolare

$z\left(n\right)=x\left(n\right)\otimes y\left(n\right)\Rightarrow Z\left(k\right)=X\left(k\right)\cdot Y\left(k\right)$

Accorgimenti

• la convoluzione circolare di due sequenze di durata $N$ campioni ha una uguale durata di $N$ campioni (a differenza della convoluzione lineare della DTFT, la cui sequenza risultante ha durata $2N-1$ campioni);
• la convoluzione circolare ${\overline {z}}\left(k\right)$ di due IDFT ${\overline {x}}\left(k\right)$ e ${\overline {y}}\left(k\right)$, entrambe di periodo $N$, è periodica di periodo $N$:

  ${\overline {z}}\left(n+N\right)=\sum _{k=0}^{N-1}{\overline {x}}\left(k\right){\overline {y}}\left(n+N-k\right)=\sum _{k=0}^{N-1}{\overline {x}}\left(k\right){\overline {y}}\left(n-k\right)={\overline {z}}\left(n\right)$

Proprietà di simmetria

Se la sequenza $x\left(n\right)$ è reale:

• $X\left(k\right)=X_{R}\left(k\right)+jX_{I}\left(k\right)$
  • se $x\left(n\right)$ è una sequenza pari: $X_{I}\left(k\right)=0$;
  • se $x\left(n\right)$ è una sequenza dispari: $X_{R}\left(k\right)=0$;
• per la DFT vale la simmetria hermitiana intorno a 0:

  $X\left(k\right)=X^{*}\left({\left|-k\right|}_{N}\right)$

• per la DFT vale la simmetria hermitiana intorno a ${\tfrac {N}{2}}$;
• il campione in $k=0$ della DFT è sempre reale:

  $X^{*}\left(0\right)=X\left(0\right)=\sum _{n=0}^{N-1}x\left(n\right)$

Trasformata di Fourier veloce (FFT)

La DFT ha complessità quadratica ($N^{2}$).

La FFT è un algoritmo di complessità linearitmica ($N\log {N}$) che implementa la DFT e la IDFT in maniera più efficiente.

DFT

La DFT può essere rappresentata in termini di matrici:

$\mathbf {X} =\mathbf {Hx} ;\;{\begin{bmatrix}X\left(0\right)\\X\left(1\right)\\X\left(2\right)\\\vdots \\X\left(N-1\right)\end{bmatrix}}={\begin{bmatrix}1&1&1&\cdots &1\\1&H_{N}^{1}&H_{N}^{2}&\cdots &H_{N}^{N-1}\\1&H_{N}^{2}&H_{N}^{4}&\cdots &H_{N}^{2\left(N-1\right)}\\\vdots &\vdots &\vdots &\ddots &\vdots \\1&H_{N}^{N-1}&H_{N}^{2\left(N-1\right)}&\cdots &H_{N}^{\left(N-1\right)\left(N-1\right)}\end{bmatrix}}\cdot {\begin{bmatrix}x\left(0\right)\\x\left(1\right)\\x\left(2\right)\\\vdots \\x\left(N-1\right)\end{bmatrix}}$

dove:

$H_{N}=e^{-j2\pi {\frac {1}{N}}}\Rightarrow H_{N}^{nk}=e^{-j2\pi n{\frac {k}{N}}}$

Algoritmo della decimazione nel tempo

L'algoritmo della decimazione nel tempo è un algoritmo di tipo "divide et impera" che riduce la complessità dell'algoritmo di DFT.

Ipotesi

$N$ è una potenza di 2

Divide

La sequenza si suddivide in due sottosequenze costituite da metà campioni:

• la sottosequenza $x_{0}\left(n\right)$ è formata dai campioni pari:

  $x_{0}\left(n\right)=x\left(2n\right)\quad n=0,\ldots ,{\frac {N}{2}}-1$

• la sottosequenza $x_{1}\left(n\right)$ è formata dai campioni dispari:

  $x_{1}\left(n\right)=x\left(2n+1\right)\quad n=0,\ldots ,{\frac {N}{2}}-1$

Ricombinazione

Ogni DFT $X\left(k\right)$, elemento del vettore $\mathbf {x}$, si ottiene dalla somma delle DFT delle singole sottosequenze $x_{0}\left(n\right)$ e $x_{1}\left(n\right)$:

$X\left(k\right)=X_{0}\left({\left|k\right|}_{\frac {N}{2}}\right)+H_{N}^{k}X_{1}\left({\left|k\right|}_{\frac {N}{2}}\right)=\sum _{n=0}^{{\frac {N}{2}}-1}x_{0}\left(n\right)H_{\frac {N}{2}}^{nk}+H_{N}^{k}\sum _{n=0}^{{\frac {N}{2}}-1}x_{1}\left(n\right)H_{\frac {N}{2}}^{nk}$

Dimostrazione

$X\left(k\right)=\sum _{n=0}^{N-1}x\left(n\right)H_{N}^{nk}=\sum _{n=0}^{{\frac {N}{2}}-1}x_{0}\left(n\right)H_{N}^{2nk}+\sum _{n=0}^{{\frac {N}{2}}-1}x_{1}\left(n\right)H_{N}^{\left(2n+1\right)k}=\sum _{n=0}^{{\frac {N}{2}}-1}x_{0}\left(n\right)H_{N}^{2nk}+H_{N}^{k}\sum _{n=0}^{{\frac {N}{2}}-1}x_{1}\left(n\right)H_{N}^{2nk}=$

Siccome $H_{N}^{2}=H_{\frac {N}{2}}$:

$=\sum _{n=0}^{{\frac {N}{2}}-1}x_{0}\left(n\right)H_{\frac {N}{2}}^{nk}+H_{N}^{k}\sum _{n=0}^{{\frac {N}{2}}-1}x_{1}\left(n\right)H_{\frac {N}{2}}^{nk}$

Complessità

Il numero totale di operazioni (somme e prodotti) complesse è pari a:

$N+{N^{2} \over 2}<N^{2}$

Si può quindi riapplicare ripetutamente il procedimento "divide et impera" sulle sottosequenze. Al passo $k$-esimo si ha una complessità pari a:

$k\cdot N+2^{k}{\left({\frac {N}{2^{k}}}\right)}^{2}\quad k=1,\ldots ,\log _{2}N$

All'ultimo passo ($k=\log _{2}N$), quando si arriva a dover valutare le DFT su 1 punto, la complessità totale è linearitmica:

$\log _{2}N\cdot N+N\approx \log _{2}N\cdot N$

Riduzione di complessità

La complessità può essere ulteriormente ridotta sfruttando le proprietà di simmetria della matrice $\mathbf {H}$.

Esempio con $N=4$ campioni

${\begin{cases}X\left(0\right)=x\left(0\right)+x\left(2\right)+x\left(1\right)+x\left(3\right)\\X\left(1\right)=x\left(0\right)+H_{2}x\left(2\right)+H_{4}\left[x\left(1\right)+H_{2}x\left(3\right)\right]\\X\left(2\right)=x\left(0\right)+x\left(2\right)+H_{4}^{2}\left[x\left(1\right)+x\left(3\right)\right]\\X\left(3\right)=x\left(0\right)+H_{2}x\left(2\right)+H_{4}^{3}\left[x\left(1\right)+H_{2}x\left(3\right)\right]\end{cases}}$

Dimostrazione

$X\left(k\right)=X_{0}\left({\left|k\right|}_{\frac {N}{2}}\right)+H_{N}^{k}X_{1}\left({\left|k\right|}_{\frac {N}{2}}\right)=X_{0}\left({\left|k\right|}_{2}\right)+H_{4}^{k}X_{1}\left({\left|k\right|}_{2}\right)\quad k=0,\ldots ,3$

${\begin{cases}X_{0}\left(k\right)=x\left(0\right)+x\left(2\right)H_{2}^{k}\\X_{1}\left(k\right)=x\left(1\right)+x\left(3\right)H_{2}^{k}\end{cases}}\Rightarrow {\begin{cases}X_{0}\left(0\right)=x\left(0\right)+x\left(2\right)\\X_{0}\left(1\right)=x\left(0\right)+x\left(2\right)H_{2}\\X_{1}\left(0\right)=x\left(1\right)+x\left(3\right)\\X_{1}\left(1\right)=x\left(1\right)+x\left(3\right)H_{2}\end{cases}}\Rightarrow {\begin{cases}X\left(0\right)=X_{0}\left(0\right)+X_{1}\left(0\right)=x\left(0\right)+x\left(2\right)+x\left(1\right)+x\left(3\right)\\X\left(1\right)=X_{0}\left(1\right)+H_{4}X_{1}\left(1\right)=x\left(0\right)+H_{2}x\left(2\right)+H_{4}\left[x\left(1\right)+H_{2}x\left(3\right)\right]\\X\left(2\right)=X_{0}\left(0\right)+H_{4}^{2}X_{1}\left(0\right)=x\left(0\right)+x\left(2\right)+H_{4}^{2}\left[x\left(1\right)+x\left(3\right)\right]\\X\left(3\right)=X_{0}\left(1\right)+H_{4}^{3}X_{1}\left(1\right)=x\left(0\right)+H_{2}x\left(2\right)+H_{4}^{3}\left[x\left(1\right)+H_{2}x\left(3\right)\right]\end{cases}}$

Schema circuitale a farfalla

IDFT

Tutte le considerazioni sulla DFT valgono anche per la IDFT: basta applicare l'algoritmo sul complesso coniugato della DFT e poi dividere per $N$:

$x\left(n\right)={\frac {1}{N}}\sum _{k=0}^{N-1}X\left(k\right)e^{j2\pi n{\frac {k}{N}}}={\frac {1}{N}}{\left[\sum _{k=0}^{N-1}X^{*}\left(k\right)e^{-j2\pi n{\frac {k}{N}}}\right]}^{*}\Rightarrow x\left(n\right)={\text{IDFT}}\left[X\left(k\right)\right]={\frac {1}{N}}{\text{DFT}}\left[X^{*}\left(k\right)\right],\quad n=0,1,2,\ldots ,N-1$