DFT, matrice di Fourier e radici n-esime dell'unità

La nota precedente (Polinomi - rappresentazione per coefficienti e per valori e convoluzioneun polinomio di grado limitato da n si rappresenta con n coefficienti o con n valori su n punti distinti (teorema di interpolazione, matrice di Vandermonde). Per coefficienti la somma costa n e il prodotto, la convoluzione lineare, n^2 prodotti e (n-1)^2 somme; per valori il prodotto costa 2n-1 prodotti ma serve una base estesa di almeno 2n-1 punti. Le conversioni generali costano troppo: valutazione n volte con Horner 2n(n-1), interpolazione con Gauss-Jordan n^3 o con Lagrange n^2. Serve una base speciale: la DFT.Polinomi - rappresentazione per coefficienti e per valori e convoluzione →) ha mostrato che per moltiplicare due polinomi in $o(n^2)$ servirebbe una base di valutazione per cui le conversioni costino poco. Qui si sceglie la base: le radici $n$-esime dell'unità. Si definisce la DFT, si dimostrano le proprietà delle radici che la rendono veloce e il teorema che la collega alla convoluzione lineare. Il calcolo veloce (FFT) è nella nota Fast Fourier Transform - algoritmo ricorsivo e trasformata inversala FFT calcola DFT_n (n potenza di due) con divide et impera: si separano i coefficienti di indice pari e dispari, p(x) = p0(x^2) + x p1(x^2), e per il lemma di dimezzamento bastano due DFT di ordine n/2; y_i = y0_i + w^i y1_i e y_{i+n/2} = y0_i - w^i y1_i. La ricorrenza T(n) = 2T(n/2) + 2n con T(1) = 0 ha soluzione esatta 2 n log2 n. L'inversa ha matrice F^{-1}_{ij} = w^{-ij}/n: si ottiene con la FFT sulle potenze negative (poi si divide per n fuori dalla ricorsione) oppure con DFT_n(y) rovesciato. LIN_CONV moltiplica due polinomi con n coefficienti in 12 n log2(2n) + 4n operazioni.Fast Fourier Transform - algoritmo ricorsivo e trasformata inversa →.

Il punto di vista dei corsi di segnali è diverso: là la DFT è la serie di Fourier di un segnale periodico (Trasformata di Fourier discreta (TFD)La TFD è la serie di Fourier dei segnali discreti periodici di periodo $N$: servono solo $N$ armoniche, $X(k)=\frac1N\sum_{n=0}^{N-1}x(n)e^{-j2\pi kn/N}$ e $x(n)=\sum_{k=0}^{N-1}X(k)e^{j2\pi kn/N}$. È un prodotto matrice-vettore in $\mathbb{C}^N$ (la matrice di Fourier); la convoluzione circolare diventa prodotto, $N X(k)Y(k)$. La FFT la calcola in $O(N\log N)$.Trasformata di Fourier discreta (TFD) →, Trasformata di Fourier discreta (DFT) e FFTUn segnale discreto periodico di periodo $NT$ ha una trasformata discreta e periodica, la DFT: $S(kF)=\sum_{n=0}^{N-1}T,s(nT)e^{-i2\pi kn/N}$, con $F=1/(NT)$, e $s(nT)=\sum_{k=0}^{N-1}F,S(kF)e^{i2\pi kn/N}$. Dipende da soli $N$ numeri, si calcola senza approssimazioni e con la FFT costa $N\log_2N$ invece di $N^2$. I coefficienti di Fourier del segnale periodico sono $S_k=F,S(kF)$. La DFT dà anche campioni della trasformata di un segnale discreto o continuo a durata limitata (con zero-padding), ma la scalatura $T$ e l'asse delle frequenze vanno gestiti con cura.Trasformata di Fourier discreta (DFT) e FFT →, Trasformata di Fourier discreta (DFT) e convoluzione circolareLa DFT di N campioni x[0..N-1] è X[k] = sum x[n] e^{-j2pi kn/N}, k = 0..N-1 (IDFT: x[n] = (1/N) sum X[k] e^{j2pi kn/N}); vede il segnale come periodico di periodo N e dà N campioni equispaziati della DTFT, cioè X(z) valutata sulle radici N-esime dell'unità. Se x ha durata M <= N non c'è aliasing temporale e lo zero-padding (N > M) infittisce solo i campioni della stessa DTFT. Proprietà: traslazione circolare, simmetria coniugata X[k] = X*[N-k], prodotto = convoluzione circolare. La convoluzione circolare coincide con la lineare solo se N >= Dx + M - 1; altrimenti la coda si ripiega sull'inizio. La FFT calcola la DFT con (N/2) log2 N moltiplicazioni invece di N^2.Trasformata di Fourier discreta (DFT) e convoluzione circolare →) e ha il segno meno nell'esponente. Qui è la valutazione di un polinomio in punti speciali, e la convenzione è quella delle dispense del corso, con il segno più: $\omega_n=e^{+2\pi i/n}$. Le due convenzioni differiscono per un'inversione degli indici (vedi più sotto); numpy.fft.fft usa il segno meno.

Numeri complessi in forma polare

Un numero complesso $z=a+ib\neq0$ ha modulo $\rho=\sqrt{a^2+b^2}$ e argomento $\theta$ (angolo con il semiasse reale positivo): $z=\rho(\cos\theta+i\sin\theta)=\rho\,e^{i\theta}$ (formula di Eulero, richiami in Numeri complessi, formula di Eulero ed esponenziali complessiUn numero complesso si scrive in forma cartesiana a+jb o polare |x|e^{jφ}; il prodotto moltiplica i moduli e somma le fasi. La formula di Eulero e^{jα}=cos α + j sin α lega esponenziali e sinusoidi e permette di trattare tutti i segnali del corso come somme di esponenziali complessi e^{(σ+jω)t}.Numeri complessi, formula di Eulero ed esponenziali complessi →). La forma esponenziale rende semplici potenze e radici: $z^n=\rho^n e^{in\theta}$.

Radici $n$-esime. Se $z_k^n=z=\rho e^{i\theta}$ allora $\rho_k^n=\rho$, cioè $\rho_k=\rho^{1/n}$, e $n\theta_k$ differisce da $\theta$ per un multiplo di $2\pi$: $n\theta_k=\theta+2\pi k$. Le $n$ radici distinte sono $$z_k=\rho^{1/n}\,e^{\,i\left(\frac{\theta}{n}+\frac{2\pi k}{n}\right)},\qquad 0\le k\le n-1.$$ Per $z=1$ ($\rho=1$, $\theta=0$) si ottengono le radici $n$-esime dell'unità $e^{2\pi ik/n}$.

Definizione (radice principale). La radice $n$-esima principale dell'unità è $\omega_n=e^{2\pi i/n}$. Le radici $n$-esime dell'unità sono le sue prime $n$ potenze, $\omega_n^0,\omega_n^1,\dots,\omega_n^{n-1}$, e formano il vettore $\Omega_n=(\omega_n^0,\dots,\omega_n^{n-1})$. Poiché $\omega_n^n=1$, vale $\omega_n^k=\omega_n^{\,k\bmod n}$ per ogni $k\in\mathbb Z$ (aritmetica modulare sugli esponenti; $\omega_n^{-1}=\omega_n^{n-1}$).

Esempio ($n=8$). $\omega_8=e^{i\pi/4}=\frac{1+i}{\sqrt2}$ e le potenze sono i vertici di un ottagono regolare inscritto nella circonferenza di raggio $1$:

Per $n=4$: $\omega_4=i$ e $\Omega_4=(1,i,-1,-i)$. Per $n=2$: $\omega_2=-1$ e $\Omega_2=(1,-1)$.

La DFT e la matrice di Fourier

Definizione (DFT). La trasformata discreta di Fourier di ordine $n$ è la funzione $\mathrm{DFT}_n:\mathbb C^n\to\mathbb C^n$ che a $\mathbf a$ (coefficienti di $p$) associa le valutazioni di $p$ sulle radici $n$-esime dell'unità: $$y_i=p(\omega_n^i)=\sum_{j=0}^{n-1}a_j\,\omega_n^{ij},\qquad 0\le i\le n-1.$$ La coppia $(\Omega_n,\mathbf y)$ è una rappresentazione per valori di $p$ (le radici sono distinte, quindi vale il teorema di interpolazione).

Come per ogni valutazione, $\mathrm{DFT}_n(\mathbf a)=F_n\mathbf a$, dove $F_n=V(\Omega_n)$ è la matrice di Fourier, la Vandermonde delle radici: $$[F_n]_{ij}=\omega_n^{\,ij},\qquad 0\le i,j\le n-1.$$ La prima riga e la prima colonna contengono solo $1$ (esponente $0$). $F_n$ è simmetrica ($ij=ji$) e invertibile (Vandermonde su punti distinti), quindi $\mathrm{DFT}_n$ è lineare e invertibile: l'inversa è l'interpolazione da $\Omega_n$.

Esempio ($n=4$). Con $\omega_4=i$ e $\omega_4^{ij}=i^{\,ij\bmod4}$: $$F_4=\begin{pmatrix}1&1&1&1\\1&i&-1&-i\\1&-1&1&-1\\1&-i&-1&i\end{pmatrix}.$$ Per $p(x)=3+4x+5x^2+6x^3\equiv(3,4,5,6)$ si ha $y_0=3+4+5+6=18$, $y_1=3+4i-5-6i=-2-2i$, $y_2=3-4+5-6=-2$, $y_3=3-4i-5+6i=-2+2i$. Quindi $\mathrm{DFT}_4(3,4,5,6)=(18,\,-2-2i,\,-2,\,-2+2i)$ (controllato con numpy, dopo aver cambiato il segno dell'esponente).

Costo con la valutazione generale. $n$ applicazioni di Horner: $2n(n-1)=\Theta(n^2)$ operazioni.

Legame con la convenzione dei segnali. Con $\bar\omega_n=\omega_n^{-1}=e^{-2\pi i/n}$ la "DFT dei segnali" è $X_k=\sum_jx_j\bar\omega_n^{\,jk}$. Poiché $\bar\omega_n^{\,k}=\omega_n^{-k}=\omega_n^{(n-k)\bmod n}$, i due vettori differiscono per la permutazione $k\mapsto(n-k)\bmod n$, che lascia fermo l'indice $0$ e rovescia gli altri. La stessa permutazione compare nella trasformata inversa (nota sulla FFT).

Proprietà delle radici dell'unità

Gruppo. $\Omega_n$ è un gruppo moltiplicativo: la chiusura è $\omega_n^i\omega_n^j=\omega_n^{(i+j)\bmod n}$; l'elemento neutro è $1=\omega_n^0$; l'inverso di $\omega_n^i$ è $\omega_n^{-i}=\omega_n^{(-i)\bmod n}$, che è $\omega_n^{n-i}$ per $1\le i\le n-1$ e $1$ per $i=0$.

Lemma di cancellazione. Per $n,d$ interi positivi e $k\in\mathbb Z$: $\omega_{dn}^{dk}=\omega_n^k$.

Dimostrazione. $\omega_{dn}^{dk}=\bigl(e^{2\pi i/(dn)}\bigr)^{dk}=e^{2\pi ik/n}=\omega_n^k$: il fattore $d$ si semplifica nell'esponente. $\square$

Corollario. Per $n>0$ pari, $\omega_n^{n/2}=\omega_2^1=-1$. Infatti $\omega_n^{n/2}=\omega_{2\cdot(n/2)}^{1\cdot(n/2)}=\omega_2$ per cancellazione, e $\omega_2=e^{i\pi}=-1$. Geometricamente: dopo metà delle potenze il vettore ha percorso mezzo giro. Quindi $\omega_n^{k+n/2}=-\omega_n^k$.

Lemma di dimezzamento. Per $n>0$ pari, i quadrati $(\omega_n^i)^2$, $0\le i\le n-1$, coincidono a coppie: per $0\le i\le n/2-1$ vale $$(\omega_n^i)^2=(\omega_n^{\,i+n/2})^2=\omega_{n/2}^{\,i}.$$ Elevando al quadrato le $n$ radici $n$-esime si ottengono quindi le $n/2$ radici $n/2$-esime, ciascuna due volte.

Dimostrazione. $(\omega_n^i)^2=\omega_n^{2i}=\omega_{2\cdot n/2}^{2i}=\omega_{n/2}^{i}$ per cancellazione. Inoltre $(\omega_n^{i+n/2})^2=\omega_n^{2i}\,\omega_n^{n}=\omega_{n/2}^i\cdot1=\omega_{n/2}^i$. $\square$

Esempio ($n=8$). Quadrando $\omega_8^0,\dots,\omega_8^7$ si ottiene $1,\ i,\ -1,\ -i,\ 1,\ i,\ -1,\ -i$: le quattro radici quarte dell'unità, ciascuna due volte.

Lemma di somma. Per $n>1$ e $k\ge0$ interi: $$\sum_{j=0}^{n-1}\bigl(\omega_n^k\bigr)^j=\begin{cases}n&\text{se }k\text{ è multiplo di }n,\\0&\text{altrimenti.}\end{cases}$$

Dimostrazione. Se $n\mid k$ allora $\omega_n^k=\omega_n^0=1$ e la somma vale $\sum_j1=n$. Altrimenti $\omega_n^k\ne1$ e si usa la somma geometrica (Serie notevoli - geometrica, telescopica, armonicaLe serie di cui si conosce il carattere e da usare come termine di paragone: geometrica (converge a 1/(1-q) se |q|<1), telescopiche (somma b_1 - lim b_n, come Mengoli), armonica generalizzata (1/n^alpha converge se e solo se alpha>1).Serie notevoli - geometrica, telescopica, armonica →, valida anche in $\mathbb C$): $$\sum_{j=0}^{n-1}(\omega_n^k)^j=\frac{(\omega_n^k)^n-1}{\omega_n^k-1}=\frac{(\omega_n^n)^k-1}{\omega_n^k-1}=\frac{1-1}{\omega_n^k-1}=0.\ \square$$

Esempio. $n=4$, $k=1$: $1+i+i^2+i^3=1+i-1-i=0$. $k=4$: $1+1+1+1=4$. $k=2$: $1+(-1)+1+(-1)=0$.

Osservazione. Tutti i risultati valgono identicamente con $\bar\omega_n=\omega_n^{-1}$ al posto di $\omega_n$ (le radici sono le stesse, in un altro ordine): serve per la trasformata inversa.

Esempio d'esame (Re e Im di $\omega_{16}^{120}$). Per periodicità $\omega_{16}^{120}=\omega_{16}^{120\bmod16}=\omega_{16}^8=\omega_2^1=-1$ (con $120=7\cdot16+8$ e cancellazione). In modo più rapido, per cancellazione ripetuta: $\omega_{16}^{120}=\omega_8^{60}=\omega_4^{30}=\omega_2^{15}=\omega_2^1=-1$ (perché $\omega_2^{15}=(-1)^{15}$). Quindi $\mathrm{Re}=-1$ e $\mathrm{Im}=0$.

Il teorema di convoluzione lineare

Teorema (convoluzione lineare). Siano $\mathbf a,\mathbf b\in\mathbb C^n$ e $\mathbf c=\mathbf a\star\mathbf b\in\mathbb C^{2n-1}$. Allora $$(\mathbf c\,|\,0)=\mathrm{DFT}^{-1}_{2n}\Bigl(\mathrm{DFT}_{2n}\bigl((\mathbf a|\mathbf 0_n)\bigr)\odot\mathrm{DFT}_{2n}\bigl((\mathbf b|\mathbf 0_n)\bigr)\Bigr),$$ dove $|$ indica la concatenazione (padding con zeri), $\odot$ il prodotto componente per componente.

Dimostrazione. Siano $p\equiv\mathbf a$, $q\equiv\mathbf b$ e $c=pq\equiv\mathbf a\star\mathbf b$. Il padding non cambia i polinomi, ma li porta a limite di grado $2n$: $(\mathbf a|\mathbf 0_n)\equiv p$ e $(\mathbf b|\mathbf 0_n)\equiv q$. Le due DFT di ordine $2n$ sono le valutazioni $p(\omega_{2n}^i)$ e $q(\omega_{2n}^i)$, $0\le i\le 2n-1$. Il prodotto componente per componente è $p(\omega_{2n}^i)q(\omega_{2n}^i)=c(\omega_{2n}^i)$: sono le valutazioni di $c$ su $\Omega_{2n}$. Il polinomio $c$ ha grado al più $2n-2<2n$, quindi $2n$ punti bastano a determinarlo e la sua rappresentazione per coefficienti con limite di grado $2n$ è $(\mathbf c|0)$. Per il teorema di interpolazione, l'antitrasformata restituisce esattamente questo vettore. $\square$

Si usa l'ordine $2n$ e non $2n-1$ perché l'algoritmo veloce (FFT) lavora su potenze di due: raddoppiare $n$ mantiene la forma giusta e basta un solo zero di padding in più.

Esempio ($n=2$). $\mathbf a=(1,2)$, $\mathbf b=(3,4)$, ordine $4$. Con $F_4$:

Il polinomio con coefficienti $(3,10,8,0)$, cioè $3+10x+8x^2$, vale $21$ in $1$, $3+10i-8=-5+10i$ in $i$, $3-10+8=1$ in $-1$ e $-5-10i$ in $-i$: sono i valori del prodotto. Quindi $\mathbf c=(3,10,8)$, che è anche $(1+2x)(3+4x)=3+10x+8x^2$.

Il costo di questa strada dipende da come si calcolano le DFT: con le valutazioni generali $\Theta(n^2)$ (nessun vantaggio), con la FFT $\Theta(n\log n)$ (vedi Fast Fourier Transform - algoritmo ricorsivo e trasformata inversala FFT calcola DFT_n (n potenza di due) con divide et impera: si separano i coefficienti di indice pari e dispari, p(x) = p0(x^2) + x p1(x^2), e per il lemma di dimezzamento bastano due DFT di ordine n/2; y_i = y0_i + w^i y1_i e y_{i+n/2} = y0_i - w^i y1_i. La ricorrenza T(n) = 2T(n/2) + 2n con T(1) = 0 ha soluzione esatta 2 n log2 n. L'inversa ha matrice F^{-1}_{ij} = w^{-ij}/n: si ottiene con la FFT sulle potenze negative (poi si divide per n fuori dalla ricorsione) oppure con DFT_n(y) rovesciato. LIN_CONV moltiplica due polinomi con n coefficienti in 12 n log2(2n) + 4n operazioni.Fast Fourier Transform - algoritmo ricorsivo e trasformata inversa →). Lo stesso legame spiega perché le matrici circolanti si diagonalizzano con $F_n$ (Matrici circolanti e supercircolanti - algoritmi divide and conquerle matrici supercircolanti (quadranti alternati ricorsivi) hanno solo n valori distinti: somma in n operazioni, prodotto in 4n log2 n + n con il trucco (A1+A2)(B1+B2) e (A1-A2)(B1-B2) invece delle 4 chiamate (Theta(n^2)); gli autovalori sono la trasformata di Hadamard della prima riga (FHT, n log2 n) e il prodotto passa dalla base diagonalizzante; le matrici circolanti sono polinomi nella matrice di rotazione R, il loro prodotto è la convoluzione ciclica delle prime righe e si diagonalizzano con la matrice di Fourier: prodotto in O(n log n) con la FFT.Matrici circolanti e supercircolanti - algoritmi divide and conquer →).

Domande d'esame

• Determinare parte reale e immaginaria di $\omega_{16}^{120}$. Traccia: $120\bmod16=8$, $\omega_{16}^8=\omega_2=-1$, quindi $(-1,0)$.
• Quanto vale $\sum_{j=0}^{7}(\omega_8^{3})^j$? E $\sum_{j=0}^{7}(\omega_8^{16})^j$? Traccia: lemma di somma; $3$ non è multiplo di $8$ quindi $0$; $16$ è multiplo di $8$ quindi $8$.
• Perché nel teorema di convoluzione lineare si usa la DFT di ordine $2n$ e non $n$? Traccia: il prodotto ha $2n-1$ coefficienti; con $n$ punti la rappresentazione per valori non lo determina (serve una base di almeno $2n-1$ punti); si arrotonda a $2n$ per avere una potenza di due.

Esercizi svolti: Esercizio - domande di prima parte su polinomi, DFT e FFT (simulazioni d'esame); Esercizio - polinomio con tutte le radici dell'unità tranne 1 e -1, valutazione e DFT dei coefficienti (compito del 26 giugno 2025, es. 1).