Esercizio - laboratorio 3, DFT di segnali periodici

Testo (laboratorio 3 del corso Teoria dei Segnali, UniPD, file lab3.m, parti 6-9). (6) Un segnale discreto periodico ha $N=12$ campioni per periodo, passo $T=5$ s (periodo $T_p=60$ s) ed è la somma di due esponenziali complessi, di frequenze $f_1=\frac1{T_p}$ e $f_2=\frac3{T_p}$ e ampiezze $1{,}5$ e $0{,}8$: rappresentarlo su 6 periodi, calcolarne la trasformata su un periodo con $\frac1N\,\text{fft}$ e rappresentarla su 4 periodi. (7) Ricostruire il segnale con l'antitrasformata. (8) Centrare lo spettro. (9) Scrivere la DFT a mano con due cicli annidati e confrontarla con la FFT per $N=1000$: tempo di calcolo ed errore.

Teoria usata: 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 →, Segnali a tempo discretoUn segnale a tempo discreto è una funzione complessa $s(nT)$ definita sui multipli interi del quanto temporale $T$ (insieme $\mathbb Z(T)$, velocità $F_p=1/T$). Le definizioni sono quelle dei segnali continui con la somma al posto dell'integrale e il quanto $T$ al posto di $dt$: area $\sum T,s(nT)$, energia $\sum T|s(nT)|^2$, convoluzione $\sum T,x(kT)y(nT-kT)$. L'impulso ideale discreto vale $1/T$ nell'origine. Esponenziali e sinusoidi discreti sono periodici solo se $f_0/F_p$ è razionale e hanno frequenza ambigua a meno di multipli di $F_p$. I segnali periodici con periodo $NT$ sono descritti da $N$ valori e si trattano al calcolatore.Segnali a tempo discreto →, Trasformata di Fourier a tempo discretoLa trasformata di Fourier di un segnale discreto è $S(f)=\sum_nT,s(nT),e^{-i2\pi fnT}$: una funzione continua e periodica in $f$ di periodo $F_p=1/T$. L'antitrasformata è l'integrale su un periodo, $s(nT)=\int_0^{F_p}S(f)e^{i2\pi fnT}df$. Le regole sono quelle del caso continuo (traslazione, convoluzione $\leftrightarrow$ prodotto, Parseval $\sum T|s|^2=\int_0^{F_p}|S|^2$), con incremento e somma corrente al posto di derivata e integrale. Per segnali reali $S(f)=S^(-f)$, quindi basta $[0,F_p/2]$. Esempi: $\delta\to1$, $1\to\delta_{F_p}$, rect $\to$ sinc periodico, $a^n\mathbf 1_0\to\frac T{1-ae^{-i2\pi fT}}$.Trasformata di Fourier a tempo discreto →, Serie di FourierUn segnale periodico di periodo $T_p$ si scrive come somma di esponenziali alle frequenze multiple della fondamentale $F=1/T_p$: $s(t)=\sum_nS_ne^{i2\pi nFt}$, con $S_n=\frac1{T_p}\int_{T_p}s(t)e^{-i2\pi nFt}dt$. Si basa sull'ortogonalità degli esponenziali su un periodo. $S_0$ è il valor medio; per segnali reali $S_{-n}=S_n^$ e si passa alla forma con coseni e seni; vale il teorema di Parseval $P=\sum|S_n|^2$. La convoluzione ciclica diventa il prodotto $T_pX_nY_n$ dei coefficienti. Le somme troncate presentano il fenomeno di Gibbs vicino ai salti.Serie di Fourier → (i coefficienti $S_k$), Complessità computazionaleCosto di un algoritmo in funzione della dimensione dell'input; caso peggiore, migliore e medio; notazione O-grande e classi di crescita; come contare i passi di cicli e ricorsioni; costo delle operazioni Python.Complessità computazionale → (parte 9).

Cosa si vuole calcolare e perché

Un segnale discreto periodico di periodo $T_p=NT$ è descritto da $N$ numeri (i campioni di un periodo); il suo spettro, per il corso, è fatto da righe alle frequenze $kF$ con $F=\frac1{T_p}$, ed è a sua volta periodico di periodo $F_p=\frac1T=NF$. La trasformata di un periodo è la DFT (discrete Fourier transform), e la FFT è l'algoritmo veloce che la calcola. Tutta la parte 6-8 serve a vedere questa doppia periodicità (periodico nel tempo $\leftrightarrow$ righe in frequenza; campionato nel tempo $\leftrightarrow$ periodico in frequenza); la parte 9 a capire perché si usa la FFT.

Le convenzioni: quale scalatura è quale

Le definizioni del corso per un segnale discreto periodico (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 →) sono

$$S(kF)=\sum_{n=0}^{N-1}T\,s(nT)\,e^{-i2\pi kFnT},\qquad s(nT)=\sum_{k=0}^{N-1}F\,S(kF)\,e^{i2\pi kFnT},$$

con $F\,T=\frac1N$ e $e^{-i2\pi kFnT}=e^{-i2\pi kn/N}$. Il comando np.fft.fft(s) calcola invece la somma senza fattori, $X_k=\sum_ns(nT)\,e^{-i2\pi kn/N}$. Quindi:

Il laboratorio usa la terza: $(1/N)\,\text{fft}$ dà i coefficienti di Fourier $S_k$ del segnale periodico, cioè le ampiezze dei singoli esponenziali nella sintesi $s(nT)=\sum_kS_k\,e^{i2\pi kFnT}$ (sono i $S_k$ della serie di Fourier, Serie di FourierUn segnale periodico di periodo $T_p$ si scrive come somma di esponenziali alle frequenze multiple della fondamentale $F=1/T_p$: $s(t)=\sum_nS_ne^{i2\pi nFt}$, con $S_n=\frac1{T_p}\int_{T_p}s(t)e^{-i2\pi nFt}dt$. Si basa sull'ortogonalità degli esponenziali su un periodo. $S_0$ è il valor medio; per segnali reali $S_{-n}=S_n^*$ e si passa alla forma con coseni e seni; vale il teorema di Parseval $P=\sum|S_n|^2$. La convoluzione ciclica diventa il prodotto $T_pX_nY_n$ dei coefficienti. Le somme troncate presentano il fenomeno di Gibbs vicino ai salti.Serie di Fourier →, applicati ai campioni). Lo si riconosce da come viene antitrasformata: N * ifft(S_p), perché ifft contiene già il fattore $\frac1N$ ($\text{ifft}(Y)_n=\frac1N\sum_kY_ke^{i2\pi kn/N}$) e per la sintesi $\sum_kS_ke^{\dots}$ serve la somma senza $\frac1N$. La trasformata $S(kF)$ del corso è $T_p$ volte più grande ($S(kF)=S_k/F=T_pS_k$): in questo esempio $S(F)=60\cdot1{,}5=90$ contro $S_1=1{,}5$. Il rapporto dà l'area dell'impulso: nella trasformata di un segnale periodico $\sum_kS_k\,\delta(f-kF)$ gli impulsi hanno area $S_k$, non $S(kF)$.

Parte 6 - il segnale e la sua trasformata su un periodo

import time as orologio
import numpy as np
import matplotlib.pyplot as plt

# PARTE 6 - segnale discreto periodico, N = 12 campioni per periodo, con due componenti
T = 5                # passo di campionamento [s]
N = 12               # campioni in un periodo
Tp = N * T           # periodo del segnale: 60 s

nTp = np.arange(N) * T            # istanti di UN periodo: 0, 5, ..., 55

f1 = 1 / Tp          # prima componente: 1 giro per periodo
f2 = 3 / Tp          # seconda: 3 giri per periodo

A1 = 1.5
A2 = 0.8

# un periodo di segnale (complesso)
s_p = A1 * np.exp(1j * 2 * np.pi * f1 * nTp) + A2 * np.exp(1j * 2 * np.pi * f2 * nTp)

# periodicita' nel tempo: si ripete il periodo
n_periodi_time = 6
s = np.tile(s_p, n_periodi_time)                 # MATLAB repmat(s_p, 1, 6)
nT = np.arange(len(s)) * T                       # 72 istanti, da 0 a 355

plt.figure()
plt.subplot(2, 1, 1)
plt.stem(nT, np.real(s), basefmt=' '); plt.grid(True); plt.xlabel('nT'); plt.ylabel('Re[s(nT)]')
plt.subplot(2, 1, 2)
plt.stem(nT, np.imag(s), basefmt=' '); plt.grid(True); plt.xlabel('nT'); plt.ylabel('Im[s(nT)]')
plt.show()

F = 1 / Tp           # passo in frequenza: 1/60 Hz
Fp = 1 / T           # periodo in frequenza (frequenza di campionamento): 0.2 Hz = N F

# un periodo di trasformata: (1/N) fft
S_p = (1 / N) * np.fft.fft(s_p)

# periodicita' in frequenza: si ripete il periodo della trasformata
n_periodi_freq = 4
S = np.tile(S_p, n_periodi_freq)
kF = np.arange(len(S)) * F                       # 48 frequenze, da 0 a 47 F

plt.figure()
plt.subplot(2, 1, 1)
plt.stem(kF, np.abs(S), linefmt='k-', markerfmt='k^', basefmt=' '); plt.grid(True); plt.xlabel('kF'); plt.ylabel('Abs[S(kF)]')
plt.subplot(2, 1, 2)
plt.stem(kF, np.angle(S), linefmt='m-', markerfmt='m^', basefmt=' '); plt.grid(True); plt.xlabel('kF'); plt.ylabel('Phase[S(kF)]')
plt.show()

print(f"Tp = {Tp} s, F = 1/Tp = {F:.5f} Hz, Fp = 1/T = {Fp} Hz = N F = {N * F:.1f} Hz")
print("modulo di (1/N) fft per k = 0..11:", np.round(np.abs(S_p), 6))
print("S_p[1] =", np.round(S_p[1], 6), "   S_p[3] =", np.round(S_p[3], 6))
print("fase delle righe non nulle:", np.round(np.angle(S_p[[1, 3]]), 6))
print("fase delle altre righe (modulo ~1e-17: e' solo rumore numerico):", np.round(np.angle(S_p[[0, 2, 4]]), 3))
print(f"numero di campioni: s = {len(s)}, S = {len(S)}")

Risultato:

Tp = 60 s, F = 1/Tp = 0.01667 Hz, Fp = 1/T = 0.2 Hz = N F = 0.2 Hz
modulo di (1/N) fft per k = 0..11: [0.  1.5 0.  0.8 0.  0.  0.  0.  0.  0.  0.  0. ]
S_p[1] = (1.5-0j)    S_p[3] = (0.8-0j)
fase delle righe non nulle: [-0. -0.]
fase delle altre righe (modulo ~1e-17: e' solo rumore numerico): [ 3.132 -2.232 -0.109]
numero di campioni: s = 72, S = 48

Il segnale. $s(nT)=1{,}5\,e^{i2\pi\,n/12}+0{,}8\,e^{i2\pi\,3n/12}$ perché $f_1nT=\frac{nT}{T_p}=\frac n{12}$ e $f_2nT=\frac{3n}{12}$: il primo termine compie un giro in 12 campioni, il secondo tre giri. Per questo l'intero segnale è periodico di periodo $12$ campioni ($60$ s). Il vettore nTp contiene gli istanti di un solo periodo ($0,5,\dots,55$), np.tile(s_p, 6) lo ripete 6 volte (72 campioni, come repmat di MATLAB). Parte reale: $1{,}5\cos\frac{2\pi n}{12}+0{,}8\cos\frac{2\pi\cdot3n}{12}$.

Lo spettro. Il modulo di $(1/N)\text{fft}$ mostra due righe sole: $S_1=1{,}5$ e $S_3=0{,}8$, tutte le altre sono zero (all'ordine di $10^{-16}$). Era prevedibile: il segnale è proprio una somma di due esponenziali, quindi i suoi coefficienti sono le loro ampiezze, nelle posizioni $k=1$ e $k=3$, cioè alle frequenze $F=\frac1{60}$ Hz e $3F=\frac1{20}$ Hz. Non ci sono righe alle frequenze negative $-F,-3F$ perché il segnale è complesso (un'esponenziale complessa ha una sola riga; un coseno reale ne ha due, in $\pm f_0$, come si vedrà nel laboratorio 4).

Le due periodicità. Nel tempo il segnale si ripete ogni $T_p=60$ s ($N=12$ campioni). In frequenza la trasformata si ripete con periodo $F_p=\frac1T=0{,}2$ Hz $=NF$, cioè ogni $N=12$ righe: np.tile(S_p, 4) ripete $S_p$ quattro volte (48 righe, da $0$ a $47F$). È la proprietà della trasformata dei segnali discreti (Trasformata di Fourier a tempo discretoLa trasformata di Fourier di un segnale discreto è $S(f)=\sum_nT,s(nT),e^{-i2\pi fnT}$: una funzione continua e periodica in $f$ di periodo $F_p=1/T$. L'antitrasformata è l'integrale su un periodo, $s(nT)=\int_0^{F_p}S(f)e^{i2\pi fnT}df$. Le regole sono quelle del caso continuo (traslazione, convoluzione $\leftrightarrow$ prodotto, Parseval $\sum T|s|^2=\int_0^{F_p}|S|^2$), con incremento e somma corrente al posto di derivata e integrale. Per segnali reali $S(f)=S^*(-f)$, quindi basta $[0,F_p/2]$. Esempi: $\delta\to1$, $1\to\delta_{F_p}$, rect $\to$ sinc periodico, $a^n\mathbf 1_0\to\frac T{1-ae^{-i2\pi fT}}$.Trasformata di Fourier a tempo discreto →): campionare nel tempo con passo $T$ ripete lo spettro con periodo $\frac1T$.

La fase. La fase delle due righe è $0$ (i coefficienti sono reali e positivi, $1{,}5$ e $0{,}8$, perché i due esponenziali partono da $n=0$ con fase nulla). Nelle altre righe il modulo è $\sim10^{-17}$ e np.angle restituisce valori qualunque ($3{,}13$; $-2{,}23$; $-0{,}11$): la fase è un'informazione priva di senso dove il modulo è nullo, e nel grafico di MATLAB le righe nulle mostrano fasi casuali. In pratica si traccia la fase solo dove $\lvert S_k\rvert$ supera una soglia.

La tabella delle scalature si verifica con questi numeri:

# Le tre "normalizzazioni" della stessa fft (convenzione del corso: S(kF) = somma di T s(nT) e^(-i 2 pi kF nT))
X = np.fft.fft(s_p)                      # somma di s(nT) e^(-i 2 pi k n / N), senza fattori
S_corso = T * X                          # S(kF) del corso
S_coeff = X / N                          # (1/N) fft = coefficienti S_k = F S(kF)
print("fft(s_p)[1]        =", np.round(X[1], 6))
print("T fft  (S(kF))[1]  =", np.round(S_corso[1], 6), "  = Tp * 1.5 =", Tp * 1.5)
print("fft/N  (S_k)[1]    =", np.round(S_coeff[1], 6))
print("F * S(kF) = S_k ?  ", np.allclose(F * S_corso, S_coeff))

fft(s_p)[1]        = (18-0j)
T fft  (S(kF))[1]  = (90-0j)   = Tp * 1.5 = 90.0
fft/N  (S_k)[1]    = (1.5-0j)
F * S(kF) = S_k ?   True

La somma grezza $X_1=\sum_ns\,e^{-i2\pi n/12}=12\cdot1{,}5=18$ (12 campioni, ciascuno dà la sua ampiezza); $T X_1=90=T_p\cdot1{,}5$ è il $S(F)$ del corso; $\frac{X_1}{12}=1{,}5$ è $S_1$; $F\,S(kF)=S_k$ è verificato su tutti i $k$.

Parte 7 - l'antitrasformata

# PARTE 7 - antitrasformata: N * ifft ricostruisce un periodo di segnale
s_p_inv = N * np.fft.ifft(S_p)

plt.figure()
plt.subplot(2, 1, 1)
plt.stem(nT, np.real(np.tile(s_p_inv, n_periodi_time)), basefmt=' '); plt.grid(True); plt.xlabel('nT'); plt.ylabel('Re[sinv(nT)]')
plt.subplot(2, 1, 2)
plt.stem(nT, np.imag(np.tile(s_p_inv, n_periodi_time)), basefmt=' '); plt.grid(True); plt.xlabel('nT'); plt.ylabel('Im[sinv(nT)]')
plt.show()

print(f"errore massimo |s_inv - s_p|: {np.max(np.abs(s_p_inv - s_p)):.2e}")

# la stessa cosa scritta con la formula del corso s(nT) = somma di S_k e^(i 2 pi k F nT)
s_sintesi = sum(S_p[k] * np.exp(1j * 2 * np.pi * k * F * nTp) for k in range(N))
print(f"errore della sintesi a mano: {np.max(np.abs(s_sintesi - s_p)):.2e}")

Risultato:

errore massimo |s_inv - s_p|: 4.44e-16
errore della sintesi a mano: 2.29e-15

N * np.fft.ifft(S_p) riporta al segnale: l'errore massimo è di $4\cdot10^{-16}$, cioè gli arrotondamenti del calcolatore. Perché il fattore $N$: la sintesi del corso (per i coefficienti $S_k$) è $s(nT)=\sum_{k=0}^{N-1}S_k\,e^{i2\pi kn/N}$, somma senza fattori; ifft calcola $\frac1N\sum_k$, quindi si moltiplica per $N$. La seconda riga del blocco fa la somma "a mano" (con le righe di $S_p$ e le esponenziali) e ritrova lo stesso segnale. I grafici (parte reale e immaginaria) sono quelli della parte 6, ripetuti su 6 periodi.

Parte 8 - centrare lo spettro

# PARTE 8 - centrare lo spettro
# (a) come nel laboratorio: stessi valori di S, nuove etichette dell'asse, da -24 F a 23 F
kF_centrato = np.arange(-len(S) // 2, len(S) // 2) * F
print("lab: S non e' stato riordinato. Vale perche' S e' periodico con periodo N = 12 bin e 24 = 2 periodi:",
      np.allclose(np.roll(S, 24), S))

plt.figure()
plt.subplot(2, 1, 1)
plt.stem(kF_centrato, np.abs(S), linefmt='k-', markerfmt='k^', basefmt=' '); plt.grid(True); plt.xlabel('kF'); plt.ylabel('Abs[S(kF)]')
plt.subplot(2, 1, 2)
plt.stem(kF_centrato, np.angle(S), linefmt='m-', markerfmt='m^', basefmt=' '); plt.grid(True); plt.xlabel('kF'); plt.ylabel('Phase[S(kF)]')
plt.show()

# (b) un solo periodo, centrato davvero: si riordinano le due meta' con fftshift
S_c = np.fft.fftshift(S_p)                         # [S_p(6..11), S_p(0..5)]
k_c = np.arange(-N // 2, N // 2)                   # -6, ..., 5
print("indici k dopo il centramento:", k_c)
print("modulo di S_c:", np.round(np.abs(S_c), 3))
print("righe nonnulle in:", [f"{k} F = {k * F:.4f} Hz" for k in k_c[np.abs(S_c) > 1e-9]])

Risultato:

lab: S non e' stato riordinato. Vale perche' S e' periodico con periodo N = 12 bin e 24 = 2 periodi: True
indici k dopo il centramento: [-6 -5 -4 -3 -2 -1  0  1  2  3  4  5]
modulo di S_c: [0.  0.  0.  0.  0.  0.  0.  1.5 0.  0.8 0.  0. ]
righe nonnulle in: ['1 F = 0.0167 Hz', '3 F = 0.0500 Hz']

Lo spettro di un segnale discreto è periodico, quindi "l'asse delle frequenze" è una scelta: si può guardare l'intervallo $[0,F_p)$ come fa fft (le frequenze "negative" compaiono alla fine, come $f-F_p$) oppure l'intervallo $[-\frac{F_p}2,\frac{F_p}2)$, più naturale.

• Modo del laboratorio. kF = (-len(S)/2 : len(S)/2-1) * F cambia solo le etichette dell'asse, da $-24F$ a $23F$, senza spostare i valori di S. È corretto soltanto perché S è fatta di 4 periodi esatti e lo scorrimento di $24$ posizioni è un multiplo di $N=12$: spostare un vettore periodico di un multiplo del periodo non lo cambia (verificato dal True stampato). Con un solo periodo non funzionerebbe.
• Modo "giusto" per un periodo. Si riordinano le due metà del vettore con np.fft.fftshift(S_p): le righe da $k=6$ a $11$ (che sono in realtà $k=-6,\dots,-1$, per la periodicità) vengono portate davanti, e l'asse diventa $k=-N/2,\dots,N/2-1$ (qui $-6,\dots,5$). Le due righe stanno in $k=1$ e $k=3$, cioè $+F=0{,}0167$ Hz e $+3F=0{,}05$ Hz: lo spettro centrato mostra le righe alle frequenze positive piccole, come atteso. Con $N$ pari la riga $k=-N/2$ è sia $-\frac{F_p}2$ che $+\frac{F_p}2$ (stessa riga per periodicità): per convenzione compare a sinistra.

Parte 9 - la DFT scritta a mano e la FFT

# PARTE 9 - la DFT scritta a mano, con due cicli, contro la FFT
def dft_manuale(x):
    """DFT con il fattore 1/N: X[k] = (1/N) somma di x[n] e^(-i 2 pi k n / N)."""
    N = len(x)
    X_dft = np.zeros(N, dtype=complex)
    for k in range(N):
        for n in range(N):
            X_dft[k] = X_dft[k] + (1 / N) * x[n] * np.exp(-1j * 2 * np.pi * k * n / N)
    return X_dft


N = 1000
T = 10
Tp = N * T
F = 1 / Tp

rng = np.random.default_rng(0)
x = 2 * rng.random(N) - 1 + 1j * (2 * rng.random(N) - 1)       # N campioni complessi casuali

t0 = orologio.perf_counter()
X_dft = dft_manuale(x)
t_dft = orologio.perf_counter() - t0                            # tempo DFT manuale

t0 = orologio.perf_counter()
X_fft = (1 / N) * np.fft.fft(x)
t_fft = orologio.perf_counter() - t0                            # tempo fft

print(f"Tempo DFT manuale: {t_dft:.4f} secondi")
print(f"Tempo FFT        : {t_fft:.6f} secondi")
errore = np.linalg.norm(X_dft - X_fft)                          # norma L2 della differenza
print(f"Errore tra dft manuale e fft: {errore:.2e}")
print(f"rapporto dei tempi: circa {t_dft / t_fft:.0f}")

Risultato (i tempi variano da un'esecuzione all'altra e da un computer all'altro):

Tempo DFT manuale: 0.3962 secondi
Tempo FFT        : 0.000166 secondi
Errore tra dft manuale e fft: 1.60e-13
rapporto dei tempi: circa 2384

La funzione dft_manuale applica la definizione: per ognuno degli $N$ valori di $k$ si sommano $N$ termini $x_n\,e^{-i2\pi kn/N}/N$. Sono $N^2=10^6$ prodotti complessi, ognuno con una esponenziale, eseguiti dall'interprete Python. La FFT calcola gli stessi $N$ numeri con circa $N\log_2N\approx10^4$ operazioni, riorganizzando la somma in modo ricorsivo (le DFT di lunghezza $N$ si spezzano in due di lunghezza $N/2$, e così via).

Errore. La differenza tra le due, in norma $L^2$, è $1{,}6\cdot10^{-13}$: sono lo stesso calcolo, differiscono solo per gli arrotondamenti. Il $\frac1N$ nel dft_manuale è lo stesso di (1/N) * fft(x).

Tempo. La DFT a mano è migliaia di volte più lenta (rapporto stampato sopra; la prima chiamata della FFT è più lenta delle successive, e con tempi ripetuti il rapporto è di decine di migliaia, perché la FFT impiega pochi microsecondi anche per $N=1000$). Il confronto giusto è come cresce il tempo con $N$:

# Quanto cresce il tempo con N: la DFT a mano e' O(N^2), la FFT e' O(N log N)
def tempo_medio(funzione, x, ripetizioni):
    t0 = orologio.perf_counter()
    for _ in range(ripetizioni):
        funzione(x)
    return (orologio.perf_counter() - t0) / ripetizioni

print("   N     N^2     N log2 N   DFT a mano [s]   FFT [s]")
precedente = None
for N in [125, 250, 500, 1000]:
    x = rng.random(N) + 1j * rng.random(N)
    t_man = tempo_medio(dft_manuale, x, 1)
    t_f = tempo_medio(np.fft.fft, x, 2000)
    print(f"{N:>5} {N**2:>8} {N * np.log2(N):>10.0f} {t_man:>14.4f} {t_f:>12.2e}"
          + ("" if precedente is None else f"   (x{t_man / precedente[0]:.1f} la DFT, x{t_f / precedente[1]:.1f} la FFT)"))
    precedente = (t_man, t_f)

   N     N^2     N log2 N   DFT a mano [s]   FFT [s]
  125    15625        871         0.0061     3.41e-06
  250    62500       1991         0.0246     4.22e-06   (x4.0 la DFT, x1.2 la FFT)
  500   250000       4483         0.0991     5.51e-06   (x4.0 la DFT, x1.3 la FFT)
 1000  1000000       9966         0.3957     8.57e-06   (x4.0 la DFT, x1.6 la FFT)

Raddoppiando $N$ la DFT a mano costa circa $4$ volte di più (cresce come $N^2$: i fattori in tabella sono vicini a $4$), la FFT molto meno ($N\log N$; per $N$ così piccoli domina il tempo fisso di una chiamata, e i fattori sono tra $1$ e $2$). Per un segnale di un milione di campioni $N^2=10^{12}$ operazioni contro $N\log_2N\approx2\cdot10^7$: la FFT rende possibile ciò che con la definizione sarebbe impraticabile (Complessità computazionaleCosto di un algoritmo in funzione della dimensione dell'input; caso peggiore, migliore e medio; notazione O-grande e classi di crescita; come contare i passi di cicli e ricorsioni; costo delle operazioni Python.Complessità computazionale →). Anche la DFT riscritta come prodotto di matrici (dft_matrice, nessun ciclo Python) fa sempre $N^2$ operazioni, soltanto più veloci perché fatte in C:

# Versione "a matrice" (senza cicli Python): sempre N^2 prodotti, ma fatti in C
def dft_matrice(x):
    N = len(x)
    k = np.arange(N).reshape(-1, 1)
    n = np.arange(N).reshape(1, -1)
    W = np.exp(-2j * np.pi * k * n / N)         # matrice N x N dei fattori e^(-i 2 pi k n / N)
    return (W @ x) / N

N = 1000
x = rng.random(N) + 1j * rng.random(N)
t0 = orologio.perf_counter()
X_mat = dft_matrice(x)
t_mat = orologio.perf_counter() - t0
print(f"DFT a matrice, N = 1000: {t_mat * 1000:.1f} ms;  errore rispetto alla fft: {np.linalg.norm(X_mat - np.fft.fft(x) / N):.2e}")

DFT a matrice, N = 1000: 29.0 ms;  errore rispetto alla fft: 1.25e-13

Qualche decina di millisecondi contro mezzo secondo del ciclo, ma sempre migliaia di volte più lenta della FFT.

Controllo

• $\frac1N$fft: righe $S_1=1{,}5$ e $S_3=0{,}8$, le altre $\approx10^{-16}$; $T\cdot$fft dà $90$ e $48$ ($T_p S_k$).
• Antitrasformata: errore $4\cdot10^{-16}$.
• Spettro periodico di periodo $F_p=\frac1T=0{,}2$ Hz $=12F$.
• DFT a mano contro FFT ($N=1000$): errore $1{,}6\cdot10^{-13}$; tempo $\propto N^2$ contro $N\log N$.