Esercizio - filtraggio di un impulso rettangolare con convoluzione e trasformate in Python (laboratorio 5)

Questa pagina non ha ancora la versione ripasso: qui sotto c'è il testo completo.

Testo (laboratorio 5, parti 1-2 e esercizio.m, corso Teoria dei Segnali, UniPD). Un impulso rettangolare $s(t)=A\,\mathrm{rect}\!\left(\frac{t}{2T'}\right)$ entra in un filtro LTI. Si chiede di calcolare l'uscita $y=s*g$ in due modi, e di confrontarli con la formula analitica:

1. con la convoluzione numerica $y(t)=\int s(u)\,g(t-u)\,du$, partendo dalla risposta impulsiva $g(t)$;
2. con le trasformate: $S=\mathcal F[s]$, $Y=G\,S$, $y=\mathcal F^{-1}[Y]$ (nello schema: fft, prodotto, ifft).

Due filtri:

• (a) esponenziale causale $g(t)=\frac1\tau e^{-t/\tau}\,1(t)$, con impulso $A=3$ V su $|t|\le5$ s (laboratorio 5, parti 1-2);
• (b) passa-basso ideale $G(f)=\mathrm{rect}\!\left(\frac{f}{2F'}\right)$ con $A=3$ V, $T'=2$ s, $F'=0{,}5$ Hz (file esercizio.m).

Teoria usata: ConvoluzioneLa convoluzione $xy(t)=\int x(u),y(t-u),du$ combina due segnali ribaltando e traslando il secondo, moltiplicando e integrando. È commutativa, associativa, lineare; l'area del risultato è il prodotto delle aree; l'estensione è la somma delle estensioni (estremo con estremo); l'impulso $\delta$ è l'elemento neutro; la convoluzione con il gradino integra. Per due segnali periodici di uguale periodo si usa la convoluzione ciclica (integrale su un periodo). È l'operazione del filtraggio: l'uscita di un filtro è la convoluzione dell'ingresso con la risposta impulsiva.Convoluzione →, Sistemi lineari tempo-invarianti e risposta impulsivaUna tf lineare e tempo-invariante (LTI, filtro) ha nucleo h(t,u) = g(t-u): l'uscita è la convoluzione y = gx con la risposta impulsiva g (uscita all'impulso ideale nell'origine). Causale se e solo se g è causale; stabile BIBO se e solo se g è assolutamente integrabile (sommabile); reale se e solo se g è reale. Cascata: g = g2g1; parallelo: g1+g2; retroazione: Ge = G/(1+HG) in frequenza.Sistemi lineari tempo-invarianti e risposta impulsiva →, Trasformata di FourierLa trasformata di Fourier $S(f)=\int s(t)e^{-i2\pi ft}dt$ associa a un segnale continuo (anche aperiodico) la sua rappresentazione in frequenza; l'antitrasformata $s(t)=\int S(f)e^{i2\pi ft}df$ lo ricostruisce, perché gli esponenziali $e^{i2\pi ft}$ sono ortogonali su tutto $\mathbb R$ ($\int e^{i2\pi ft}dt=\delta(f)$). Per un segnale reale $S(-f)=S^(f)$. Si calcola per i segnali notevoli (rect $\leftrightarrow$ sinc, $e^{-\alpha t}\mathbf 1(t)\leftrightarrow\frac1{\alpha+i2\pi f}$, gaussiana, $\delta\leftrightarrow1$, $1\leftrightarrow\delta$, gradino) e per i segnali periodici, la cui trasformata è un treno di impulsi di area $S_n$ in $nF$.Trasformata di Fourier →, Proprietà della trasformata di FourierLe regole della trasformata di Fourier trasformano operazioni sui segnali in operazioni sulle trasformate: linearità, ribaltamento, coniugio, traslazione nel tempo ($\times e^{-i2\pi ft_0}$) e in frequenza, convoluzione $\leftrightarrow$ prodotto, cambio di scala $s(at)\to\frac1{|a|}S(f/a)$, derivazione ($\times i2\pi f$), integrazione, regola di simmetria ($S(t)\to s(-f)$). Area $S(0)=\int s$, teorema di Parseval $\int|s|^2=\int|S|^2$. Durata e banda sono inversamente legati e un segnale non può avere durata e banda entrambe limitate; la banda del prodotto è la somma delle bande. Con queste regole si ricavano quasi tutte le trasformate senza integrare.Proprietà della trasformata di Fourier →, 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 →, Risposta in frequenza e filtriLa risposta in frequenza G(f) = F[g] di un filtro LTI dà Y = G·X. Gli esponenziali complessi sono autofunzioni (autovalore G(f)), quindi un ingresso sinusoidale esce sinusoidale con ampiezza moltiplicata per |G(f0)| e fase aumentata di arg G(f0). Il filtro è reale se G è hermitiana, invertibile se G non si annulla. I filtri ideali sono rect in frequenza; non distorsione secondo Heaviside: |G| costante e fase lineare.Risposta in frequenza e filtri →, Segnali notevoli - gradino, rect, tri, sinc ed esponenzialiI segnali di uso più frequente sono la costante, la sinusoide $A_0\cos(2\pi f_0t+\varphi_0)$ e l'esponenziale complesso $Ae^{i2\pi f_0t}$ (periodici, a potenza finita), il gradino $\mathbf 1(t)$ e il segno, e gli impulsi a energia finita: $\operatorname{rect}$ (area $D$), $\operatorname{tri}$, $\operatorname{sinc}$ (area $1$), la gaussiana $e^{-\pi t^2}$ e gli esponenziali smorzati. Per ciascuno si sanno a memoria forma, area ed energia; gli altri segnali si ottengono da questi con traslazioni, scalature, somme e differenze.Segnali notevoli - gradino, rect, tri, sinc ed esponenziali →.

(1) Dal continuo alla griglia: perché compare $dt$

Un computer lavora su campioni, cioè su un dominio discreto $\mathbb Z(dt)$ (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 →): si sceglie il passo $dt$, il numero di campioni $N$ e si ottiene l'asse dei tempi $t_n=n\,dt-T_1$, $n=0,\dots,N-1$, con $dt=2T_1/N$. Un integrale su $\mathbb R$ diventa una somma pesata con $dt$ (integrale di Haar su $\mathbb Z(dt)$: area $=\sum_n dt\,x(t_n)$), ed è proprio questo peso che trasforma la convoluzione di sequenze del computer in una buona approssimazione di quella continua: $$y(t_n)=\int s(u)g(t_n-u)\,du\ \approx\ dt\sum_k s(t_k)\,g(t_n-t_k).$$ np.convolve(s, g) calcola la somma senza il peso: va moltiplicata per $dt$ (in Matlab dt*conv(s,g)). Senza il fattore, l'uscita risulterebbe $1/dt$ volte troppo grande.

Asse dei tempi della convoluzione. Se $s$ comincia in $t_0=-T_1$ e $g$ comincia in $-T_1$, il primo campione di np.convolve (indici $0+0$) corrisponde a $t_0+t_0=-2T_1$; poi si avanza di $dt$. Il risultato ha $2N-1$ campioni, quindi l'asse è t_conv = 2*t[0] + np.arange(2*N-1)*dt (da $-2T_1$ a $2T_1-2\,dt$). Usare l'asse t di $s$ è un errore classico: il picco compare traslato.

(2) Filtro esponenziale: convoluzione numerica

Impulso $s=A$ per $|t|\le T_2=5$ s, filtro $g=\frac1\tau e^{-t/\tau}1(t)$ con $\tau=1{,}5$ s (area $1$, guadagno in continua $1$). Si prende $N=2000$, $T_1=10$ s, quindi $dt=0{,}01$ s.

Formula analitica. Il filtro vale $0$ per $t<0$, quindi in $y(t)=A\int_{t-T_2}^{t+T_2}g(u)\,du$ conta solo la parte con $u\ge0$:

• $t<-T_2$: l'intervallo è tutto a sinistra di $0$, $y=0$;
• $-T_2\le t\le T_2$ (carica): $y=A\int_0^{t+T_2}\frac1\tau e^{-u/\tau}du=A\left(1-e^{-(t+T_2)/\tau}\right)$;
• $t>T_2$ (scarica): $y=A\left(e^{-(t-T_2)/\tau}-e^{-(t+T_2)/\tau}\right)=A\left(e^{T_2/\tau}-e^{-T_2/\tau}\right)e^{-t/\tau}$.

È la risposta di un circuito RC a un impulso di durata $2T_2$ (Risposta in frequenza e filtriLa risposta in frequenza G(f) = F[g] di un filtro LTI dà Y = G·X. Gli esponenziali complessi sono autofunzioni (autovalore G(f)), quindi un ingresso sinusoidale esce sinusoidale con ampiezza moltiplicata per |G(f0)| e fase aumentata di arg G(f0). Il filtro è reale se G è hermitiana, invertibile se G non si annulla. I filtri ideali sono rect in frequenza; non distorsione secondo Heaviside: |G| costante e fase lineare.Risposta in frequenza e filtri →): sale esponenzialmente verso $A$, poi decade con costante di tempo $\tau$.

(3) Filtro esponenziale: con le trasformate

Si passa a $N=4000$, $T_1=100$ s (finestra più lunga, perché il prodotto delle trasformate calcolato con la FFT equivale a una convoluzione circolare: la coda dell'esponenziale non deve "riavvolgersi" sull'inizio della finestra), $dt=0{,}05$ s, $\tau=2$ s. Il passo in frequenza è fissato dalla durata della finestra, $df=\frac1{N\,dt}=\frac1{2T_1}=0{,}005$ Hz, e l'asse è $f_k=k\,df$ per $k=-N/2,\dots,N/2-1$.

Ogni fattore ha un motivo preciso:

• ifftshift prima della fft. fft presuppone che il campione di indice $0$ sia il tempo $t=0$. Il nostro vettore ha $t=0$ a metà (indice $N/2$): ifftshift lo riporta all'inizio (Matlab: s = [s(N/2+1:N), s(1:N/2)]). Senza questa mossa la trasformata acquista una fase spuria $e^{-i2\pi f\,T_1}$.
• Fattore $dt$ nella trasformata. $S(f_k)=\int s(t)e^{-i2\pi f_kt}dt\approx dt\sum_n s_n e^{-i2\pi kn/N}$: la somma è esattamente fft, il peso $dt$ è l'integrale di Haar. Controllo: $S(0)$ deve essere l'area di $s$, $A\cdot2T_2=30$ V s.
• fftshift dopo la fft. Porta $f=0$ al centro, in modo da avere l'asse $f_k$ simmetrico (Matlab: S = [S(N/2+1:N), S(1:N/2)]).
• Fattore $N\,df$ nell'inversa. $y(t_n)=\int Y(f)e^{i2\pi ft_n}df\approx df\sum_kY_ke^{i2\pi kn/N}$, ma ifft divide per $N$: si moltiplica per $N\,df=\frac1{dt}$ (Matlab: y = N*df*ifft(Y)).
• Prodotto $Y=GS$. Il filtraggio diventa un prodotto punto per punto (ConvoluzioneLa convoluzione $x*y(t)=\int x(u),y(t-u),du$ combina due segnali ribaltando e traslando il secondo, moltiplicando e integrando. È commutativa, associativa, lineare; l'area del risultato è il prodotto delle aree; l'estensione è la somma delle estensioni (estremo con estremo); l'impulso $\delta$ è l'elemento neutro; la convoluzione con il gradino integra. Per due segnali periodici di uguale periodo si usa la convoluzione ciclica (integrale su un periodo). È l'operazione del filtraggio: l'uscita di un filtro è la convoluzione dell'ingresso con la risposta impulsiva.Convoluzione →, Proprietà della trasformata di FourierLe regole della trasformata di Fourier trasformano operazioni sui segnali in operazioni sulle trasformate: linearità, ribaltamento, coniugio, traslazione nel tempo ($\times e^{-i2\pi ft_0}$) e in frequenza, convoluzione $\leftrightarrow$ prodotto, cambio di scala $s(at)\to\frac1{|a|}S(f/a)$, derivazione ($\times i2\pi f$), integrazione, regola di simmetria ($S(t)\to s(-f)$). Area $S(0)=\int s$, teorema di Parseval $\int|s|^2=\int|S|^2$. Durata e banda sono inversamente legati e un segnale non può avere durata e banda entrambe limitate; la banda del prodotto è la somma delle bande. Con queste regole si ricavano quasi tutte le trasformate senza integrare.Proprietà della trasformata di Fourier →); con il filtro esponenziale $G(f)=\frac1{1+i2\pi f\tau}$.

import numpy as np
from numpy.fft import fft, ifft, fftshift, ifftshift
import matplotlib.pyplot as plt

A, T2 = 3.0, 5.0                                   # impulso: A volt su |t| <= T2 (durata 10 s)

def y_analitica(t, tau):                           # s * g con g = (1/tau) e^{-t/tau} 1(t)
    y = np.zeros_like(t)
    on = (t >= -T2) & (t <= T2)
    y[on] = A * (1 - np.exp(-(t[on] + T2) / tau))                       # fase di carica
    y[t > T2] = A * (np.exp(T2 / tau) - np.exp(-T2 / tau)) * np.exp(-t[t > T2] / tau)   # scarica
    return y

# ---------- parte 1: convoluzione numerica ----------
N, T1, tau = 2000, 10.0, 1.5
dt = 2 * T1 / N                                    # passo: 0.01 s
t = np.arange(N) * dt - T1                         # asse: -10 ... 9.99 s
s = A * (np.abs(t) <= T2)
g = (1 / tau) * np.exp(-t / tau) * (t >= 0)

y = dt * np.convolve(s, g)                         # MATLAB: dt*conv(s, g); il fattore dt e' il peso dell'integrale
t_conv = 2 * t[0] + np.arange(2 * N - 1) * dt      # asse della convoluzione: parte da t[0]+t[0] = -20 s
print("lunghezze:", len(y), "asse da", t_conv[0], "a", round(t_conv[-1], 2))
print("area di s, g, y:", dt * s.sum(), dt * g.sum(), dt * y.sum(), "(attese 30, 1, 30)")
print("errore massimo rispetto alla formula:", np.abs(y - y_analitica(t_conv, tau)).max())
for tt in (-5.0, 0.0, 5.0, 10.0):
    i = np.argmin(np.abs(t_conv - tt))
    print(f"  y({tt:5.1f}) numerica {y[i]:.4f}  analitica {y_analitica(t_conv[i:i+1], tau)[0]:.4f}")

# ---------- parte 2: stessa cosa con le trasformate (tau = 2 s) ----------
N, T1, tau = 4000, 100.0, 2.0
dt = 2 * T1 / N                                    # 0.05 s
df = 1 / (N * dt)                                  # 0.005 Hz: passo in frequenza
t = np.arange(N) * dt - T1
f = np.arange(-N // 2, N // 2) * df
s = A * (np.abs(t) <= T2)
g = (1 / tau) * np.exp(-t / tau) * (t >= 0)

S = fftshift(dt * fft(ifftshift(s)))               # MATLAB: S = dt*fft([s(N/2+1:N), s(1:N/2)]) poi fftshift
G = fftshift(dt * fft(ifftshift(g)))
Y = S * G                                          # filtraggio = prodotto in frequenza
y2 = fftshift(N * df * ifft(ifftshift(Y)))         # antitrasformata: fattore N*df = 1/dt
print("S(0) =", S[N // 2].real, "(= area di s = 30);  G(0) =", G[N // 2].real, "(= 1)")
print("G rispetto a 1/(1+i2 pi f tau), |f|<5:", np.abs(G - 1 / (1 + 2j * np.pi * f * tau))[np.abs(f) < 5].max())
print("parte immaginaria di y:", np.abs(y2.imag).max())
print("errore massimo rispetto alla formula:", np.abs(y2.real - y_analitica(t, tau)).max())
yc = dt * np.convolve(s, g)                        # la stessa uscita via convoluzione
tc = 2 * t[0] + np.arange(2 * N - 1) * dt
i0 = np.argmin(np.abs(tc - t[0]))                  # indice di yc corrispondente a t[0]
print("differenza FFT - conv:", np.abs(y2.real - yc[i0:i0 + N]).max())

# l'errore e' di ordine dt: dimezzando il passo si dimezza
for Nn in (2000, 4000, 8000):
    dtn = 2 * T1 / Nn; tn = np.arange(Nn) * dtn - T1
    yn = dtn * np.convolve(A * (np.abs(tn) <= T2), (1 / tau) * np.exp(-tn / tau) * (tn >= 0))
    print(f"dt = {dtn:.4f}: errore massimo {np.abs(yn - y_analitica(2 * tn[0] + np.arange(2 * Nn - 1) * dtn, tau)).max():.4f}")

fig, ax = plt.subplots(3, 1)
ax[0].plot(t, s); ax[0].set_xlim(-10, 10)
ax[1].plot(f, np.abs(S)); ax[1].set_xlim(-2, 2)
ax[2].plot(t, y2.real, label="FFT"); ax[2].plot(t, y_analitica(t, tau), "r--", label="formula"); ax[2].set_xlim(-10, 10)
plt.show()

Cosa si vede e valori ottenuti.

• Figura 1: il rettangolo di ampiezza $3$ V su $[-5,5]$ s; $|S(f)|$ è un seno cardinale $2T_2A\,\mathrm{sinc}(2T_2f)=30\,\mathrm{sinc}(10f)$ con picco $30$ V s in $f=0$ e primo zero in $f=0{,}1$ Hz; $|G(f)|$ parte da $1$ e scende come $1/\sqrt{1+(2\pi f\tau)^2}$.
• Figura 2 (uscita): sale da $0$ in $t=-5$ s verso $3$ V con costante di tempo $\tau$, arriva a poco meno di $3$ V in $t=5$ s e poi decade; la curva numerica e quella analitica sono sovrapposte.
• Parte 1 ($\tau=1{,}5$, $dt=0{,}01$): 3999 campioni, asse da $-20$ a $19{,}98$ s; area di $y$ $=30{,}09$ (area di $s$ per area di $g$: $30{,}03\times1{,}002$); errore massimo $0{,}02$. Valori: $y(0)=2{,}9033$ contro $2{,}8930$ analitico; $y(5)=3{,}0062$ contro $2{,}9962$; $y(10)=0{,}1035$ contro $0{,}1069$.
• Parte 2 ($\tau=2$, $dt=0{,}05$): $S(0)=30{,}15$, $G(0)=1{,}0126$; $y$ calcolata con la FFT è reale a meno di $4\times10^{-16}$, coincide con la convoluzione a $10^{-15}$ (le due strade sono lo stesso conto) e dista dalla formula al massimo $0{,}075$.

Perché c'è un errore e come si riduce. L'errore non viene dalla FFT ma dalla discretizzazione di segnali con salti ($s$ ha due salti, $g$ salta in $t=0$): il campione sul salto pesa $dt$ interamente invece di $dt/2$ (regola dei rettangoli e non dei trapezi), per questo $S(0)=30{,}15$ e non $30$ e $G(0)=1+\frac{dt}{2\tau}$. L'errore è proporzionale a $dt$: con $dt=0{,}1;\ 0{,}05;\ 0{,}025$ vale $0{,}150;\ 0{,}075;\ 0{,}0375$. Restano trascurabili la troncatura della coda di $g$ ($e^{-T_1/\tau}$) e l'avvolgimento circolare, grazie a $T_1=100$ s molto più grande di $\tau$.

Nota di lettura sul file del corso: la parte analitica riusa la variabile Tau rimasta dalla parte 2 ($\tau=2$ s), quindi serve a confrontarsi con l'uscita calcolata con le trasformate; qui y_analitica prende $\tau$ come argomento per evitare l'equivoco.

(4) Filtro passa-basso ideale ($A=3$ V, $T'=2$ s, $F'=0{,}5$ Hz)

Risposta impulsiva. Dalla coppia $\mathrm{sinc}(t/D)\leftrightarrow D\,\mathrm{rect}(Df)$ (Trasformata di FourierLa trasformata di Fourier $S(f)=\int s(t)e^{-i2\pi ft}dt$ associa a un segnale continuo (anche aperiodico) la sua rappresentazione in frequenza; l'antitrasformata $s(t)=\int S(f)e^{i2\pi ft}df$ lo ricostruisce, perché gli esponenziali $e^{i2\pi ft}$ sono ortogonali su tutto $\mathbb R$ ($\int e^{i2\pi ft}dt=\delta(f)$). Per un segnale reale $S(-f)=S^*(f)$. Si calcola per i segnali notevoli (rect $\leftrightarrow$ sinc, $e^{-\alpha t}\mathbf 1(t)\leftrightarrow\frac1{\alpha+i2\pi f}$, gaussiana, $\delta\leftrightarrow1$, $1\leftrightarrow\delta$, gradino) e per i segnali periodici, la cui trasformata è un treno di impulsi di area $S_n$ in $nF$.Trasformata di Fourier →) con $D=\frac1{2F'}=1$ s: $g(t)=2F'\,\mathrm{sinc}(2F't)=\mathrm{sinc}(t)$, che vale $1$ in $t=0$, si annulla in tutti gli istanti interi non nulli e non è causale (si estende anche per $t<0$).

Calcolo numerico. Qui $T_1=100$ s e $N=5000$ ($dt=0{,}04$ s, $df=0{,}005$ Hz). Si definisce $G$ direttamente in frequenza, si ottiene $g$ con l'antitrasformata numerica (stesso schema ifftshift / fattore $N\,df$ / fftshift) e poi si filtra in entrambi i modi. Per la convoluzione l'asse è $t_0+t_{g,0}$, cioè $-100-100=-200$ s.

Formula analitica. Si integra il seno cardinale con la funzione seno integrale $\mathrm{Si}(x)=\int_0^x\frac{\sin u}u\,du$: $$y(t)=A\int_{t-T'}^{t+T'}2F'\,\mathrm{sinc}(2F'u)\,du=\frac A\pi\left[\mathrm{Si}\big(2\pi F'(t+T')\big)-\mathrm{Si}\big(2\pi F'(t-T')\big)\right].$$

import numpy as np
from numpy.fft import fft, ifft, fftshift, ifftshift
from scipy.special import sici
import matplotlib.pyplot as plt

A, Tp, Fp = 3.0, 2.0, 0.5                          # s = A rect(t/(2Tp)), G = rect(f/(2Fp)): A = 3 V, Tp = 2 s, Fp = 0.5 Hz
N, T1 = 5000, 100.0
dt = 2 * T1 / N                                    # 0.04 s
df = 1 / (N * dt)                                  # 0.005 Hz
t = np.arange(N) * dt - T1
f = np.arange(-N // 2, N // 2) * df
s = A * (np.abs(t) <= Tp)
G = (np.abs(f) <= Fp).astype(float)                # passa-basso ideale, guadagno 1

# risposta impulsiva: antitrasformata numerica di G (fattore N*df), poi riordino con fftshift
g = fftshift(N * df * ifft(ifftshift(G))).real     # g(t) = 2 Fp sinc(2 Fp t) = sinc(t)
tg = np.arange(-N // 2, N // 2) * dt               # asse dei tempi di g
print("g(0) =", g[N // 2], "(teorico 2*Fp = 1)")

# (a) filtraggio per convoluzione
y_conv = dt * np.convolve(s, g)
t_y = np.arange(2 * N - 1) * dt + t[0] + tg[0]     # asse della convoluzione: t[0] + tg[0] = -200 s
# (b) filtraggio per trasformate
S = fftshift(dt * fft(ifftshift(s)))
y_fft = fftshift(N * df * ifft(ifftshift(S * G))).real

# formula analitica: y(t) = (A/pi) [Si(2 pi Fp (t+Tp)) - Si(2 pi Fp (t-Tp))]
Si = lambda x: sici(x)[0]
y_an = lambda tt: (A / np.pi) * (Si(2 * np.pi * Fp * (tt + Tp)) - Si(2 * np.pi * Fp * (tt - Tp)))
print("errori massimi: conv", np.abs(y_conv - y_an(t_y)).max(), " FFT", np.abs(y_fft - y_an(t)).max())
print("y(0) =", y_fft[N // 2], " analitico 2A/pi*Si(2pi) =", 2 * A / np.pi * Si(2 * np.pi))
tt = np.linspace(-6, 6, 120001); v = y_an(tt)
print("massimo", v.max(), "in t =", tt[v.argmax()], " (sovraelongazione", v.max() / A - 1, ")")
print("valori nei punti -4..4:", [round(float(y_an(float(k))), 4) for k in range(-4, 5)])

# miglioramento: peso 1/2 sui punti di bordo dei rettangoli (regola dei trapezi)
def rect_half(x, a):                               # 1 dentro, 1/2 esattamente sul bordo, 0 fuori
    return (np.abs(x) < a - 1e-9) + 0.5 * (np.abs(np.abs(x) - a) < 1e-9)
s2, G2 = A * rect_half(t, Tp), rect_half(f, Fp)
y2 = fftshift(N * df * ifft(ifftshift(fftshift(dt * fft(ifftshift(s2))) * G2))).real
print("errore con peso 1/2 ai bordi:", np.abs(y2 - y_an(t)).max(), "; aree:", dt * s2.sum(), G2.sum() * df)

fig, ax = plt.subplots(2, 1)
ax[0].plot(t, s); ax[0].plot(t, y_fft, "k"); ax[0].set_xlim(-10, 10)
ax[1].plot(tg, g); ax[1].set_xlim(-10, 10)
plt.show()

Cosa si vede e valori ottenuti.

• $g(t)$ è un seno cardinale di picco $1{,}005$ (teorico $1$: il piccolo eccesso viene da $G$ che ha $201$ campioni a $1$, cioè banda $1{,}005$ Hz invece di $1$ Hz, perché entrambi gli estremi $\pm0{,}5$ cadono sulla griglia).
• L'uscita non è un rettangolo: i fronti sono addolciti (banda finita) e intorno all'impulso compaiono oscillazioni (fenomeno di Gibbs). I valori esatti della formula sono: $y(\pm4)=0{,}0954$, $y(\pm3)=-0{,}2081$, $y(\pm2)=1{,}4249$, $y(\pm1)=3{,}3677$, $y(0)=2{,}7085$.
• Il massimo è $3{,}3677$ V in $t=\pm1$ s: sovraelongazione del 12,3% rispetto ai $3$ V dell'ingresso. In $t=0$ l'uscita scende a $2{,}7085$ (numerico $2{,}7096$): con una banda di $0{,}5$ Hz e un impulso largo $4$ s, la banda non basta per un plateau piatto. Nei bordi $t=\pm2$ si ha $1{,}4249$ e non $A/2=1{,}5$, perché anche il bordo opposto dell'impulso contribuisce.
• Convoluzione e FFT distano entrambe dalla formula al massimo $0{,}061$, e tra loro $0{,}005$ (la coda del seno cardinale decade lentamente, come $1/t$, e i due metodi la troncano o la avvolgono in modo diverso). L'errore principale è sempre dovuto ai bordi: pesando $1/2$ i campioni sui bordi ($s$ ha area $12{,}0$ e $G$ banda $1{,}0$) l'errore scende a $0{,}0007$.

Confronto dei due metodi

Errori tipici.

• Dimenticare $dt$ nella convoluzione o nella trasformata (risultati scalati di $1/dt$), oppure dimenticare $N\,df$ nell'antitrasformata.
• Non applicare ifftshift prima della fft (fase spuria) o non usare fftshift per l'asse delle frequenze.
• Sovrapporre l'uscita della convoluzione all'asse di $s$ invece che a quello che parte da $t_0+t_{g,0}$.

Vedi anche: Esercizio - filtri FIR e IIR e risposta in frequenza in Python (laboratorio 5), 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 e sistemi in Python - campioni, convoluzione e filtriAl calcolatore un segnale continuo è un vettore di campioni con un asse dei tempi. Con numpy si definiscono i segnali, si approssimano area ed energia con somme moltiplicate per il passo, si calcola la convoluzione continua con np.convolve(x, y)*dt (asse = somma degli istanti iniziali), si simulano filtri con equazioni alle differenze e si verificano numericamente linearità e tempo-invarianza.Segnali e sistemi in Python - campioni, convoluzione e filtri →.