Esercizio - banda di un prodotto e convoluzione di sinc

Testo (dispense del corso Teoria dei Segnali, UniPD, Es. 3.4I, 3.4L).

• Es. 3.4I. Provare che il prodotto $s(t)=x(t)\,y(t)$ di due segnali (reali) a banda rispettivamente $B_x$ e $B_y$ è un segnale a banda limitata con $B_s=B_x+B_y$, e quindi, in particolare, che la banda di $x^2(t)$ è $2B_x$.
• Es. 3.4L. Calcolare la convoluzione dei segnali $x(t)=\operatorname{sinc}\left(\frac t{T_1}\right)$ e $y(t)=\operatorname{sinc}\left(\frac t{T_2}\right)$.

Teoria usata: 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 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 →, 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 →, 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 →, 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 →.

Convenzioni: $S(f)=\int s(t)e^{-i2\pi ft}dt$. Un segnale reale ha banda $B$ (limitata) se la sua trasformata è nulla fuori da $[-B,B]$, cioè $E[s]=[-B,B]$ è l'estensione spettrale (il più piccolo intervallo simmetrico fuori dal quale $S(f)=0$).

Es. 3.4I - la banda di un prodotto

Idea. Il prodotto nel tempo è la convoluzione delle trasformate: $s=xy\ \leftrightarrow\ S=X*Y$, $$S(f)=\int_{-\infty}^{+\infty}X(\lambda)\,Y(f-\lambda)\,d\lambda.$$ Il problema diventa: dove è non nulla la convoluzione di due funzioni nulle fuori da $[-B_x,B_x]$ e $[-B_y,B_y]$?

Dimostrazione. Fissiamo $f$ con $|f|>B_x+B_y$. Nell'integrale l'integrando $X(\lambda)Y(f-\lambda)$ è non nullo solo se contemporaneamente $|\lambda|\le B_x$ (altrimenti $X(\lambda)=0$) e $|f-\lambda|\le B_y$ (altrimenti $Y(f-\lambda)=0$). Ma per la disuguaglianza triangolare $$|f|=|\lambda+(f-\lambda)|\le|\lambda|+|f-\lambda|\le B_x+B_y,$$ che contraddice $|f|>B_x+B_y$. Dunque per quei $f$ l'integrando è nullo per ogni $\lambda$ e $S(f)=0$. Quindi $s$ è a banda limitata e $B_s\le B_x+B_y$.

Uguaglianza. Vediamo perché non si può fare meglio. Supponiamo che $B_x$ e $B_y$ siano le bande vere, cioè gli spettri non siano nulli subito dentro i bordi (dentro $(B_x-\varepsilon,B_x)$ e $(B_y-\varepsilon,B_y)$ per ogni $\varepsilon>0$), e che in quei tratti $X$ e $Y$ non cambino segno (come avviene per i segnali di uso pratico: rettangoli, triangoli, esponenziali...). Per $f=B_x+B_y-\eta$ con $\eta>0$ piccolo, l'integrando è non nullo solo per $\lambda\in[B_x-\eta,B_x]$, dove $\lambda$ è a distanza al più $\eta$ dal bordo di $X$ e $f-\lambda\in[B_y-\eta,B_y]$ è a distanza al più $\eta$ dal bordo di $Y$: l'integrando non cambia segno e non è nullo, quindi $S(f)\ne0$. Perciò $S$ è non nulla arbitrariamente vicino a $B_x+B_y$ e $$\boxed{B_s=B_x+B_y.}$$ (Per spettri molto particolari, con segni che cambiano in modo da cancellarsi, potrebbe valere solo $\le$; la dispensa enuncia l'uguaglianza, vera nei casi normali.)

Caso particolare. Per $y=x$ si ottiene $B_{x^2}=2B_x$: quadrare un segnale raddoppia la banda. Si noti che $x^2$ è reale anche quando lo è $x$. L'interpretazione è che l'elevamento al quadrato crea frequenze somma, come $\cos^2\theta=\frac{1+\cos2\theta}2$ ha frequenza doppia.

Interpretazione sul tempo (dualità). Con lo stesso ragionamento, la durata di una convoluzione nel tempo è la somma delle durate; qui è la banda di un prodotto a essere la somma delle bande.

Esempi.

1. Due coseni. $x=\cos2\pi f_1t$ ha banda $f_1$, $y=\cos2\pi f_2t$ ha banda $f_2$ (con $f_1>f_2$). Il prodotto è $\cos(2\pi f_1t)\cos(2\pi f_2t)=\frac12\cos2\pi(f_1+f_2)t+\frac12\cos2\pi(f_1-f_2)t$ (formule di prostaferesi): la frequenza massima è $f_1+f_2=B_x+B_y$ ✓. Per $f_1=3$, $f_2=2$ Hz: componenti a $5$ e $1$ Hz, banda $5$ Hz.
2. Il quadrato di un sinc. $x=\operatorname{sinc}(t)$ ha $X=\operatorname{rect}(f)$, $B_x=\frac12$. $x^2=\operatorname{sinc}^2(t)$ ha trasformata $\operatorname{tri}(f)=\operatorname{rect}(f)*\operatorname{rect}(f)$, nulla fuori da $[-1,1]$: banda $1=2B_x$ ✓ (vedi Esercizio - dimensioni di una trasformata, modulazione e sinc quadro). Il rettangolo di larghezza $1$ convoluto con sé stesso dà un triangolo di base $2$.
3. Rettangolo e triangolo. Se $X=\operatorname{tri}(f)$ ($B_x=1$) e $Y=\operatorname{rect}(f)$ ($B_y=\frac12$), la convoluzione numerica ha supporto $[-1{,}5;\,1{,}5]$: $B_s=1{,}5=B_x+B_y$ ✓.

Conseguenza per il campionamento. Per ricostruire $x$ dai suoi campioni occorre (teorema del campionamento, 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 →) una frequenza di campionamento $F_c\ge2B_x$. Per ricostruire $x^2$ servono $F_c\ge2\cdot2B_x=4B_x$: operare un'elevazione al quadrato richiede di raddoppiare la frequenza di campionamento, altrimenti si ha aliasing.

Es. 3.4L - convoluzione di due sinc

Trasformate. Dalla coppia $\operatorname{sinc}\left(\frac tT\right)\leftrightarrow T\operatorname{rect}(fT)$ (la trasformata di rect con la regola di simmetria e di scala; area $T$): $$X(f)=T_1\operatorname{rect}(fT_1),\qquad Y(f)=T_2\operatorname{rect}(fT_2).$$ Si noti che il calcolo diretto dell'integrale di convoluzione $\int\operatorname{sinc}\left(\frac u{T_1}\right)\operatorname{sinc}\left(\frac{t-u}{T_2}\right)du$ è difficile (primitive non elementari); in frequenza è un prodotto.

Prodotto in frequenza. Per la regola di convoluzione $x*y\leftrightarrow XY$: $$X(f)Y(f)=T_1T_2\operatorname{rect}(fT_1)\operatorname{rect}(fT_2).$$ $\operatorname{rect}(fT)$ vale $1$ per $|f|<\frac1{2T}$. Il prodotto di due rect centrati in $0$ vale $1$ dove valgono $1$ entrambi, cioè per $|f|<\min\left(\frac1{2T_1},\frac1{2T_2}\right)=\frac1{2\max(T_1,T_2)}$, perché il più piccolo tra due inversi è l'inverso del più grande: $$\operatorname{rect}(fT_1)\operatorname{rect}(fT_2)=\operatorname{rect}\big(f\,T_{\max}\big),\qquad T_{\max}=\max(T_1,T_2).$$ (Il rect più stretto "vince": è quello di banda minore, cioè con $T$ maggiore.)

Antitrasformata. Dalla stessa coppia $\operatorname{sinc}\left(\frac tT\right)\leftrightarrow T\operatorname{rect}(fT)$ con $T=T_{\max}$, il segnale che ha trasformata $\operatorname{rect}(fT_{\max})$ è $\frac1{T_{\max}}\operatorname{sinc}\left(\frac t{T_{\max}}\right)$. Quindi $$x*y(t)=T_1T_2\cdot\frac1{T_{\max}}\operatorname{sinc}\left(\frac t{T_{\max}}\right)=\boxed{T_{\min}\,\operatorname{sinc}\left(\frac t{T_{\max}}\right)},\qquad T_{\min}=\min(T_1,T_2),$$ perché $\frac{T_1T_2}{T_{\max}}=T_{\min}$ (il prodotto $T_1T_2=T_{\min}T_{\max}$).

Controlli e interpretazione.

• Area: l'area di $\operatorname{sinc}\left(\frac tT\right)$ è $T$ ($=T\operatorname{rect}(0)$); l'area della convoluzione deve essere il prodotto delle aree: $T_1T_2$. Dal risultato: area $=T_{\min}\cdot T_{\max}=T_1T_2$ ✓.
• Caso $T_1=T_2=T$: $x*x=T\operatorname{sinc}\left(\frac tT\right)=T\,x$. Il sinc è, a meno del fattore $T$, "idempotente" rispetto alla convoluzione.
• Lettura da filtro: $x(t)=\operatorname{sinc}\left(\frac t{T_1}\right)$ è la risposta impulsiva di un filtro passa-basso ideale con guadagno $T_1$ e banda $\frac1{2T_1}$. Se $y$ ha banda minore ($T_2>T_1$, banda $\frac1{2T_2}<\frac1{2T_1}$) passa intatto, moltiplicato per $T_1$: $x*y=T_1\,y$ ✓ (il risultato con $T_{\min}=T_1$ e $T_{\max}=T_2$). Se invece $y$ ha banda maggiore, il filtro la taglia a $\frac1{2T_1}$ e il risultato è proporzionale a $x$: $T_2x$.
• Valore in $0$: $x*y(0)=T_{\min}$. Per $T_1=2$, $T_2=3$: $2$.

Controllo

import numpy as np
from scipy.integrate import quad

# 3.4I: X = tri(f), Y = rect(f); la convoluzione ha supporto [-1.5, 1.5] = Bx + By
df = 0.001
f = np.arange(-4, 4 + df/2, df)
X = np.maximum(0, 1 - abs(f))
Y = (abs(f) < 0.5).astype(float)
S = np.convolve(X, Y, mode="same") * df
print(f[abs(S) > 1e-9].min(), f[abs(S) > 1e-9].max())           # -1.498 1.499

# x = sinc(t): lo spettro di x^2 (FFT su una griglia fitta) è tri(f), nullo fuori da [-1, 1]
dt, N = 0.05, 2**18
t = (np.arange(N) - N//2) * dt
fr = np.fft.fftshift(np.fft.fftfreq(N, dt))
X2 = np.fft.fftshift(np.fft.fft(np.fft.ifftshift(np.sinc(t)**2))) * dt
m = abs(fr) < 3
print(np.max(abs(X2[m].real - np.maximum(0, 1 - abs(fr[m])))))   # 1.5e-05

# 3.4L: integrale di convoluzione calcolato direttamente (T1 = 2, T2 = 3), a confronto con T_min sinc(t/T_max)
T1, T2 = 2.0, 3.0
u = np.arange(-6000, 6000 + 0.005, 0.01)
for t0 in (0, 1, 1.5, 3, 4.5):
    num = np.trapezoid(np.sinc(u/T1) * np.sinc((t0 - u)/T2), u)
    print(t0, round(num, 5), round(min(T1, T2)*np.sinc(t0/max(T1, T2)), 5))

Risultati: il supporto del prodotto è $[-1{,}498;\,1{,}499]$ (atteso $\pm1{,}5$); lo spettro di $\operatorname{sinc}^2(t)$ differisce da $\operatorname{tri}(f)$ per meno di $2\cdot10^{-5}$ ed è nullo (entro $10^{-3}$) per $|f|>1$; il prodotto di due coseni coincide con la somma di prostaferesi a $10^{-14}$. Per le convoluzioni (T1 = 2, T2 = 3) il confronto tra integrale numerico e formula dà $2{,}00000;\ 1{,}65399;\ 1{,}27324;\ 0{,}00000;\ -0{,}42441$ per $t=0;\,1;\,1{,}5;\,3;\,4{,}5$, valori identici. Controllati anche i casi $(T_1,T_2)=(1,2),(1,1),(0{,}5,3)$.