Esercizio - interpolatori in cascata e interpolazione parabolica a tratti

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

Testo (dispense del corso, Esercizi 13.3A e 13.3B).

• 13.3A. Verificare che la cascata costituita da un filtro interpolatore $\mathbb Z(T)\to\mathbb R$ con risposta impulsiva $g_1(t)$ seguito da un filtro ordinario su $\mathbb R$ con risposta impulsiva $g_2(t)$ è equivalente a un filtro interpolatore con risposta impulsiva $$g(t)=g_1*g_2(t),\qquad t\in\mathbb R.$$
• 13.3B. Studiare lo schema di interpolazione composto da un interpolatore con risposta impulsiva $g_1$ seguito da un filtro con risposta impulsiva $g_2$, dove $$g_1(t)=A_1\,\operatorname{tri}\!\left(\frac t{D_1}\right),\qquad g_2(t)=A_2\operatorname{rect}\!\left(\frac t{D_2}\right),\qquad D_1=\frac T2,\quad D_2=T,\quad A_1A_2=2.$$ Verificare che la corrispondente funzione interpolatrice $g_0$ è data dalla (13.22) delle dispense: $$g_0(t)=\begin{cases}1-2\left(\dfrac tT\right)^2&\lvert t\rvert\le\dfrac T2\\[8pt]2\left(1-\dfrac{\lvert t\rvert}T\right)^2&\dfrac T2<\lvert t\rvert\le T\\[8pt]0&\lvert t\rvert>T.\end{cases}$$

Teoria usata: Interpolazione LTI e teorema del campionamentoIl campionatore R→Z(T) ripete lo spettro con periodo Fc = 1/T (Sc(f) = Σ S(f - kFc), senza fattore 1/T). Il filtro interpolatore Z(T)→R ha y(t) = Σ x(nT) g0(t-nT) con g0 = T g e in frequenza Y = G·X. Se S è nulla fuori da (-B,B) e Fc ≥ 2B, con g0(t) = sinc(Fc t) si ricostruisce esattamente s(t) dai campioni. Altrimenti c'è un errore (in banda per l'aliasing, fuori banda per la parte tagliata), ridotto da un prefiltro anti-aliasing.Interpolazione LTI e teorema del campionamento →, 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 →, 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 →, 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 →, 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 →.

Convenzioni delle dispense: un interpolatore $\mathbb Z(T)\to\mathbb R$ con risposta impulsiva $g$ opera così $$\tilde s(t)=\sum_{n}T\,g(t-nT)\,s(nT)=\sum_ns(nT)\,g_0(t-nT),\qquad g_0(t)=T\,g(t)\ \ \text{(funzione interpolatrice)},$$ e in frequenza $\tilde S(f)=G(f)\,S(f)$ con $S$ periodica e $G$ no. Qui $\operatorname{tri}(x)=1-\lvert x\rvert$ per $\lvert x\rvert<1$, $\operatorname{rect}(x)=1$ per $\lvert x\rvert<\frac12$.

(1) Cascata di interpolatore e filtro (Es. 13.3A)

Nel tempo. L'ingresso è $x(nT)$. Il primo blocco (interpolatore con $g_1$) dà $$x_0(t)=\sum_nT\,g_1(t-nT)\,x(nT),\qquad t\in\mathbb R.$$ Il secondo blocco è un filtro ordinario su $\mathbb R$: $y=g_2*x_0$, cioè $y(t)=\int g_2(t-u)\,x_0(u)\,du$. Per la linearità del filtro si può portare la convoluzione dentro la somma (lo scalare $T\,x(nT)$ esce dall'operazione): $$y(t)=\sum_nT\,x(nT)\,\big[g_2*g_1(\cdot-nT)\big](t).$$ Per l'invarianza temporale, convolvere con $g_2$ una funzione ritardata di $nT$ dà la convoluzione ritardata di $nT$: $\left[g_2*g_1(\cdot-nT)\right](t)=(g_1*g_2)(t-nT)$ (la convoluzione è commutativa, 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 →). Quindi $$y(t)=\sum_nT\,(g_1*g_2)(t-nT)\,x(nT),$$ che è proprio la relazione ingresso-uscita di un interpolatore con risposta impulsiva $g=g_1*g_2$. $\blacksquare$

In frequenza (controllo): il primo blocco dà $X_0(f)=G_1(f)X(f)$ con $X$ periodica; il filtro dà $Y(f)=G_2(f)X_0(f)=G_1(f)G_2(f)X(f)$. La risposta in frequenza della cascata è il prodotto $G=G_1G_2$, che è la trasformata di $g_1*g_2$ ✓. In termini di funzione interpolatrice: $g_0=T\,(g_1*g_2)$.

(2) Il caso triangolo + rettangolo (Es. 13.3B)

Ampiezze. Il triangolo $g_1(t)=A_1\operatorname{tri}(t/D_1)$ con $D_1=T/2$ vale $A_1$ in $t=0$ e si annulla per $\lvert t\rvert\ge T/2$ (area $A_1D_1=A_1T/2$). Il rettangolo $g_2(t)=A_2\operatorname{rect}(t/T)$ vale $A_2$ per $\lvert t\rvert<T/2$ (area $A_2T$). Si pone $C=A_1A_2$.

Convoluzione. $g(t)=g_1*g_2(t)=C\displaystyle\int\operatorname{tri}\!\left(\frac\tau{D_1}\right)\operatorname{rect}\!\left(\frac{t-\tau}{T}\right)d\tau$. Il rettangolo è $1$ per $\tau\in\left(t-\frac T2,\,t+\frac T2\right)$: la convoluzione calcola l'area del triangolo che cade dentro questa "finestra" di lunghezza $T$. Il triangolo ha supporto $\left[-\frac T2,\frac T2\right]$, quindi $g$ ha supporto $[-T,T]$ ed è pari (convoluzione di due funzioni pari). Si calcola per $t\ge0$.

• Caso $0\le t\le\frac T2$. La finestra parte da $a=t-\frac T2\in\left[-\frac T2,0\right]$ e comprende la fine del triangolo: l'area è quella totale $D_1$ meno la parte a sinistra di $a$. Su $\left[-D_1,0\right]$ il triangolo vale $1+\frac\tau{D_1}$, quindi $$\int_{-D_1}^{a}\left(1+\frac\tau{D_1}\right)d\tau=\frac{(a+D_1)^2}{2D_1}=\frac{t^2}{2D_1}=\frac{t^2}{T},$$ perché $a+D_1=t-\frac T2+\frac T2=t$ e $2D_1=T$. Perciò $$g(t)=C\left[\frac T2-\frac{t^2}T\right]=C\,\frac T2\left[1-2\left(\frac tT\right)^2\right].$$
• Caso $\frac T2\le t\le T$. Ora $a=t-\frac T2\in\left[0,\frac T2\right]$: la finestra taglia via tutta la parte sinistra del triangolo e parte da $a$. Su $[0,D_1]$ il triangolo vale $1-\frac\tau{D_1}$, e $$\int_a^{D_1}\left(1-\frac\tau{D_1}\right)d\tau=\frac{(D_1-a)^2}{2D_1}=\frac{(T-t)^2}{T}\quad\text{perché }D_1-a=T-t,$$ quindi $$g(t)=C\,\frac{(T-t)^2}{T}=C\,\frac T2\cdot2\left(1-\frac tT\right)^2.$$
• Caso $t>T$: la finestra non tocca più il triangolo, $g=0$.

Funzione interpolatrice. $g_0=T\,g$ vale in $t=0$: $g_0(0)=C\,\dfrac{T^2}2$. Poiché la funzione interpolatrice deve valere $1$ nell'origine (condizione di corretta interpolazione, $g_0(0)=1$) si deve avere $$C=A_1A_2=\frac{2}{T^2},$$ che con il quanto $T=1$ della dispensa dà $A_1A_2=2$, il dato del testo. (Per un quanto qualsiasi la condizione dimensionale è $A_1A_2=2/T^2$: $g$ ha le dimensioni di $1/\text{tempo}$, $g_0$ è adimensionale.) Con questo $C$, i fattori $C\frac{T}{2}\cdot T=1$ e $$g_0(t)=\begin{cases}1-2\left(\dfrac tT\right)^2&\lvert t\rvert\le\dfrac T2\\[8pt]2\left(1-\dfrac{\lvert t\rvert}T\right)^2&\dfrac T2<\lvert t\rvert\le T\\[8pt]0&\lvert t\rvert>T,\end{cases}$$ cioè esattamente la (13.22). $\blacksquare$

Controlli sulla $g_0$ ottenuta.

• Continuità in $\pm T/2$: dal primo ramo $1-2\cdot\frac14=\frac12$; dal secondo $2\left(1-\frac12\right)^2=\frac12$ ✓.
• Continuità della derivata: in $t=\frac T2$ il primo ramo ha pendenza $-\frac{4t}{T^2}=-\frac2T$ e il secondo $-\frac4T\left(1-\frac tT\right)=-\frac2T$ ✓; in $t=T$ il secondo ramo ha pendenza $0$ (raccordo con lo zero) ✓. In $t=0$ la pendenza è $0$ per simmetria. È la proprietà dei "tre archi di parabola ottenuti imponendo la continuità della derivata".
• Corretta interpolazione: $g_0(0)=1$ e $g_0(\pm T)=2(1-1)^2=0$, $g_0(nT)=0$ per $\lvert n\rvert\ge1$ ✓ (Interpolazione LTI e teorema del campionamentoIl campionatore R→Z(T) ripete lo spettro con periodo Fc = 1/T (Sc(f) = Σ S(f - kFc), senza fattore 1/T). Il filtro interpolatore Z(T)→R ha y(t) = Σ x(nT) g0(t-nT) con g0 = T g e in frequenza Y = G·X. Se S è nulla fuori da (-B,B) e Fc ≥ 2B, con g0(t) = sinc(Fc t) si ricostruisce esattamente s(t) dai campioni. Altrimenti c'è un errore (in banda per l'aliasing, fuori banda per la parte tagliata), ridotto da un prefiltro anti-aliasing.Interpolazione LTI e teorema del campionamento →): l'uscita passa per i campioni.
• Somma delle traslate: per $t\in[0,T]$ contribuiscono $g_0(t)$ e $g_0(t-T)$. Per $t\in\left[0,\frac T2\right]$: $\left[1-2\frac{t^2}{T^2}\right]+2\left(1-\frac{T-t}{T}\right)^2=1-2\frac{t^2}{T^2}+2\frac{t^2}{T^2}=1$ ✓ (per $t\in\left[\frac T2,T\right]$ si ha $2\left(1-\frac tT\right)^2+\left[1-2\left(1-\frac tT\right)^2\right]=1$). Quindi $\sum_ng_0(t-nT)=1$ per ogni $t$: l'interpolatore riproduce i segnali costanti.
• Risposta in frequenza: $G_1(f)=A_1D_1\operatorname{sinc}^2(D_1f)$, $G_2(f)=A_2T\operatorname{sinc}(Tf)$, quindi con $C=2/T^2$ $$G(f)=G_1G_2=\frac{2}{T^2}\cdot\frac T2\cdot T\operatorname{sinc}^2\!\left(\frac{Tf}2\right)\operatorname{sinc}(Tf)=\operatorname{sinc}^2\!\left(\frac{Tf}2\right)\operatorname{sinc}(Tf).$$ Vale $G(0)=1$ (area di $g$ uguale a $1$) e il fattore $\operatorname{sinc}(Tf)$ si annulla in $f=kF_c$, $k\ne0$, con $F_c=1/T$: cancella il centro delle repliche dello spettro periodico $S(f)$, e il fattore $\operatorname{sinc}^2\left(\frac{Tf}{2}\right)$ le attenua ulteriormente (zeri doppi in $f=2kF_c$).

(3) Verifica numerica

Si discretizza il tempo con passo $\Delta t=0{,}0005$ (con $T=1$), si prende una sequenza casuale di 9 campioni, si fa passare nel primo interpolatore ($g_1=2\operatorname{tri}(2t)$), poi si convolve con $g_2=\operatorname{rect}(t)$, e si confronta con l'interpolatore diretto con la $g_0$ della (13.22):

import numpy as np
T, dt = 1.0, 0.0005
t = np.arange(-4, 13, dt)
g1 = lambda x: 2 * np.maximum(0, 1 - np.abs(x) / (T / 2))   # A1 = 2, D1 = T/2
k = np.arange(-T / 2, T / 2 + dt / 2, dt)                      # asse dei tempi di g2
g2 = np.where(np.abs(k) < T / 2, 1.0, 0.0)                     # A2 = 1, D2 = T  (A1*A2 = 2)

def g0(x):                                                     # funzione interpolatrice (13.22)
    a = np.abs(x) / T
    return np.where(a <= .5, 1 - 2 * a**2, np.where(a <= 1, 2 * (1 - a)**2, 0.0))

rng = np.random.default_rng(1)
s = rng.normal(size=9)                                         # 9 campioni casuali
x0 = sum(T * g1(t - n * T) * s[n] for n in range(9))            # uscita del primo interpolatore
y = np.convolve(x0, g2, mode="same") * dt                       # poi il filtro g2: convoluzione
y_dir = sum(s[n] * g0(t - n * T) for n in range(9))             # interpolatore diretto con g0
print(np.abs(y - y_dir).max(), np.abs(y_dir).max())            # 0.0007  1.303

La moltiplicazione per dt serve a trasformare la somma discreta nell'integrale di convoluzione. L'errore massimo è $7\cdot10^{-4}$ su un'ampiezza di $1{,}30$ (è l'errore di discretizzazione dei bordi del rettangolo, proporzionale a $\Delta t$). Sempre numericamente, la convoluzione $g_1*g_2$ vale in $t=0,\ 0{,}25,\ 0{,}5,\ 0{,}75,\ 1$: $T g=1;\ 0{,}8747;\ 0{,}4995;\ 0{,}1247;\ 0$, contro la (13.22) $1;\ 0{,}875;\ 0{,}5;\ 0{,}125;\ 0$ ✓. Con $A_1A_2=2$ e un quanto diverso da $1$ (provati $T=0{,}25$ e $T=2$) il confronto fallisce (scarto $0{,}94$ e $3{,}0$), mentre con $A_1A_2=2/T^2$ ($32$ e $0{,}5$) riesce a $10^{-9}$: conferma che il dato del testo $A_1A_2=2$ è scritto per $T=1$.

Errori tipici.

• Dimenticare che $g_0=T\,g$: la risposta impulsiva $g=g_1*g_2$ ha area $1$ e valore $C\,T/2=1/T$ nell'origine, non $1$.
• Calcolare male i casi della convoluzione: la finestra del rettangolo ha larghezza $T$ e centro in $t$; i due rami si scambiano in $\lvert t\rvert=T/2$ perché il triangolo ha semi-base $T/2$.
• Credere che l'interpolatore di (13.22) sia quello "parabolico" della (13.20), $g_0=\left[1-(t/T)^2\right]\operatorname{rect}(t/2T)$: quella ha un solo arco di parabola e in $\lvert t\rvert=T$ arriva a zero con pendenza $\mp2/T$ (spigolo), mentre la (13.22) si raccorda con derivata continua.

Vedi anche: Esercizio - area ed energia di un segnale campionato e, per la materia gemella, Teorema del campionamento, interpolazione e aliasingTeorema di Shannon: un segnale a banda limitata $\omega_M$ si ricostruisce esattamente dai campioni se $T_c<\pi/\omega_M$ (frequenza di campionamento maggiore di quella di Nyquist $2f_{\max}$), con la formula di interpolazione ideale $x(t)=\sum_nx(nT_c)\operatorname{sinc}\left(\frac{t-nT_c}{T_c}\right)$. Sotto Nyquist c'è aliasing: le frequenze alte si confondono con quelle basse e l'informazione è persa.Teorema del campionamento, interpolazione e aliasing →.