Esercizio - DFT di un segnale discreto periodico

Testo (dispense del corso Teoria dei Segnali, UniPD, §6.5 e §8.4; argomento della lezione 12). Un segnale discreto $s(nT)$ periodico di periodo $T_p=NT$ è assegnato dai suoi $N$ valori in un periodo. Si svolgono:

• (a) la DFT a mano per $N=4$ e per $N=8$ (matrice di Fourier, valori di $e^{-i2\pi/N}$), con l'inversione e la forma unitaria;
• (b) la DFT di un coseno campionato che ha un numero intero di periodi nella finestra e di uno che ne ha un numero non intero (perdita spettrale, leakage), con np.fft.fft e la convenzione del corso;
• (c) la relazione tra la DFT e i coefficienti della serie di Fourier di un segnale continuo periodico campionato (aliasing), con l'esempio dell'onda quadra;
• (d) il teorema di Parseval nella forma $\sum T|s|^2=\sum F|S|^2$;
• (e) il legame con la trasformata di un segnale discreto a durata limitata (zero padding).

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 →, 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 →, Campionamento e ricostruzioneIl campionamento $s_c(nT)=s(nT)$ trasforma un segnale continuo in uno discreto e, in frequenza, ripete lo spettro con periodo $F_c=1/T$: $S_c(f)=\sum_kS(f-kF_c)$. Se le repliche si sovrappongono si ha aliasing e il segnale non è recuperabile. Un interpolatore $\mathbb Z(T)\to\mathbb R$ con risposta impulsiva $g$ produce $\tilde s(t)=\sum_nT g(t-nT)s(nT)$ e $\tilde S=G,S_c$. Teorema del campionamento: se $s$ ha banda $B$ e $F_c\ge2B$, l'interpolatore ideale $G=\operatorname{rect}(f/F_c)$ ricostruisce esattamente $s(t)=\sum s(nT)\operatorname{sinc}(F_c(t-nT))$. Se le ipotesi non valgono c'è un errore, in banda e fuori banda, riducibile con un prefiltro anti-aliasing.Campionamento e ricostruzione →, Numeri complessiI numeri complessi estendono i reali introducendo l'unità immaginaria $i$ ($i^2 = -1$) e possono essere rappresentati in forma algebrica, trigonometrica o polare. Tramite la formula di Eulero e le proprietà del modulo e dell'argomento, è possibile calcolare agilmente prodotti, potenze e radici ennesime.Numeri complessi →, SommatorieIl simbolo di sommatoria, le sue proprietà (linearità, additività, cambio di indice) e le somme notevoli di Gauss e geometrica.Sommatorie →.

La DFT con la convenzione del corso

Sia $s(nT)$ periodico di periodo $T_p=NT$ (cioè $N$ campioni per periodo). Si pongono

$$F=\frac1{T_p}=\frac1{NT}\ \ (\text{quanto in frequenza}),\qquad F_p=\frac1T=NF\ \ (\text{periodo in frequenza}),\qquad FT=\frac1N .$$

(Nel §6.5 del testo il simbolo $F_p$ è usato anche per $\frac1{NT}$; qui si tiene $F$ per $\frac1{T_p}$ come nel §8.4.) La trasformata e la sua inversa (8.27) sono

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

(si è usato $kFnT=\frac{kn}N$). Sono entrambe periodiche di periodo $N$: $S((k+N)F)=S(kF)$ e $s((n+N)T)=s(nT)$, perché $e^{\pm i2\pi N n/N}=1$. Quindi $S$ è una sequenza di $N$ valori (alle frequenze $0,F,\dots,(N-1)F$) e occupa un periodo $F_p=NF$.

Perché l'inversa torna. Si sostituisce $S(kF)$ nella seconda formula: $$\sum_kF\sum_mT\,s(mT)e^{-i2\pi k(m-n)/N}=FT\sum_m s(mT)\sum_{k=0}^{N-1}e^{-i2\pi k(m-n)/N}.$$ La somma interna è geometrica (SommatorieIl simbolo di sommatoria, le sue proprietà (linearità, additività, cambio di indice) e le somme notevoli di Gauss e geometrica.Sommatorie →) di ragione $z=e^{-i2\pi(m-n)/N}$: vale $N$ se $m=n$ (con $m,n\in\{0,\dots,N-1\}$, $z=1$) e $\frac{1-z^N}{1-z}=0$ altrimenti (perché $z^N=1$, $z\ne1$). Resta $FT\cdot N\,s(nT)=s(nT)$, essendo $FT\,N=1$. ✓

Come si calcola con NumPy. np.fft.fft(s) calcola $X_k=\sum_ns[n]e^{-i2\pi kn/N}$ (niente $T$, niente $\frac1N$). Quindi, con la convenzione del corso:

$$S(kF)=T\cdot X_k=T\cdot\texttt{np.fft.fft(s)}[k],\qquad S_h=F\,S(hF)=\frac{X_h}N .$$

I coefficienti $S_h$ del §6.5. La trasformata ordinaria $S(f)$ del segnale periodico è un treno di impulsi alle frequenze $hF$; l'area dell'impulso in $hF$ è $S_h=F\,S(hF)$ e si ha $s(nT)=\sum_{h=0}^{N-1}S_h\,e^{i2\pi hn/N}$: sono i coefficienti della serie di Fourier discreta ($S_h=\frac1{T_p}\sum_nT\,s(nT)e^{-i2\pi hn/N}=\frac1N\sum_ns(nT)e^{-i2\pi hn/N}$, valor medio $S_0$, ampiezza della componente a frequenza $hF$). Attenzione a una discrepanza del testo: la (6.21) scrive $S_h=\frac1N\sum_{k}T\,s(kT)e^{-i2\pi hk/N}$ con un fattore $T$ in più. Nella derivazione si moltiplica per $F_p=\frac1{NT}$ che già contiene $\frac1T$: $F_p\,T\,s=\frac sN$. Prova: per $s(nT)\equiv1$ l'impulso di $S(f)$ in $0$ deve avere area $1$ (la trasformata del segnale costante è $\delta_{F_p}(f)$, (6.14)) e la formula con la $T$ darebbe $T$. Qui si usa la forma corretta $S_h=F\,S(hF)=\frac1N\sum_n s(nT)e^{-i2\pi hn/N}$, coerente con la (8.27) e con la (8.24).

Forma matriciale e forma unitaria. Con $W=e^{-i2\pi/N}$ (radice $N$-esima dell'unità, $W^N=1$) si mettano i campioni in un vettore $\mathbf s=(s(0),\dots,s((N-1)T))^T$ e i valori di $S$ in $\mathbf S$. Allora

$$\mathbf S=T\,\mathbf W\,\mathbf s,\qquad \mathbf W=[W^{kn}]_{k,n=0}^{N-1}\quad(\text{matrice di Fourier}),\qquad \mathbf s=F\,\mathbf W^H\mathbf S .$$

($\mathbf W^H$ è la trasposta coniugata, con elementi $W^{-kn}$.) La matrice $\mathbf U=\frac1{\sqrt N}\mathbf W$ è unitaria: $\mathbf U^H\mathbf U=I$ (è proprio l'ortogonalità $\sum_ke^{-i2\pi k(m-n)/N}=N\delta_{mn}$). Se si rinormalizzano i vettori, $\hat{\mathbf s}=\sqrt T\,\mathbf s$ e $\hat{\mathbf S}=\sqrt F\,\mathbf S$, la trasformata diventa la rotazione $\hat{\mathbf S}=\mathbf U\hat{\mathbf s}$ (infatti $\sqrt F\,T\,\mathbf W\mathbf s=\sqrt{FT}\,\mathbf W\,\hat{\mathbf s}=\frac1{\sqrt N}\mathbf W\hat{\mathbf s}$). Una matrice unitaria conserva la norma: $\|\hat{\mathbf S}\|=\|\hat{\mathbf s}\|$, cioè il teorema di Parseval (si veda (d)).

(a) DFT a mano

$N=4$, $T=\frac12$, $s=(1,\ 2,\ 0,\ -1)$

Dati. $T_p=NT=2$, $F=\frac1{T_p}=\frac12$ Hz, $F_p=NF=2$ Hz. Le frequenze sono $0,\ 0{,}5,\ 1,\ 1{,}5$ Hz.

La matrice. $W=e^{-i2\pi/4}=e^{-i\pi/2}=-i$. Le potenze si ripetono con periodo $4$: $W^0=1$, $W^1=-i$, $W^2=-1$, $W^3=i$, $W^4=1$. L'elemento $(k,n)$ è $W^{kn\bmod4}$:

$$\mathbf W=\begin{pmatrix}1&1&1&1\\1&-i&-1&i\\1&-1&1&-1\\1&i&-1&-i\end{pmatrix}.$$

Il calcolo $X=\mathbf W\mathbf s$ riga per riga (con $\mathbf s=(1,2,0,-1)^T$):

• $X_0=1+2+0-1=2$;
• $X_1=1\cdot1+(-i)\cdot2+(-1)\cdot0+i\cdot(-1)=1-2i-i=1-3i$;
• $X_2=1-2+0+1=0$;
• $X_3=1+i\cdot2+(-1)\cdot0+(-i)\cdot(-1)=1+2i+i=1+3i$.

Perciò $X=(2,\ 1-3i,\ 0,\ 1+3i)$ e, moltiplicando per $T=\frac12$,

$$S(kF)=\left(1,\ \tfrac12-\tfrac32i,\ 0,\ \tfrac12+\tfrac32i\right).$$

Il segnale è reale, quindi $X_{N-k}=\overline{X_k}$ ($X_3=\overline{X_1}$) e $X_0$, $X_2$ sono reali: bastano $X_0,X_1,X_2$. $S(0)=1$ è l'area di un periodo: $\sum_nT\,s=\frac12\cdot2=1$ ✓.

Inversione, $s(nT)=\sum_kF\,S(kF)\,e^{i2\pi kn/4}$, cioè $s=\frac14\sum_kX_ki^{kn}$ (perché $FT=\frac1N=\frac14$). Per $n=0$: $\frac14(2+1-3i+0+1+3i)=\frac44=1$ ✓. Per $n=1$ (i fattori sono $i^k=1,i,-1,-i$): $\frac14[2+(1-3i)i+0+(1+3i)(-i)]=\frac14[2+(i+3)+(-i+3)]=\frac84=2$ ✓. Anche $n=2$ ($i^{2k}=1,-1,1,-1$): $\frac14[2-(1-3i)+0-(1+3i)]=\frac04=0$ ✓ e $n=3$: $-1$ ✓.

Parseval. $\sum T|s|^2=\frac12(1+4+0+1)=3$ e $\sum F|S|^2=\frac12\left(1+\frac52+0+\frac52\right)=3$ ✓ (con $\left|\frac12-\frac32i\right|^2=\frac14+\frac94=\frac52$).

I coefficienti $S_h$ $=\frac{X_h}N=\left(\frac12,\ \frac14-\frac34i,\ 0,\ \frac14+\frac34i\right)$: $S_0=\frac12$ è il valor medio $\frac{1+2+0-1}4$.

$N=8$, $T=\frac18$, $s=(1,1,1,1,0,0,0,0)$

Dati. $T_p=1$ s, $F=1$ Hz, $F_p=8$ Hz. Il segnale vale $1$ per mezzo periodo e $0$ per l'altro mezzo (un'onda quadra campionata).

Le potenze di $W=e^{-i2\pi/8}=e^{-i\pi/4}=\frac{1-i}{\sqrt2}$:

con $\frac1{\sqrt2}=0{,}7071$.

Il calcolo. Solo $n=0,\dots,3$ hanno $s=1$, quindi $X_k=1+W^k+W^{2k}+W^{3k}$ (somma geometrica: $\frac{1-W^{4k}}{1-W^k}$ con $W^{4k}=(-1)^k$, quindi $X_k=0$ per $k$ pari non nullo e $X_k=\frac2{1-W^k}$ per $k$ dispari).

• $X_0=4$.
• $X_1=1+(0{,}7071-0{,}7071i)+(-i)+(-0{,}7071-0{,}7071i)=1-2{,}4142\,i$ (dove $2{,}4142=1+\sqrt2$).
• $X_2=1+(-i)+(-1)+(i)=0$, e anche $X_4=1-1+1-1=0$, $X_6=0$.
• $X_3=1+W^3+W^6+W^9=1+(-0{,}7071-0{,}7071i)+i+(0{,}7071-0{,}7071i)=1-0{,}4142\,i$ (con $W^9=W$).
• Per simmetria hermitiana: $X_5=\overline{X_3}=1+0{,}4142i$, $X_7=\overline{X_1}=1+2{,}4142i$.

Controllo. $\sum|X_k|^2=16+2\cdot(1+5{,}8284)+2\cdot(1+0{,}1716)=32=N\sum|s|^2=8\cdot4$; nelle unità del corso $\sum T|s|^2=\frac48=\frac12$ e $\sum F|S|^2=\frac{32}{64}=\frac12$ ($S(kF)=\frac{X_k}8$) ✓. Le armoniche pari mancano, come in ogni onda quadra a metà periodo.

(b) Un coseno campionato: periodo intero e non intero

Dati. $T=0{,}01$ s ($F_p=100$ Hz), $N=32$ campioni, quindi $T_p=0{,}32$ s e $F=\frac1{T_p}=3{,}125$ Hz. Il segnale è $s(nT)=\cos(2\pi f_0nT)$, $n=0,\dots,31$, in due casi:

• caso A: $f_0=4F=12{,}5$ Hz: nella finestra ci sono esattamente $4$ periodi del coseno, il segnale è periodico di periodo $NT$ (e la finestra è un periodo del segnale);
• caso B: $f_0=4{,}5F=14{,}0625$ Hz: nella finestra ci sono $4{,}5$ periodi. Il segnale discreto $\cos(2\pi f_0nT)$ è periodico di periodo $2NT$ ma non di periodo $NT$: ripetere la finestra di $32$ campioni crea uno scalino tra l'ultimo campione e il successivo.

Caso A. Per Eulero $s=\frac12e^{i2\pi\cdot4n/N}+\frac12e^{-i2\pi\cdot4n/N}$. Per l'ortogonalità (la somma geometrica di sopra) tutto si concentra in $k=4$ e $k=N-4=28$ (la frequenza $-4F$ ripiegata):

$$X_4=X_{28}=\frac N2=16,\qquad X_k=0\ \text{altrove}.$$

Con la convenzione del corso $S(4F)=S(28F)=T\frac N2=0{,}16$ (che è $\frac{T_p}2$), mentre i coefficienti sono $S_4=S_{28}=\frac{X_4}N=\frac12$: i due coefficienti $\frac12$ del coseno ($\cos=\frac12e^{+}+\frac12e^{-}$), come nella serie di Fourier di un coseno continuo. Tutta la potenza ($\frac12$, cioè $\frac14+\frac14$) è nei due bin.

Caso B. Ora $f_0=4{,}5F$ non coincide con nessun bin ($k$ intero). La DFT assume che i $32$ campioni siano un periodo; il segnale effettivamente periodico è la ripetizione della finestra, che non è il coseno originale. Il coseno "vero" spaccato in due esponenziali ha $\nu=4{,}5$ nel mezzo dei bin $4$ e $5$: ognuno dei due esponenziali dà un contributo che, nel bin $h$, ha modulo $\frac12\left|\frac{\sin\pi(h-\nu)}{N\sin(\pi(h-\nu)/N)}\right|$ (è la funzione $\operatorname{sinc}_N$ della (6.16)), che per $h-\nu=\pm\frac12$ vale circa $\frac1\pi=0{,}318$ e non $\frac12$, e decade lentamente (come $\frac1{|h-\nu|}$) lontano. I numeri ottenuti con np.fft.fft:

e la stessa cosa è simmetrica per $h\to N-h$ (segnale reale). Il massimo $|S_h|$ è $0{,}331$ invece di $0{,}5$ (ampiezza sottostimata di circa un terzo), e nei quattro bin centrali $4,5,27,28$ c'è solo l'$81{,}2\%$ della potenza (nel caso A il $100\%$). La parte restante sbava su tutti i bin: è la perdita spettrale (leakage). La somma dei quadrati resta $\frac12$ (Parseval vale sempre, vedi (d)): la potenza si ridistribuisce, non si perde.

Morale. La risoluzione in frequenza è $F=\frac1{NT}$: la si migliora osservando più a lungo (più campioni a parità di $T$), non campionando più in fretta. E se la finestra non contiene un numero intero di periodi la frequenza non è "su un bin" e compare il leakage (in pratica si attenua con una finestra di pesatura, argomento oltre questo esercizio).

(c) DFT e serie di Fourier di un segnale continuo periodico campionato

Il fatto generale. Sia $x(t)$ continuo e periodico di periodo $T_p$, con serie di Fourier $x(t)=\sum_{h\in\mathbb Z}X_h\,e^{i2\pi hFt}$, $F=\frac1{T_p}$ (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 →). Lo si campiona con $T=\frac{T_p}N$ (esattamente $N$ campioni per periodo): $s(nT)=x(nT)=\sum_hX_h\,e^{i2\pi hn/N}$, perché $FnT=\frac nN$. Gli esponenziali $e^{i2\pi hn/N}$ con $h$ che differisce per un multiplo di $N$ coincidono sui campioni ($e^{i2\pi kN n/N}=1$). Si raggruppano gli indici $h=r+kN$, $r=0,\dots,N-1$, $k\in\mathbb Z$:

$$s(nT)=\sum_{r=0}^{N-1}\left(\sum_{k=-\infty}^{+\infty}X_{r+kN}\right)e^{i2\pi rn/N}\qquad\Longrightarrow\qquad\boxed{S_r=\sum_{k=-\infty}^{+\infty}X_{r+kN}}$$

perché i coefficienti della rappresentazione di $s$ con $N$ esponenziali sono unici. I coefficienti della versione campionata (cioè $\frac{X^{\text{DFT}}_r}N$) sono la somma dei coefficienti del continuo a distanza $N$: le armoniche alte del continuo si "ripiegano" sulle basse. È lo stesso fenomeno del Teorema 6.1, $\operatorname{rep}_{F_p}$ di uno spettro a righe. Se $x$ ha banda limitata, $X_h=0$ per $|h|\ge N/2$, le righe non si sovrappongono e $S_r=X_r$ per $0\le r<N/2$ ($S_{N-r}=X_{-r}$): nessun aliasing (è la condizione di Nyquist $N>2h_{\max}$).

Esempio: onda quadra. $x(t)=+1$ per $|t|<\frac{T_p}4$ e $-1$ per $\frac{T_p}4<|t|\le\frac{T_p}2$, ripetuta con periodo $T_p$. Coefficienti: $X_0=0$ e, per $h\ne0$,

$$X_h=\frac1{T_p}\int_{-T_p/4}^{T_p/4}e^{-i2\pi hFt}dt-\frac1{T_p}\int_{T_p/4<|t|<T_p/2}e^{-i2\pi hFt}dt=\frac{2\sin(\pi h/2)}{\pi h}=\operatorname{sinc}\left(\frac h2\right),$$

cioè $X_{\pm1}=\frac2\pi=0{,}6366$, $X_{\pm3}=-\frac2{3\pi}=-0{,}2122$, $X_{\pm5}=\frac2{5\pi}$, e $X_h=0$ per $h$ pari. Decadono come $\frac1h$ (il segnale ha salti): l'aliasing sarà importante.

Campionamento con $N=6$. I campioni cadono in $t=\frac{n}{6}T_p$ per $n=0,\dots,5$ (cioè $n=-2,\dots,3$): i salti stanno in $t=\pm\frac{T_p}4=\pm1{,}5$ campioni, tra due campioni, quindi non si campiona mai su un salto. I valori sono $s=(1,\ 1,\ -1,\ -1,\ -1,\ 1)$ ($n=0$: $t=0$; $n=\pm1$: $\frac{T_p}6<\frac{T_p}4$, vale $+1$; $n=\pm2$ e $n=3$: $-1$). Il segnale è pari, quindi la DFT è reale:

$$S_r=\frac16\Bigl[s_0+2s_1\cos\tfrac{2\pi r}6+2s_2\cos\tfrac{4\pi r}6+s_3\cos\pi r\Bigr]=\frac16\Bigl[1+2\cos\tfrac{\pi r}3-2\cos\tfrac{2\pi r}3-(-1)^r\Bigr].$$

Valori: $r=0$: $\frac16(1+2-2-1)=0$; $r=1$: $\frac16(1+1+1+1)=\frac23$; $r=2$: $\frac16(1-1+1-1)=0$; $r=3$: $\frac16(1-2-2+1)=-\frac13$; poi $S_4=S_2=0$, $S_5=S_1=\frac23$. Confronto con le armoniche del continuo:

Non sono uguali: $S_1=X_1+X_{-5}+X_7+X_{-11}+X_{13}+\dots=0{,}6366+0{,}1273-0{,}0909-0{,}0579+0{,}0490+\dots$, e $S_3=X_3+X_{-3}+X_9+X_{-9}+\dots=-0{,}2122-0{,}2122+0{,}0707+0{,}0707-\dots$. I termini decadono come $\frac1h$ e la serie converge lentamente: a $4$ milioni di termini la somma è $0{,}66667$ per $S_1$ e $-0{,}33333$ per $S_3$, uguali alla DFT (verificato nel Controllo).

Più campioni, meno aliasing. L'errore decade come $\frac1{N^2}$ (le armoniche ripiegate partono da $h=N\pm r$ e i coefficienti scendono come $\frac1h$):

(per $N\equiv2\pmod4$ nessun campione cade su un salto). Con una forma d'onda continua senza salti i coefficienti scendono molto più in fretta e l'aliasing con pochi campioni è già trascurabile.

(d) Parseval nella forma $\sum T|s|^2=\sum F|S|^2$

Enunciato (8.29 per segnali discreti periodici, $T_p=NT$, $F=\frac1{T_p}$):

$$\sum_{n=0}^{N-1}T\,|s(nT)|^2=\sum_{k=0}^{N-1}F\,|S(kF)|^2 .$$

Dimostrazione con la forma unitaria. Se $\hat{\mathbf s}=\sqrt T\,\mathbf s$ e $\hat{\mathbf S}=\sqrt F\,\mathbf S$, allora $\hat{\mathbf S}=\mathbf U\hat{\mathbf s}$ con $\mathbf U$ unitaria, quindi $\|\hat{\mathbf S}\|^2=\hat{\mathbf s}^H\mathbf U^H\mathbf U\hat{\mathbf s}=\|\hat{\mathbf s}\|^2$. I due membri sono proprio $\|\hat{\mathbf s}\|^2$ e $\|\hat{\mathbf S}\|^2$. In termini di $X$: $\sum|X_k|^2=N\sum|s|^2$.

In termini di coefficienti. Con $S(kF)=T\,N\,S_k$ si ha $F|S(kF)|^2=\frac1{NT}T^2N^2|S_k|^2=TN|S_k|^2=T_p|S_k|^2$. Dividendo per $T_p$: la potenza $P=\frac1{T_p}\sum_nT|s|^2=\frac1N\sum_n|s|^2$ è $\sum_k|S_k|^2$ (Parseval della serie di Fourier, $P=\sum|S_k|^2$).

Verifiche sugli esempi precedenti.

• $N=4$, $T=\frac12$: $\sum T|s|^2=3$ e $\sum F|S|^2=3$. In potenza: $P=\frac3{T_p}=1{,}5$ e $\sum|S_h|^2=\frac14+\frac58+0+\frac58=1{,}5$ ✓ (con $|S_1|^2=\frac1{16}+\frac9{16}=\frac58$).
• $N=8$: $\frac12=\frac12$ (visto sopra).
• Coseno, caso A: $\sum T|s|^2=T\frac N2=0{,}16$ e $\sum F|S|^2=2\cdot F\cdot0{,}16^2=2\cdot3{,}125\cdot0{,}0256=0{,}16$ ✓; $P=\frac{0{,}16}{0{,}32}=\frac12=\left(\frac12\right)^2+\left(\frac12\right)^2$. Caso B: stessa energia $0{,}16$ (calcolata in entrambi i modi), anche se distribuita su tutti i bin.
• Onda quadra, $N=6$: $P=\frac16\sum1=1=\frac49+\frac19+\frac49$ ✓. Il continuo ha la stessa potenza ($P=1=\sum|X_h|^2$): l'aliasing non cambia la potenza, la ridistribuisce spostando le righe alte sulle basse.

(e) Segnale a durata limitata: la DFT campiona la trasformata

Se un segnale aperiodico $s_2$ (per esempio i tre campioni $A_0$ in $n=-1,0,1$ dell'Es. 6.3A, $S_2(f)=A_0T[1+2\cos(2\pi fT)]$) si estende per $L\le N$ campioni, lo si "periodicizza" con periodo $NT$ (aggiungendo $N-L$ zeri, zero padding) e si calcola la DFT: $S(kF)=\sum_{n}T\,s(nT)e^{-i2\pi kn/N}$ è esattamente la somma che definisce la trasformata $S_2(f)$ valutata in $f=kF$, perché per $f=kF$ si ha $e^{-i2\pi fnT}=e^{-i2\pi kn/N}$. Quindi la DFT dà i campioni della trasformata aperiodica nei punti $kF=\frac k{NT}$, e aumentare $N$ con zeri infittisce i punti. (Per i campioni in $n=-1,0,1$ basta ricordare che l'indice $n=-1$ equivale a $n=N-1$ per la periodicità.) Lo si vede nel Controllo con $N=8$ e $N=32$.

Controllo

import numpy as np
np.set_printoptions(precision=4, suppress=True, linewidth=120)

# Convenzione del corso per N campioni su un periodo Tp = N T:
#   S(kF) = somma_{n=0}^{N-1} T s(nT) e^{-i 2 pi k n/N}   =  T * np.fft.fft(s)[k]
#   s(nT) = somma_{k=0}^{N-1} F S(kF) e^{+i 2 pi k n/N}   (F = 1/(N T), F T = 1/N)
def dft_corso(s, T):
    return T * np.fft.fft(s)

def idft_corso(S, T):
    N = len(S); F = 1 / (N * T)
    n = np.arange(N); k = np.arange(N)
    return np.array([np.sum(F * S * np.exp(2j * np.pi * k * m / N)) for m in n])

# ---------- (a) N = 4, T = 1/2, s = (1, 2, 0, -1)
T, N = 0.5, 4
F = 1 / (N * T)
s = np.array([1, 2, 0, -1.0])
W = np.exp(-2j * np.pi / N)                       # = -i
M = np.array([[W**(k * n) for n in range(N)] for k in range(N)])
print(np.round(M, 3))
X = M @ s                                         # DFT "senza T"
S = T * X
print("X =", X, "  S(kF) =", S)
print("uguale a T*fft:", np.allclose(S, dft_corso(s, T)), " inversa:", np.allclose(idft_corso(S, T), s))
print("Parseval:", np.sum(T * abs(s)**2), np.sum(F * abs(S)**2))
U = M / np.sqrt(N)
print("U unitaria:", np.allclose(U.conj().T @ U, np.eye(N)), " S_norm = U s_norm:", np.allclose(U @ (np.sqrt(T) * s), np.sqrt(F) * S))

# N = 8, T = 1/8, s = (1,1,1,1,0,0,0,0)
T, N = 1 / 8, 8
s = np.array([1, 1, 1, 1, 0, 0, 0, 0.0])
W = np.exp(-2j * np.pi / N)
X = np.array([sum(W**(k * n) for n in range(4)) for k in range(N)])
print("N=8: X =", X)
print("uguale a fft:", np.allclose(X, np.fft.fft(s)), " sum|X|^2 =", np.sum(abs(X)**2), " N sum|s|^2 =", N * np.sum(s**2))

# ---------- (b) coseno campionato: periodo intero / non intero di campioni
T, N = 0.01, 32
F = 1 / (N * T)                                    # 3.125 Hz; Fp = 100 Hz
n = np.arange(N)
for rapporto in (4.0, 4.5):                        # f0 = rapporto * F
    f0 = rapporto * F
    s = np.cos(2 * np.pi * f0 * n * T)
    X = np.fft.fft(s)
    Sh = X / N                                     # coefficienti S_h = F S(hF)
    S = T * X                                      # S(kF), convenzione del corso
    p = abs(Sh)**2
    vicini = sorted({int(np.floor(rapporto)), int(np.ceil(rapporto)), N - int(np.floor(rapporto)), N - int(np.ceil(rapporto))})
    print("f0 = %.4f Hz (%.1f F)  max|S_h| = %.4f  max|S(kF)| = %.4f" % (f0, rapporto, abs(Sh).max(), abs(S).max()),
          " potenza nei bin", vicini, "= %.1f%%" % (100 * p[vicini].sum() / p.sum()),
          " Parseval:", np.sum(T * s**2), np.sum(F * abs(S)**2))

# ---------- (c) onda quadra campionata: S_r = somma dei coefficienti X_{r+kN} (aliasing)
def X_quadra(h):                                   # onda +-1, pari, +1 per |t| < Tp/4
    h = np.asarray(h, float); out = np.zeros_like(h); nz = h != 0
    out[nz] = 2 * np.sin(np.pi * h[nz] / 2) / (np.pi * h[nz]); return out

for N in (6, 10, 18, 34):                          # N = 2 mod 4: nessun campione cade su un salto
    t = np.arange(N) / N                           # tempo in unita' di Tp
    tt = (t + 0.5) % 1 - 0.5
    s = np.where(np.abs(tt) < 0.25, 1.0, -1.0)
    Sh = np.fft.fft(s) / N
    K = np.arange(-2_000_000, 2_000_001)           # serie dell'aliasing, troncata
    alias = [np.sum(X_quadra(r + K * N)) for r in (1, 3)]
    print("N = %2d  S_1 = %.5f (alias %.5f, X_1 = %.5f)  S_3 = %.5f (alias %.5f, X_3 = %.5f)  potenza = %.4f"
          % (N, Sh[1].real, alias[0], X_quadra([1])[0], Sh[3].real, alias[1], X_quadra([3])[0], np.sum(abs(Sh)**2)))

# ---------- (e) segnale finito: la DFT campiona la trasformata a tempo discreto
T, A0 = 0.5, 2.0
for N in (8, 32):
    s = np.zeros(N); s[[0, 1, N - 1]] = A0         # tre campioni in n = -1, 0, 1
    F = 1 / (N * T); k = np.arange(N)
    print("N =", N, " DFT = S2(kF):", np.allclose(T * np.fft.fft(s), A0 * T * (1 + 2 * np.cos(2 * np.pi * k * F * T))))

Risultati ottenuti eseguendolo: (a) $X=(2,\,1-3i,\,0,\,1+3i)$, $S(kF)=(1,\,0{,}5-1{,}5i,\,0,\,0{,}5+1{,}5i)$, inversa e Parseval vere ($3=3$), $\mathbf U$ unitaria; per $N=8$, $X=(4,\,1-2{,}4142i,\,0,\,1-0{,}4142i,\,0,\,1+0{,}4142i,\,0,\,1+2{,}4142i)$ e $\sum|X|^2=32$. (b) caso A: massimo $|S_h|=0{,}5$, $|S(kF)|=0{,}16$, $100\%$ della potenza nei bin $4$ e $28$; caso B: massimo $0{,}3311$, $|S(kF)|\le0{,}106$, $81{,}2\%$ nei quattro bin centrali; Parseval $0{,}16=0{,}16$ in entrambi. (c) $S_1=0{,}66667,\ 0{,}64721,\ 0{,}63986,\ 0{,}63753$ e $S_3=-0{,}33333,\ -0{,}24721,\ -0{,}22222,\ -0{,}21495$ per $N=6,10,18,34$, uguali alle somme di aliasing; potenza $1$. (e) la DFT coincide con $S_2(kF)$ per $N=8$ e $N=32$.

Errori comuni

• Scambiare la DFT di NumPy con quella del corso: manca il fattore $T$ ($S(kF)=T\,\texttt{fft}$) e, per i coefficienti, il fattore $\frac1N$ ($S_h=\texttt{fft}/N$).
• Dimenticare che $F=\frac1{NT}$ dipende dalla durata della finestra, non dalla velocità: aumentare $N$ a parità di $T$ migliora la risoluzione.
• Credere che i coefficienti della versione campionata di un segnale continuo siano uguali a quelli del continuo: sono la somma $\sum_kX_{r+kN}$ (aliasing).
• Concludere che un coseno a frequenza non multipla di $F$ "ha la DFT sbagliata": la DFT è corretta per il segnale periodico che ripete la finestra; il leakage dipende dalla scelta della finestra.
• Dimenticare che l'ultima metà dei $k$ ($k>N/2$) rappresenta le frequenze negative: $X_{N-k}=\overline{X_k}$ per segnali reali.