Esercizio - trasformata di un sinc discreto e antitrasformata del coseno al cubo

Testo (dispense del corso Teoria dei Segnali, UniPD, Es. 6.3C e 6.3D).

• 6.3C Dimostrare la relazione $$\operatorname{sinc}(nF_0T)\ \xrightarrow{\ \mathcal F\ }\ \frac1{F_0}\operatorname{rep}_{F_p}\operatorname{rect}\left(\frac f{F_0}\right),\qquad F_p=\frac1T,$$ illustrata nel testo per $F_0T=\frac12$. Guida: conviene dimostrare che il secondo membro ha per antitrasformata il primo membro.
• 6.3D Calcolare l'antitrasformata della funzione $S(f)=A_0\cos^3(2\pi fT)$.

Teoria usata: 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 → (definizione, antitrasformata su un periodo, trasformata dell'impulso ideale traslato), 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 →, 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 → (aliasing), 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 → (Eulero), Binomio di NewtonLa formula per sviluppare (a+b)^n con i coefficienti binomiali.Binomio di Newton →.

Promemoria

Per un segnale discreto su $\mathbb Z(T)$ ($F_p=1/T$):

$$S(f)=\sum_nT\,s(nT)e^{-i2\pi fnT},\qquad s(nT)=\int_{f_0}^{f_0+F_p}S(f)\,e^{i2\pi fnT}\,df\ \ \text{(un qualsiasi periodo)}.$$

La ripetizione periodica è $\operatorname{rep}_{F_p}g(f)=\sum_{k=-\infty}^{+\infty}g(f-kF_p)$. L'impulso ideale discreto vale $\frac1T$ in $n=0$ e $0$ altrove; la sua trasformata è $\delta(nT-n_0T)\xrightarrow{\mathcal F}e^{-i2\pi fn_0T}$ (regola (6.13)). Siccome $S$ è periodica, non ci sono due segnali diversi con la stessa $S$: se un'antitrasformata dà $s$, allora $S$ è la trasformata di $s$.

Es. 6.3C: il sinc discreto

Sia $S(f)=\dfrac1{F_0}\operatorname{rep}_{F_p}\operatorname{rect}\left(\dfrac f{F_0}\right)$. Si calcola la sua antitrasformata e si verifica che sia $\operatorname{sinc}(nF_0T)$.

Caso $F_0\le F_p$ (nessuna sovrapposizione). Su un periodo scelto come $[-\frac{F_p}2,\frac{F_p}2]$ le repliche del rect (larghe $F_0$, centrate in $kF_p$) non si toccano e della ripetizione resta solo la replica centrale, $\frac1{F_0}\operatorname{rect}\left(\frac f{F_0}\right)$, che vale $\frac1{F_0}$ per $|f|<\frac{F_0}2$. Allora

$$s(nT)=\int_{-F_p/2}^{F_p/2}S(f)e^{i2\pi fnT}df=\frac1{F_0}\int_{-F_0/2}^{F_0/2}e^{i2\pi fnT}df=\frac1{F_0}\cdot\frac{e^{i\pi F_0nT}-e^{-i\pi F_0nT}}{i2\pi nT}=\frac{\sin(\pi F_0nT)}{\pi F_0nT}.$$

(Ho usato $\int e^{i2\pi fa}df=\frac{e^{i2\pi fa}}{i2\pi a}$ con $a=nT$, e $\frac{e^{i\varphi}-e^{-i\varphi}}{2i}=\sin\varphi$.) Questo è $\operatorname{sinc}(nF_0T)$, con $\operatorname{sinc}(x)=\frac{\sin\pi x}{\pi x}$. Per $n=0$ l'integrale vale $\frac1{F_0}\cdot F_0=1=\operatorname{sinc}(0)$. $\checkmark$

Dimostrazione valida per ogni $F_0$ (si "srotola" la ripetizione). Si scrive $g(f)=\frac1{F_0}\operatorname{rect}\left(\frac f{F_0}\right)$ e $S=\operatorname{rep}_{F_p}g$. Siccome per $n$ intero $e^{i2\pi(u+kF_p)nT}=e^{i2\pi unT}e^{i2\pi kn}=e^{i2\pi unT}$, ogni replica, tradotta di nuovo nell'origine con $u=f-kF_p$, dà lo stesso integrando:

$$\int_0^{F_p}\sum_kg(f-kF_p)\,e^{i2\pi fnT}df=\sum_k\int_{-kF_p}^{(1-k)F_p}g(u)\,e^{i2\pi unT}du=\int_{-\infty}^{+\infty}g(u)\,e^{i2\pi unT}du ,$$

perché gli intervalli $[-kF_p,(1-k)F_p]$ coprono tutto $\mathbb R$ senza sovrapporsi. L'ultimo integrale è il calcolo di sopra: $\operatorname{sinc}(nF_0T)$. È la stessa idea della dimostrazione del Teorema 6.1 (campionare = ripetere lo spettro).

Il caso di figura, $F_0T=\frac12$. Allora $F_0=\frac{F_p}2$: il rect ha larghezza metà del periodo e altezza $\frac1{F_0}=2T$. Spettro a rettangoli alti $2T$ larghi $\frac{F_p}2$, separati da "buchi" larghi $\frac{F_p}2$; il segnale è $\operatorname{sinc}\left(\frac n2\right)$, che vale $1$ in $n=0$, si annulla per $n$ pari non nullo e vale $\pm\frac2{\pi n}$ per $n$ dispari (per esempio $\frac2\pi=0{,}6366$ in $n=1$). Nel grafico $T=1$ ($F_p=1$, $F_0=\frac12$, altezza $2$).

Cosa succede al variare di $F_0$ rispetto a $F_p$.

Quando $F_0>F_p$ la dimostrazione "srotolata" vale ancora, ma lo spettro non è più un rettangolo. Esempio $F_0=\frac32F_p$ (con $T=1$: $F_0=1{,}5$): la replica centrale occupa $|f|<0{,}75$ e quelle vicine $(0{,}25,1{,}75)$ e $(-1{,}75,-0{,}25)$; in $|f|<0{,}25$ passa solo la replica centrale (valore $\frac1{1{,}5}=0{,}667$) e in $0{,}25<|f|<0{,}5$ ne passano due (valore $2\cdot0{,}667=1{,}333$). Interpretazione: $\operatorname{sinc}(nF_0T)$ è il campionamento del segnale continuo $\operatorname{sinc}(F_0t)$, la cui banda ha estensione $F_0$ (da $-\frac{F_0}2$ a $\frac{F_0}2$): per non avere aliasing bisogna campionare con $F_p\ge F_0$ (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 →). Caso estremo: $F_0=2F_p$ dà $\operatorname{sinc}(2n)=0$ per ogni $n\ne0$, cioè lo stesso segnale discreto di $F_0=F_p$: dai campioni non si distinguono.

Es. 6.3D: antitrasformata di $A_0\cos^3(2\pi fT)$

Idea. $S$ è un polinomio trigonometrico in $e^{\pm i2\pi fT}$: se si scrive come somma di esponenziali $e^{-i2\pi fn_0T}$, ciascuno è la trasformata di un impulso ideale in $n_0T$ (regola (6.13)).

1. Sviluppo con Eulero. Con $\varphi=2\pi fT$, $\cos\varphi=\frac{e^{i\varphi}+e^{-i\varphi}}2$ e, per il cubo del binomio (Binomio di NewtonLa formula per sviluppare (a+b)^n con i coefficienti binomiali.Binomio di Newton →) $(a+b)^3=a^3+3a^2b+3ab^2+b^3$ con $a=e^{i\varphi}$, $b=e^{-i\varphi}$ (e $ab=1$):

$$\cos^3\varphi=\frac18\left(e^{i3\varphi}+3e^{i\varphi}+3e^{-i\varphi}+e^{-i3\varphi}\right).$$

(Controllo: per $\varphi=0$ si ha $\frac{1+3+3+1}8=1=\cos^30$.)

2. Dalle esponenziali agli impulsi. Dalla (6.13), $\delta(nT-n_0T)\xrightarrow{\mathcal F}e^{-i2\pi fn_0T}=e^{-in_0\varphi}$. Quindi:

e per linearità

$$s(nT)=\frac{A_0}8\Big[\delta(nT+3T)+3\,\delta(nT+T)+3\,\delta(nT-T)+\delta(nT-3T)\Big].$$

3. Valori dei campioni. L'impulso ideale vale $\frac1T$ nel suo punto di applicazione e $0$ altrove. Quindi il segnale ha quattro campioni non nulli:

$$\boxed{s(\pm T)=\frac{3A_0}{8T},\qquad s(\pm3T)=\frac{A_0}{8T},\qquad s(nT)=0\ \text{altrove (anche in }n=0)}.$$

Il fattore $\frac1T$ viene dal fatto che $\delta(0)=\frac1T$ (area unitaria con larghezza $T$); gli impulsi hanno area (coefficiente) $\frac{A_0}8$ e $\frac{3A_0}8$. Tenere a mente che il coseno ha già prodotto un $\frac12$ per ogni impulso: se si ragiona con $\cos^3\varphi=\frac{3\cos\varphi+\cos3\varphi}4$, ciascuna $\cos k\varphi$ vale $\frac12[\delta\text{ in }kT+\delta\text{ in }-kT]$ e si ritrovano $\frac34\cdot\frac12=\frac38$ e $\frac14\cdot\frac12=\frac18$ (non $\frac34$ e $\frac14$).

4. Controlli.

• Con l'integrale. Per $n=1$: $s(T)=\int_0^{F_p}A_0\cos^3\varphi\,e^{i\varphi}df$; per la parità del coseno resta la parte reale $A_0\int_0^{F_p}\cos^4\varphi\,df=A_0F_p\cdot\overline{\cos^4}=A_0F_p\cdot\frac38=\frac{3A_0}{8T}$ (il valor medio di $\cos^4$ su un periodo è $\frac38$). ✓
• In $n=0$: $\int_0^{F_p}\cos^3\varphi\,df=0$ (per simmetria, $\cos^3$ ha media nulla): $s(0)=0$. ✓
• Parseval: $\sum_nT|s|^2=T\cdot2\left[\left(\frac{3A_0}{8T}\right)^2+\left(\frac{A_0}{8T}\right)^2\right]=\frac{20A_0^2}{64T}=\frac{5A_0^2}{16T}$, e $\int_0^{F_p}A_0^2\cos^6\varphi\,df=A_0^2F_p\cdot\frac5{16}=\frac{5A_0^2}{16T}$ ($\overline{\cos^6}=\frac5{16}$). ✓
• Simmetria: $S$ reale e pari $\Rightarrow$ $s$ reale e pari; durata limitata ($\{-3T,\dots,3T\}$) perché $S$ è un polinomio trigonometrico di grado $3$ nelle esponenziali $e^{\pm i2\pi fT}$ (ogni esponenziale $e^{-i2\pi fn_0T}$ è un campione in $n_0T$).

Controllo

import numpy as np

# ---------------- Es. 6.3C: sinc(n F0 T) <-> (1/F0) rep_Fp rect(f/F0)
T = 1.0
Fp = 1 / T

def S_rep(f, F0, K=6):
    """(1/F0) * somma_k rect((f - k Fp)/F0), con la regola dell'emivalore"""
    f = np.atleast_1d(f).astype(float)
    tot = np.zeros_like(f)
    for k in range(-K, K + 1):
        x = np.abs((f - k * Fp) / F0)
        tot += np.where(x < 0.5, 1.0, np.where(np.isclose(x, 0.5), 0.5, 0.0))
    return tot / F0

# (1) antitrasformata numerica su un periodo: deve dare sinc(n F0 T) per ogni F0
def antitrasformata(F0, n, N=400001):
    f = np.linspace(-Fp / 2, Fp / 2, N)
    S = S_rep(f, F0)
    return np.array([np.trapezoid(S * np.exp(2j * np.pi * f * m * T), f) for m in n])

n = np.arange(-6, 7)
for F0 in (0.5, 1.0, 1.5, 2.7):
    s_num = antitrasformata(F0, n)
    print("F0 =", F0, " antitrasformata == sinc(n F0 T):",
          np.allclose(s_num.real, np.sinc(n * F0 * T), atol=2e-5))

# (2) trasformata come somma della serie: S(f) = somma di T sinc(n F0 T) cos(2 pi f n T)
# (la serie converge lentamente, come 1/n: si usa una finestra triangolare che attenua Gibbs)
def S_serie(f, F0, Nmax=200000):
    m = np.arange(-Nmax, Nmax + 1)
    w = 1 - np.abs(m) / (Nmax + 1)
    return np.array([np.sum(w * T * np.sinc(m * F0 * T) * np.cos(2 * np.pi * ff * m * T)) for ff in f])

for F0, fs in ((0.5, [0.0, 0.1, 0.4, 0.8]), (1.5, [0.1, 0.4, 0.45])):
    fs = np.array(fs)
    print("F0 =", F0, " serie:", np.round(S_serie(fs, F0), 4), " formula:", np.round(S_rep(fs, F0), 4))

# energia: somma di T sinc^2 = 1/F0 solo se F0 <= Fp
for F0 in (0.5, 1.0, 1.5):
    m = np.arange(-400000, 400001)
    print("F0 =", F0, " somma T sinc^2 =", np.sum(T * np.sinc(m * F0 * T)**2), " 1/F0 =", 1 / F0)

# F0 = 2 Fp: stessi campioni di F0 = Fp
m = np.arange(-6, 7)
print("sinc(2n) == sinc(n):", np.allclose(np.sinc(m * 2.0), np.sinc(m * 1.0), atol=1e-12))

# ---------------- Es. 6.3D: antitrasformata di A0 cos^3(2 pi f T)
A0, T = 3.0, 0.2
Fp = 1 / T
f = np.linspace(0, Fp, 200001)
S = A0 * np.cos(2 * np.pi * f * T)**3
n = np.arange(-5, 6)
s = np.array([np.trapezoid(S * np.exp(2j * np.pi * f * m * T), f) for m in n])
for m, v in zip(n, s):
    print(m, np.round(v.real, 6), np.round(v.imag, 8))
print("attesi:  n=+-1:", 3 * A0 / (8 * T), "  n=+-3:", A0 / (8 * T))
print("Parseval:", np.sum(T * np.abs(s)**2), np.trapezoid(S**2, f), 5 / 16 * A0**2 / T)

Risultati (eseguito): l'antitrasformata coincide con $\operatorname{sinc}(nF_0T)$ per $F_0=0{,}5;\ 1;\ 1{,}5;\ 2{,}7$ (anche con sovrapposizione); la serie dà, per $F_0=0{,}5$, $S=2,\,2,\,0,\,2$ in $f=0;\ 0{,}1;\ 0{,}4;\ 0{,}8$ e, per $F_0=1{,}5$, $0{,}6667,\ 1{,}3333,\ 1{,}3333$ in $f=0{,}1;\ 0{,}4;\ 0{,}45$, come previsto; $\sum T\operatorname{sinc}^2=2{,}000,\ 1{,}000,\ 1{,}111$ per $F_0=0{,}5;\ 1;\ 1{,}5$ (cioè $\frac1{F_0}$ solo quando $F_0\le F_p$; per $F_0=1{,}5$ viene $1{,}111$ e non $0{,}667$, effetto dell'aliasing); con $A_0=3$, $T=0{,}2$ i campioni sono $5{,}625=\frac{3A_0}{8T}$ in $n=\pm1$ e $1{,}875=\frac{A_0}{8T}$ in $n=\pm3$, zero altrove; Parseval: $14{,}0625$ in tutti e tre i modi.

Errori comuni

• Scrivere un solo rect invece della ripetizione periodica: la trasformata di un segnale discreto è sempre periodica.
• Dimenticare che per $F_0>F_p$ i rect si sommano: lo spettro non è più un rettangolo e il segnale non coincide con il segnale a banda $F_0$ da cui è stato campionato.
• In 6.3D, scrivere il coefficiente $\frac14$ invece di $\frac18$ (e $\frac34$ invece di $\frac38$): il coseno porta un $\frac12$ per ogni impulso.
• Dimenticare il $\frac1T$: i campioni valgono $\frac{3A_0}{8T}$ e $\frac{A_0}{8T}$, non $\frac{3A_0}8$ e $\frac{A_0}8$; i coefficienti $\frac{3A_0}8$, $\frac{A_0}8$ sono le aree degli impulsi.