Fast Fourier Transform - algoritmo ricorsivo e trasformata inversa

La FFT è un algoritmo divide and conquer (vedi Divide and conquer - proprietà di sottostruttura e albero delle chiamateUn algoritmo divide and conquer nasce da una proprietà di sottostruttura (funzioni Divide e Conquer, base n0) con sottoistanze di taglia strettamente minore; l'albero delle chiamate ne descrive l'esecuzione, e il D&C non ha memoria (ripete le sottoistanze uguali). Esempio completo: MAX di un array, con correttezza per induzione, ricorrenza e conteggio esatto n-1 tramite il teorema sugli alberi binari pieni.Divide and conquer - proprietà di sottostruttura e albero delle chiamate →) per calcolare la DFT (definita in DFT, matrice di Fourier e radici n-esime dell'unitàla DFT di ordine n valuta un polinomio di grado limitato da n sulle n radici n-esime dell'unità, y_i = p(ω^i) con ω = e^{2πi/n}; è il prodotto per la matrice di Fourier F_n, una Vandermonde invertibile. Le radici formano un gruppo moltiplicativo e valgono tre lemmi: cancellazione (ω_{dn}^{dk} = ω_n^k, ω_n^{n/2} = -1), dimezzamento (i quadrati delle n radici sono le n/2 radici n/2-esime, ciascuna due volte) e somma (la somma delle potenze di ω^k vale n se k è multiplo di n, 0 altrimenti). Teorema di convoluzione lineare: (c|0) = DFT^{-1}{2n}(DFT{2n}(a|0) ⊙ DFT_{2n}(b|0)).DFT, matrice di Fourier e radici n-esime dell'unità →) in tempo $\Theta(n\log n)$ invece di $\Theta(n^2)$. Qui il punto di vista è quello dell'algoritmo: sottostruttura, correttezza per induzione, ricorrenza con soluzione esatta. Gli stessi passi, letti come fattorizzazione della matrice di Fourier e come analisi spettrale, sono in 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 → e FFT e zero-padding - TF, TFtd e asse delle pulsazioniLa FFT dei campioni di un segnale a durata finita, moltiplicata per il passo $T_c$, approssima la trasformata di Fourier: $X(\omega_k)\approx T_c,\mathtt{fft}(x,M)[k]$ con $\omega_k=\frac{2\pi k}{MT_c}$ (e fattore di fase $e^{-j\omega t_0}$ se l'asse parte da $t_0$). Lo zero-padding ($M>N$) infittisce i punti della stessa TFtd senza aggiungere informazione; la risoluzione dipende dalla durata osservata. Per un segnale reale $|X|$ è simmetrico: il picco in $k$ ha un gemello in $M-k$.FFT e zero-padding - TF, TFtd e asse delle pulsazioni →.

D'ora in avanti $n$ è una potenza di due. La restrizione non pesa per la convoluzione: si può sempre fare il padding fino alla potenza di due successiva.

La proprietà di sottostruttura

Caso base. $n=1$: il polinomio è costante, $p(x)=a_0$, e la sua valutazione in $\omega_1^0=1$ è $a_0$. Quindi $\mathrm{DFT}_1(\mathbf a)=\mathbf a$ (si restituisce comunque un vettore, di una componente).

Divisione. Per $n>1$ si dividono i coefficienti in base alla parità dell'indice: $$\mathbf a^{[0]}=(a_0,a_2,\dots,a_{n-2}),\qquad\mathbf a^{[1]}=(a_1,a_3,\dots,a_{n-1}),$$ che sono i coefficienti dei polinomi di grado limitato da $n/2$ $$p^{[0]}(x)=\sum_{j=0}^{n/2-1}a_{2j}x^j,\qquad p^{[1]}(x)=\sum_{j=0}^{n/2-1}a_{2j+1}x^j.$$ (Non è una divisione in due metà contigue, come nella ricerca del massimo o nei blocchi delle matrici.)

Identità polinomiale. Per ogni $x\in\mathbb C$: $\;p(x)=p^{[0]}(x^2)+x\,p^{[1]}(x^2).$

Dimostrazione. $p^{[0]}(x^2)=\sum_ja_{2j}x^{2j}$ è la somma dei monomi di grado pari di $p$; $x\,p^{[1]}(x^2)=\sum_ja_{2j+1}x^{2j}\cdot x=\sum_ja_{2j+1}x^{2j+1}$ è la somma dei monomi di grado dispari. $\square$

Esempio. $p(x)=3+4x+5x^2+6x^3\equiv(3,4,5,6)$ dà $p^{[0]}=3+5x\equiv(3,5)$ e $p^{[1]}=4+6x\equiv(4,6)$, e infatti $(3+5x^2)+x(4+6x^2)=3+4x+5x^2+6x^3$.

Teorema (sottostruttura). Siano $\mathbf y=\mathrm{DFT}_n(\mathbf a)$, $\mathbf y^{[0]}=\mathrm{DFT}_{n/2}(\mathbf a^{[0]})$, $\mathbf y^{[1]}=\mathrm{DFT}_{n/2}(\mathbf a^{[1]})$. Per $0\le i\le n/2-1$: $$y_i=y^{[0]}_i+\omega_n^i\,y^{[1]}_i,\qquad y_{i+n/2}=y^{[0]}_i-\omega_n^i\,y^{[1]}_i.$$

Dimostrazione. Per le prime $n/2$ componenti: $$\begin{aligned}y_i&=p(\omega_n^i)=p^{[0]}\bigl((\omega_n^i)^2\bigr)+\omega_n^i\,p^{[1]}\bigl((\omega_n^i)^2\bigr)&&\text{(identità polinomiale)}\\&=p^{[0]}(\omega_{n/2}^i)+\omega_n^i\,p^{[1]}(\omega_{n/2}^i)&&\text{(dimezzamento)}\\&=y^{[0]}_i+\omega_n^i\,y^{[1]}_i.\end{aligned}$$ Per le ultime $n/2$, con $0\le i\le n/2-1$: $$\begin{aligned}y_{i+n/2}&=p(\omega_n^{i+n/2})=p^{[0]}\bigl(\omega_n^{2i+n}\bigr)+\omega_n^{i+n/2}\,p^{[1]}\bigl(\omega_n^{2i+n}\bigr)&&\text{(identità polinomiale)}\\&=p^{[0]}(\omega_{n/2}^i)-\omega_n^i\,p^{[1]}(\omega_{n/2}^i)&&\text{(}\omega_n^{n}=1,\ \omega_n^{n/2}=-1\text{)}\\&=y^{[0]}_i-\omega_n^i\,y^{[1]}_i.\ \square\end{aligned}$$

Ogni coppia di valutazioni sullo stesso punto $\omega_{n/2}^i$ produce, con le due combinazioni $y^{[0]}_i\pm\omega_n^iy^{[1]}_i$, due valori distinti di $p$: i termini comuni si calcolano una volta sola. Questa coppia di formule si chiama butterfly (farfalla).

L'algoritmo

FFT(a)
n <- a.len
if n = 1 then return a                      { base }
a0 <- (a[0], a[2], ..., a[n-2])             { divide: indici pari }
a1 <- (a[1], a[3], ..., a[n-1])             { divide: indici dispari }
y0 <- FFT(a0);  y1 <- FFT(a1)               { recurse }
ON <- e^(i 2 pi / n)                        { conquer: ON = w_n }
ONI <- 1                                    { all'inizio dell'iterazione i, ONI = w_n^i }
for i <- 0 to n/2 - 1
    t <- ONI * y1[i]
    y[i] <- y0[i] + t
    y[i + n/2] <- y0[i] - t
    ONI <- ONI * ON
return y

Correttezza. Induzione sull'esponente $k$ con $n=2^k$. Per $k=0$ è il caso base. Per $k\ge1$, la taglia delle sottoistanze è $2^{k-1}$ e per ipotesi induttiva y0 e y1 sono $\mathrm{DFT}_{n/2}(\mathbf a^{[0]})$ e $\mathrm{DFT}_{n/2}(\mathbf a^{[1]})$. Il ciclo calcola, con ONI $=\omega_n^i$ (invariante: alla fine dell'iterazione $i$ vale $\omega_n^{i+1}$), le due formule del teorema di sottostruttura. $\square$

Costo esatto. Ogni iterazione del ciclo fa $4$ operazioni tra complessi: due prodotti (t e ONI) e due somme/differenze. Il ciclo ha $n/2$ iterazioni: lavoro di ricombinazione $4\cdot\frac n2=2n$. Quindi $$T(1)=0,\qquad T(n)=2\,T\!\left(\tfrac n2\right)+2n\quad(n\ge2).$$ Soluzione esatta: $T(n)=2n\log_2n$. Dimostrazione per induzione. $T(1)=0=2\cdot1\cdot\log_21$. Se $T(n/2)=2\cdot\frac n2\log_2\frac n2=n(\log_2n-1)$, allora $T(n)=2n(\log_2n-1)+2n=2n\log_2n$. $\square$ (Verificato contando le operazioni di un'implementazione per $n=2^0,\dots,2^{10}$.)

Lo stesso risultato segue dal Master Theorem - enunciato, casi e applicazioniMaster Theorem per T(n) = aT(n/b) + w(n): si confronta w(n) con la funzione soglia n^(log_b a). Caso 1 (w minore di una potenza polinomiale): Theta(n^(log_b a)); caso 2 (w della stessa taglia): Theta(n^(log_b a) log n); caso 3 (w maggiore e regolare): Theta(w(n)). Prova completa dai livelli dell'albero, versione estesa (parti intere, n0 qualsiasi), esempi, e l'uso scorretto quando nessun caso si applica (w = n log n).Master Theorem - enunciato, casi e applicazioni →: $s=2$ sottoistanze di taglia $n/2$, soglia $n^{\log_22}=n$, lavoro $w(n)=2n$ dello stesso ordine: caso 2, $T(n)=\Theta(n\log n)$. L'albero delle chiamate ha $\log_2n+1$ livelli, $n$ foglie, e a ogni livello il lavoro totale è $2n$ (Albero della ricorrenza e formula generale per le ricorrenzeL'albero della ricorrenza ha la forma dell'albero delle chiamate dell'istanza peggiore e T(n) è la somma dei lavori dei nodi. Per la famiglia T(n) = s(n)T(f(n)) + w(n) con base T0 per n <= n0, contando i nodi per livello con le iterate di f e il livello f*(n,n0) si ottiene la formula chiusa T(n) = somma_l [prodotto_j s(f^(j)(n))] w(f^(l)(n)) + [prodotto s] T0, verificata per induzione e applicata a esempi (merge sort, Hanoi, ricorrenze con radice).Albero della ricorrenza e formula generale per le ricorrenze →; richiami sull'ordine di grandezza in Notazione asintoticaDefinizioni di O, Omega, Theta e o piccolo con le costanti c ed n0; esempi con costanti esplicite; proprietà (polinomi, esponenziali, logaritmi, somme, implicazioni tra notazioni); sommatorie notevoli; terminologia (logaritmica, lineare, polinomiale, esponenziale).Notazione asintotica →).

Osservazione. Se le potenze $\omega_n^i$ sono precalcolate in una tabella, ogni iterazione fa un prodotto e due somme: $3\cdot\frac n2$ operazioni per livello, cioè $\frac32n\log_2n$ in tutto. Il costo $2n\log_2n$ del corso conta anche l'aggiornamento di ONI.

Esempio $n=4$

$\mathbf a=(3,4,5,6)$, $\omega_4=i$. Le sottoistanze sono $(3,5)$ e $(4,6)$; con $\omega_2=-1$, $\mathrm{FFT}(u,v)=(u+v,\,u-v)$:

• $\mathbf y^{[0]}=\mathrm{FFT}(3,5)=(8,-2)$ e $\mathbf y^{[1]}=\mathrm{FFT}(4,6)=(10,-2)$.
• $i=0$ ($\omega_4^0=1$): $y_0=8+10=18$, $y_2=8-10=-2$.
• $i=1$ ($\omega_4^1=i$): $t=i\cdot(-2)=-2i$, $y_1=-2-2i$, $y_3=-2+2i$.

Risultato $(18,-2-2i,-2,-2+2i)$, uguale alla DFT calcolata dalla definizione.

Esempio $n=8$

$\mathbf a=(1,2,3,4,5,6,7,8)$. Le foglie dell'albero, lette da sinistra, sono $a_0,a_4,a_2,a_6,a_1,a_5,a_3,a_7$ (indici in ordine di bit invertiti).

All'ultimo livello $\omega_8=\frac{1+i}{\sqrt2}$ e $\mathbf y^{[0]}=(16,-4-4i,-4,-4+4i)$, $\mathbf y^{[1]}=(20,-4-4i,-4,-4+4i)$:

Quindi $\mathrm{DFT}_8(1,\dots,8)=\bigl(36,\ -4-9{,}657i,\ -4-4i,\ -4-1{,}657i,\ -4,\ -4+1{,}657i,\ -4+4i,\ -4+9{,}657i\bigr)$. Si noti la simmetria coniugata ($y_{n-i}=\overline{y_i}$) che vale perché $\mathbf a$ è reale.

La trasformata inversa

L'inversa $\mathrm{DFT}^{-1}_n(\mathbf y)=F_n^{-1}\mathbf y$ è l'interpolazione da $\Omega_n$ ai coefficienti.

Teorema (matrice inversa di Fourier). $\bigl[F_n^{-1}\bigr]_{ij}=\dfrac1n\,\omega_n^{-ij}$ per $0\le i,j\le n-1$.

Dimostrazione. Sia $G_{ij}=\frac1n\omega_n^{-ij}$. L'elemento $(i,j)$ di $GF_n$ è $$\sum_{k=0}^{n-1}G_{ik}[F_n]_{kj}=\frac1n\sum_{k=0}^{n-1}\omega_n^{-ik}\omega_n^{kj}=\frac1n\sum_{k=0}^{n-1}\bigl(\omega_n^{\,j-i}\bigr)^k.$$ Se $i=j$ la ragione è $\omega_n^0=1$ e la somma vale $n$: l'elemento è $\frac1n\cdot n=1$. Se $i\ne j$, allora $j-i\in[-(n-1),n-1]\setminus\{0\}$, che non contiene multipli di $n$: per il lemma di somma la somma è $0$. Quindi $GF_n=I_n$, e per matrici quadrate una inversa sinistra è anche destra (Matrice inversaL'inversa di una matrice quadrata A è la matrice A⁻¹ con A A⁻¹ = A⁻¹ A = I; esiste se e solo se rango(A) = n e si calcola con Gauss-Jordan riducendo (A | I) fino a (I | A⁻¹).Matrice inversa →). $\square$

$F_n^{-1}$ è quindi una Vandermonde sulle potenze negative $V(\omega_n^0,\omega_n^{-1},\dots,\omega_n^{-(n-1)})$ divisa per $n$ (le sue righe sono una permutazione di quelle di $F_n$). Esplicitamente $$a_i=\frac1n\sum_{k=0}^{n-1}y_k\,\omega_n^{-ik}=\frac1n\,q\bigl(\omega_n^{-i}\bigr),\qquad q(x)=\sum_ky_kx^k.$$ Cioè l'interpolazione da $\Omega_n$ è una valutazione del polinomio $q$ che ha come coefficienti i valori $\mathbf y$, nei punti $\bar\omega_n^{\,i}=\omega_n^{-i}$, seguita da una divisione per $n$.

Primo metodo: FFT sulle potenze negative. I lemmi di cancellazione, dimezzamento e somma valgono anche con $\bar\omega_n=\omega_n^{-1}$ come radice principale (le radici sono le stesse), quindi la FFT funziona identica con $\bar\omega_n$ al posto di $\omega_n$:

FFT_NEG(y)
n <- y.len
if n = 1 then return y
y0 <- (y[0], y[2], ..., y[n-2]);  y1 <- (y[1], y[3], ..., y[n-1])
a0 <- FFT_NEG(y0);  a1 <- FFT_NEG(y1)
ON_NEG <- e^(-i 2 pi / n)                   { = w_n^(-1) }
ONI <- 1
for i <- 0 to n/2 - 1
    t <- ONI * a1[i]
    a[i] <- a0[i] + t
    a[i + n/2] <- a0[i] - t
    ONI <- ONI * ON_NEG
return a

INV_FFT(y)
n <- y.len
a <- FFT_NEG(y)
for i <- 0 to n-1
    a[i] <- a[i] / n
return a

La divisione per $n$ va fatta fuori dalla ricorsione, una volta sola: se fosse dentro, ogni livello dividerebbe per la propria taglia e si accumulerebbe il fattore sbagliato. Il costo è $T_{\text{INV}}(n)=2n\log_2n+n=\Theta(n\log n)$.

Secondo metodo: stessa FFT e rovesciamento. Poiché $\bar\omega_n^{\,i}=\omega_n^{-i}=\omega_n^{(n-i)\bmod n}$, la valutazione di $q$ nei punti $\omega_n^{-i}$ è la DFT di $\mathbf y$ letta in ordine rovesciato: $$\bigl[\mathrm{DFT}^{-1}_n(\mathbf y)\bigr]_i=\frac1n\,\bigl[\mathrm{DFT}_n(\mathbf y)\bigr]_{(n-i)\bmod n},\quad\text{cioè}\quad\mathrm{DFT}^{-1}_n(\mathbf y)=\frac1n\,\mathrm{DFT}_n(\mathbf y)^R,$$ dove per un vettore $\mathbf x$ si pone $\mathbf x^R=(x_0,x_{n-1},x_{n-2},\dots,x_1)$ (l'indice $0$ resta al suo posto, gli altri si ribaltano).

INV_FFT(y)
n <- y.len
a' <- FFT(y)
a[0] <- a'[0] / n
for i <- 1 to n-1
    a[i] <- a'[n-i] / n
return a

Esempio. $\mathbf y=(18,-2-2i,-2,-2+2i)$ (dall'esempio di prima). $\mathrm{FFT}(\mathbf y)=(12,24,20,16)$: per esempio $z_1=18+i(-2-2i)+(-1)(-2)+(-i)(-2+2i)=18+(2-2i)+2+(2+2i)=24$. Rovesciando: $(12,16,20,24)$; dividendo per $4$: $(3,4,5,6)$, i coefficienti originali.

Prodotto di polinomi con la FFT

LIN_CONV(a, b)
n <- a.len                                   { n potenza di due }
ya <- FFT((a | 0_n));  yb <- FFT((b | 0_n))  { valutazione su Omega_2n }
for i <- 0 to 2n-1
    yc[i] <- ya[i] * yb[i]
return INV_FFT(yc)                           { restituisce (c | 0) }

Costo esatto. Tre trasformate di ordine $2n$, ciascuna $2\cdot2n\log_2(2n)=4n\log_2(2n)$ operazioni; in più $2n$ prodotti $\odot$ e $2n$ moltiplicazioni per $\frac1{2n}$ nell'inversa: $$T_{\text{LIN\_CONV}}(n)=12\,n\log_2(2n)+4n=\Theta(n\log n).$$ Verificato con un'implementazione che conta le operazioni, per $n=2^0,\dots,2^7$, e confrontando il risultato con numpy.convolve.

Per la convoluzione per definizione il costo è $2n^2-2n+1$. Il confronto per qualche valore:

Per piccoli $n$ la definizione è più economica (costanti più basse): è lo stesso fenomeno che motiva gli algoritmi ibridi (Algoritmo di Strassen e ibridazioneStrassen calcola il prodotto n x n con 7 chiamate ricorsive invece di 8, usando 10 somme/sottrazioni nel divide e 8 nel conquer: T(n) = 7T(n/2) + (9/2)n^2, Theta(n^(log2 7)) con soluzione esatta 7n^(log2 7) - 6n^2; la costante 7 è più che tripla di quella dell'iterativo (2) e Strassen batte l'iterativo solo da n = 655 (per le potenze di 2, da 1024). L'ibridazione cambia il caso base (iterativo per n <= n0): costo n^(log2 7)(2n0+5)/n0^(log2 7 - 2) - 6n^2, minimo per n0 = 8 con costante ~3.92.Algoritmo di Strassen e ibridazione →).

Taglie generiche. Se $n$ non è una potenza di due basta portare entrambi i vettori alla lunghezza $2^{\lceil\log_2n\rceil+1}$, la prima potenza di due almeno pari a $2n$ (e minore di $4n$). Per vettori di lunghezze diverse $m$ e $n$ si porta tutto alla prima potenza di due $\ge m+n-1$. Il costo resta $\Theta(n\log n)$.

Nota numerica. Con numeri in virgola mobile le radici $\omega_n^i$ sono approssimate e il risultato ha piccoli errori (ad esempio parti immaginarie dell'ordine di $10^{-15}$ al posto di $0$); con coefficienti interi si arrotonda. Un uso diverso della stessa idea, per convolvere segnali molto lunghi a blocchi, è in Convoluzione a blocchi - overlap-add e overlap-saveLa convoluzione lineare y = x * h (h lungo M) si può calcolare con la FFT scegliendo N >= Dx + M - 1 (potenza di 2): costo per campione circa 3 k log2 N contro k M del calcolo diretto, quindi la FFT conviene per M abbastanza grande (con Dx circa M, già da M = 16). Per segnali lunghi o infiniti si divide x in blocchi. Overlap-add: blocchi disgiunti di lunghezza L, ognuno convoluto con h (lunghezza L+M-1 con FFT di N >= L+M-1) e le code di M-1 campioni si sommano. Overlap-save: blocchi di N campioni sovrapposti di M-1, convoluzione circolare, si scartano i primi M-1 campioni di ogni uscita (corrotti dall'aliasing). Stesso costo.Convoluzione a blocchi - overlap-add e overlap-save →.

Domande d'esame

• Scrivere e risolvere la ricorrenza del numero di operazioni tra complessi della FFT, contando 4 operazioni per ogni iterazione del ciclo di ricombinazione. Traccia: $T(1)=0$, $T(n)=2T(n/2)+2n$, soluzione $T(n)=2n\log_2n$, oppure caso 2 del Master Theorem.
• Perché la divisione per $n$ della trasformata inversa si esegue fuori dalla ricorsione? Traccia: l'inversa è $\frac1n$ per la valutazione sulle potenze negative; la FFT ricorsiva su $\bar\omega$ non ha il fattore, che va applicato una sola volta al risultato finale.
• Dato $\mathrm{DFT}_4(\mathbf a)=(18,-2-2i,-2,-2+2i)$, ricavare $\mathbf a$ usando soltanto la FFT diretta. Traccia: $\mathrm{DFT}_4(\mathbf y)=(12,24,20,16)$, si rovescia in $(12,16,20,24)$ e si divide per $4$: $\mathbf a=(3,4,5,6)$.

Esercizi svolti: Esercizio - matrice di Fourier per una matrice generica a sinistra e a destra e per aI + bU (compito del 27 gennaio 2025, es. 2); Esercizio - convoluzione con il vettore delle potenze di due in tempo lineare (compito del 5 settembre 2025, es. 2); Esercizio - domande di prima parte su polinomi, DFT e FFT (simulazioni d'esame).