Esercizio - laboratorio 4, trasformata di segnali discreti con la FFT

Testo (laboratorio 4 del corso Teoria dei Segnali, UniPD, file lab4.m, parti 1-8). Un impulso rettangolare discreto di ampiezza $A=2$ V ha $N_1+1=4$ campioni non nulli ($n=0,\dots,3$) con $T=1$ ms. (1-2) Rappresentarlo nel tempo e scrivere la sua trasformata analitica, $X(f)=A\,T\,\dfrac{e^{-i2\pi fT(N_1+1)}-1}{e^{-i2\pi fT}-1}$. (3-5) Calcolarla con la FFT: prima in modo sbagliato (fft nuda e asse in indici), poi corretta con $X=T\cdot\text{fft}(x)$ e $f=kF$, $F=F_p/N$, prima con $N=10$ e poi con $N=1000$. (6-8) Un rect discreto di $11$ campioni centrato nell'origine, $N=1000$ campioni da $-N/2$ a $N/2-1$: calcolare la trasformata con la FFT, osservare la fase e correggerla scambiando le due metà del vettore.

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 →, 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 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 →, Serie notevoli - geometrica, telescopica, armonicaLe serie di cui si conosce il carattere e da usare come termine di paragone: geometrica (converge a 1/(1-q) se |q|<1), telescopiche (somma b_1 - lim b_n, come Mengoli), armonica generalizzata (1/n^alpha converge se e solo se alpha>1).Serie notevoli - geometrica, telescopica, armonica → (la somma dei termini di una progressione geometrica), 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 → (traslazione nel tempo).

Cosa si vuole calcolare e perché

Il segnale è discreto su $\mathbb Z(T)$, la sua trasformata (nel corso, 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 →) è

$$X(f)=\sum_nT\,x(nT)\,e^{-i2\pi fnT},$$

funzione di $f$ continua e periodica di periodo $F_p=\frac1T=1000$ Hz. Il calcolatore non può calcolarla per ogni $f$: la calcola su una griglia $f=kF$. La FFT, applicata a $N$ campioni, dà proprio $X(kF)$ per $k=0,\dots,N-1$ con

$$F=\frac{F_p}{N}=\frac1{NT}\qquad\text{(passo in frequenza = inverso della durata } NT\text{ dei campioni osservati)},$$

purché il segnale sia nullo fuori dai $N$ campioni (altrimenti la FFT calcola la trasformata del segnale periodicizzato ogni $NT$, vedi 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 →). Il lavoro del laboratorio è tutto sulle due scalature: il fattore $T$ davanti alla somma e il significato dell'asse delle frequenze; più la fase nel caso di un segnale non partente da $t=0$.

Parti 1-2 - il segnale e la trasformata analitica

import numpy as np
import matplotlib.pyplot as plt

# PARTE 1 - rect discreto: N = 10 campioni, i primi N1 + 1 = 4 valgono A
T = 1e-3          # quanto temporale [s]
N = 10            # campioni totali
N1 = 3            # l'ultimo campione acceso e' n = N1
A = 2

n = np.arange(N)                  # indici 0..9
t = n * T                         # istanti nT: 0, 1 ms, ..., 9 ms
x = A * (n <= N1)                 # MATLAB: A*(t <= N1*T); con l'indice intero si evita il confronto tra decimali

plt.figure()
plt.stem(t, x, basefmt=' ')
plt.grid(True); plt.ylim(-0.5, 2.5); plt.xlabel('nT [s]')
plt.show()

print("x =", x)

# PARTE 2 - trasformata analitica: X(f) = A T (e^(-i 2 pi f T (N1+1)) - 1) / (e^(-i 2 pi f T) - 1)
Fp = 1 / T                        # periodo della trasformata: 1000 Hz
df = 1
fa = np.arange(0, int(Fp) + 1, df)         # 0, 1, ..., 1000 Hz (1001 punti)

with np.errstate(divide='ignore', invalid='ignore'):
    Xa = A * T * (np.exp(-1j * 2 * np.pi * fa * T * (N1 + 1)) - 1) / (np.exp(-1j * 2 * np.pi * fa * T) - 1)
# in f = 0 (e Fp) il rapporto e' 0/0: si sostituisce con il limite A T (N1+1)
Xa[fa % Fp == 0] = A * T * (N1 + 1)

plt.figure()
plt.plot(fa, np.real(Xa), 'r', label='Re X')
plt.plot(fa, np.imag(Xa), 'm', label='Im X')
plt.xlabel('f [Hz]'); plt.ylabel('X(f)'); plt.legend(); plt.grid(True)
plt.show()

print(f"X(0) = {Xa[0].real:.4f}  (A T (N1+1) = {A * T * (N1 + 1):.4f},  area del segnale)")
print(f"|X| in 250, 500, 750 Hz: {abs(Xa[250]):.1e}, {abs(Xa[500]):.1e}, {abs(Xa[750]):.1e}  (zeri: multipli di Fp/(N1+1) = 250 Hz)")
print(f"X(100 Hz) = {Xa[100]:.5f},  X(900 Hz) = {Xa[900]:.5f}  (coniugati: il segnale e' reale)")

Risultato:

x = [2 2 2 2 0 0 0 0 0 0]
X(0) = 0.0080  (A T (N1+1) = 0.0080,  area del segnale)
|X| in 250, 500, 750 Hz: 3.5e-19, 4.9e-19, 1.0e-18  (zeri: multipli di Fp/(N1+1) = 250 Hz)
X(100 Hz) = 0.00362-0.00498j,  X(900 Hz) = 0.00362+0.00498j  (coniugati: il segnale e' reale)

Il segnale. n <= N1 è vero per $n=0,1,2,3$ (quattro campioni a $A=2$) e falso per $n=4,\dots,9$. Il file MATLAB scrive t <= N1*T, che vale lo stesso ma confronta numeri decimali; con gli indici interi si evita il rischio di un arrotondamento sbagliato. La durata dell'impulso è $4T=4$ ms.

La trasformata. Per l'impulso $x(nT)=A$ per $n=0,\dots,N_1$ la somma è una progressione geometrica di ragione $r=e^{-i2\pi fT}$ con $N_1+1$ termini (Serie notevoli - geometrica, telescopica, armonicaLe serie di cui si conosce il carattere e da usare come termine di paragone: geometrica (converge a 1/(1-q) se |q|<1), telescopiche (somma b_1 - lim b_n, come Mengoli), armonica generalizzata (1/n^alpha converge se e solo se alpha>1).Serie notevoli - geometrica, telescopica, armonica →):

$$X(f)=A\,T\sum_{n=0}^{N_1}r^n=A\,T\,\frac{1-r^{N_1+1}}{1-r}=A\,T\,\frac{e^{-i2\pi fT(N_1+1)}-1}{e^{-i2\pi fT}-1},$$

la formula del laboratorio (si moltiplicano numeratore e denominatore per $-1$). È periodica di periodo $F_p=1000$ Hz e, essendo il segnale reale, $X(F_p-f)=X^*(f)$ (i valori a $100$ Hz e a $900$ Hz sono complessi coniugati, come stampato). Valori utili:

• $X(0)=A\,T\,(N_1+1)=0{,}008$: è l'area del segnale ($\sum_nT\,x(nT)=1\text{ ms}\cdot4\cdot2$). Nel codice è il limite del rapporto $0/0$ (che numpy darebbe nan): lo si sostituisce a mano in $f=0$ e $f=F_p$;
• zeri nei multipli di $\frac{F_p}{N_1+1}=250$ Hz (si annulla il numeratore senza che si annulli il denominatore): $|X|\sim10^{-19}$ a $250$, $500$, $750$ Hz.

Parte 3 - la prima versione, sbagliata

# PARTE 3 - prima versione: fft nuda, tracciata in funzione dell'indice (SBAGLIATA)
X = np.fft.fft(x)

plt.figure()
plt.plot(np.real(X), label='Re X')
plt.plot(np.imag(X), label='Im X')
plt.grid(True); plt.xlabel('f [Hz]'); plt.ylabel('X(f)'); plt.legend()
plt.show()

print("fft(x)[0]  =", X[0], "  contro  X(0) =", Xa[0].real)
print("fft(x)[1]  =", np.round(X[1], 4), " (indice 1: a quale frequenza?)")

Risultato:

fft(x)[0]  = (8+0j)   contro  X(0) = 0.008
fft(x)[1]  = (3.618-4.9798j)  (indice 1: a quale frequenza?)

Due errori in uno. (a) Scala delle ampiezze. np.fft.fft(x) è la somma $\sum_nx_n\,e^{-i2\pi kn/N}$, senza il passo $T$: $X[0]=\sum x=8$ contro i $0{,}008$ del grafico analitico, mille volte di più. Nella definizione del corso ogni campione pesa $T$ (stesso discorso di area ed energia). (b) Asse delle frequenze. fft restituisce $N$ numeri indicati dall'indice $k=0,\dots,N-1$ e il grafico mostra le ascisse $0,\dots,9$; ma il valore di indice $k$ corrisponde alla frequenza $kF$ ($F=100$ Hz): per esempio l'indice $1$ è $100$ Hz, non "1". Il risultato è una figura con ordinate e ascisse sbagliate.

Parte 4 - correzione 1: fattore $T$ e asse in Hz

# PARTE 4 - correzione 1: fattore T e asse delle frequenze f = k F con F = Fp / N
F = Fp / N                        # passo in frequenza: 100 Hz
f = np.arange(N) * F              # 0, 100, ..., 900 Hz
X = T * np.fft.fft(x)

plt.figure()
plt.plot(f, np.real(X), 'r', label='Re X')
plt.plot(f, np.imag(X), 'm', label='Im X')
plt.grid(True); plt.xlabel('f [Hz]'); plt.ylabel('X(f)'); plt.legend()
plt.show()

print("  k   f [Hz]   T*fft(x)[k]            formula analitica        differenza")
for k in range(N):
    fk = int(f[k])
    print(f"{k:>3}  {fk:>6}   {X[k].real:+.5f}{X[k].imag:+.5f}i   {Xa[fk].real:+.5f}{Xa[fk].imag:+.5f}i   {abs(X[k] - Xa[fk]):.1e}")

Risultato:

  k   f [Hz]   T*fft(x)[k]            formula analitica        differenza
  0       0   +0.00800+0.00000i   +0.00800+0.00000i   0.0e+00
  1     100   +0.00362-0.00498i   +0.00362-0.00498i   4.3e-19
  2     200   -0.00062-0.00190i   -0.00062-0.00190i   8.7e-19
  3     300   +0.00138+0.00045i   +0.00138+0.00045i   3.9e-19
  4     400   +0.00162-0.00118i   +0.00162-0.00118i   2.2e-19
  5     500   +0.00000-0.00000i   -0.00000-0.00000i   2.7e-19
  6     600   +0.00162+0.00118i   +0.00162+0.00118i   6.9e-19
  7     700   +0.00138-0.00045i   +0.00138-0.00045i   8.7e-19
  8     800   -0.00062+0.00190i   -0.00062+0.00190i   1.3e-18
  9     900   +0.00362+0.00498i   +0.00362+0.00498i   3.2e-18

Con $X=T\cdot\text{fft}(x)$ e $f=kF$ (in Python np.arange(N) * F) i 10 punti coincidono con la formula analitica fino a $10^{-18}$: per $k=0,\dots,9$ si leggono $X(kF)$ esatti. La fft non approssima niente: dà esattamente la trasformata $X(f)$ nei punti $f=kF$, perché il segnale ha solo 4 campioni non nulli, tutti dentro i $10$ osservati. Il grafico, però, ha solo 10 punti (uno ogni $100$ Hz): non si vede la forma di $X(f)$, e per esempio lo zero a $250$ Hz resta nascosto tra $200$ e $300$. Si noti anche la simmetria $X(F_p-f)=X^*(f)$: $k=9$ è il coniugato di $k=1$, $k=8$ di $k=2$ e così via, e $k=5$ ($500$ Hz $=F_p/2$) vale $0$.

Parte 5 - correzione 2: più campioni (zeri aggiunti) per più punti in frequenza

# PARTE 5 - correzione 2: stessi 4 campioni, ma N = 1000 (si aggiungono zeri)
T = 1e-3
N = 1000
N1 = 3
A = 2

n = np.arange(N)
t = n * T                          # 0 ... 0.999 s
x = A * (n <= N1)

plt.figure()
plt.stem(t, x, basefmt=' ')
plt.grid(True); plt.ylim(-0.5, 2.5); plt.xlim(0, 0.05); plt.xlabel('nT [s]')
plt.show()

F = Fp / N                         # 1 Hz
f = np.arange(N) * F
X = T * np.fft.fft(x)

plt.figure()
plt.plot(f, np.real(X), 'r', label='Re X')
plt.plot(f, np.imag(X), 'm', label='Im X')
plt.grid(True); plt.xlabel('f [Hz]'); plt.ylabel('X(f)'); plt.legend()
plt.show()

print(f"N = {N}: F = {F} Hz, {N} righe da 0 a {f[-1]:.0f} Hz")
print(f"differenza massima con la formula analitica (stesse frequenze): {np.max(np.abs(X - Xa[:-1])):.1e}")

Risultato:

N = 1000: F = 1.0 Hz, 1000 righe da 0 a 999 Hz
differenza massima con la formula analitica (stesse frequenze): 3.8e-17

Il segnale ha ancora 4 campioni non nulli, ma ora il vettore ne ha $N=1000$: si sono aggiunti $996$ zeri (zero padding). Il passo in frequenza diventa $F=\frac{F_p}N=1$ Hz, cioè la durata osservata è $NT=1$ s. La trasformata del segnale è la stessa (gli zeri non cambiano la somma $\sum T\,x\,e^{\dots}$), ma ora è calcolata in $1000$ punti invece di $10$: sono esattamente i punti della formula analitica (differenza $4\cdot10^{-17}$) e il grafico appare come una curva continua. Morale: il passo in frequenza è l'inverso della durata della finestra di osservazione; per vedere meglio $X(f)$ si allunga la finestra, non si cambia $T$. L'asse arriva sempre a $F_p-F$ ($999$ Hz): la metà alta dell'asse è la parte a frequenze negative $f-F_p$.

Parti 6-7 - rect centrato nell'origine: modulo giusto, fase strana

# PARTE 6 - rect discreto centrato nell'origine: 11 campioni a A, N = 1000 campioni da -N/2 a N/2 - 1
T = 1e-3
N = 1000
N1 = 11                            # numero di campioni accesi
A = 2

n = np.arange(-N // 2, N // 2)     # -500, ..., 499
t = n * T
T1 = N1 * T                        # durata dell'impulso: 11 ms
x = A * (np.abs(n) < N1 / 2)       # 1 per n = -5, ..., 5 (|nT| < T1/2)

plt.figure()
plt.stem(t, x, basefmt=' ')
plt.xlabel('nT [s]'); plt.axis([-3 * T1, 3 * T1, -0.1 * A, 1.1 * A]); plt.grid(True)
plt.show()

print(f"campioni accesi: {int((x > 0).sum())}  da n = {n[x > 0][0]} a n = {n[x > 0][-1]};  primo campione dell'array: t = {t[0]} s")

campioni accesi: 11  da n = -5 a n = 5;  primo campione dell'array: t = -0.5 s

Un rect discreto simmetrico rispetto all'origine: $x(nT)=A$ per $|n|\le5$ (11 campioni). Nel file MATLAB N1 = 11 è il numero di campioni accesi (durata $T_1=N_1T=11$ ms), mentre nelle parti 1-5 N1 = 3 era l'indice dell'ultimo (quindi i campioni erano $N_1+1$): attenzione al diverso significato. Il vettore va da $n=-500$ a $n=499$, cioè da $t=-0{,}5$ s a $0{,}499$ s. Teoricamente (somma simmetrica di coseni, formula del seno)

$$X(f)=A\,T\sum_{n=-5}^{5}e^{-i2\pi fnT}=A\,T\,\frac{\sin(\pi fTN_1)}{\sin(\pi fT)},$$

una funzione reale e pari (segnale reale e pari), con $X(0)=A\,T\,N_1=0{,}022$ e primo zero in $f=\frac{F_p}{N_1}=90{,}9$ Hz.

# PARTE 7 - trasformata con fft: modulo giusto, fase strana
Fp = 1 / T
F = Fp / N
f = np.arange(N) * F
X = T * np.fft.fft(x)

# trasformata analitica del rect centrato: reale, A T sin(pi f T N1) / sin(pi f T)
with np.errstate(divide='ignore', invalid='ignore'):
    X_vera = A * T * np.sin(np.pi * f * T * N1) / np.sin(np.pi * f * T)
X_vera[0] = A * T * N1

plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.plot(f, np.abs(X), color='r'); plt.xlabel('f [Hz]'); plt.ylabel('|X(f)| '); plt.grid(True)
plt.subplot(1, 2, 2)
plt.plot(f, np.angle(X), color='b'); plt.xlabel('f [Hz]'); plt.ylabel('arg(X(f))'); plt.ylim(-1.5 * np.pi, 1.5 * np.pi); plt.grid(True)
plt.show()

print("modulo:  max | |X| - |X_vera| | =", f"{np.max(np.abs(np.abs(X) - np.abs(X_vera))):.1e}")
print("X[0..5]       =", np.round(X[:6].real, 5), " (parte immaginaria massima:", f"{np.max(np.abs(X.imag)):.1e})")
print("X_vera[0..5]  =", np.round(X_vera[:6], 5))
print("X / X_vera per k = 0..7:", np.round((X / X_vera)[:8].real, 6))
print("fase di X per k = 0..7 :", np.round(np.angle(X[:8]), 4))

Risultato:

modulo:  max | |X| - |X_vera| | = 1.1e-15
X[0..5]       = [ 0.022   -0.022    0.02198 -0.02196  0.02193 -0.02189]  (parte immaginaria massima: 2.2e-18)
X_vera[0..5]  = [0.022   0.022   0.02198 0.02196 0.02193 0.02189]
X / X_vera per k = 0..7: [ 1. -1.  1. -1.  1. -1.  1. -1.]
fase di X per k = 0..7 : [ 0.      3.1416 -0.     -3.1416  0.      3.1416  0.      3.1416]

Il modulo è giusto (differenza con la formula $10^{-15}$): il modulo non risente di traslazioni nel tempo. La fase no. Il vero $X$ è reale e positivo vicino a $f=0$, quindi dovrebbe avere fase $0$; la fft invece dà, per $k=0,1,2,\dots$, valori che alternano segno: $X[k]=(-1)^kX_{vera}(kF)$, come dice il rapporto stampato, e np.angle restituisce $0,\pi,0,-\pi,\dots$ (alternanza fra $0$ e $\pm\pi$).

Il perché. La fft suppone che il primo elemento del vettore sia l'istante $0$: essa calcola $\sum_{m=0}^{N-1}x[m]\,e^{-i2\pi km/N}$ con $m$ che parte da $0$. Ma qui il primo elemento (m = 0) è l'istante $t=-0{,}5$ s, e l'istante $0$ è m = N/2 = 500: il segnale visto dalla fft è quello vero ritardato di $N/2$ campioni, cioè di $\frac{NT}2=0{,}5$ s. Per la regola della traslazione (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 →) un ritardo $t_0$ moltiplica la trasformata per $e^{-i2\pi ft_0}$, una fase a rampa (retta in $f$) di pendenza $-2\pi t_0$. Sui punti $f=kF=\frac{k}{NT}$ con $t_0=\frac{NT}2$ il fattore è $e^{-i2\pi\frac k{NT}\frac{NT}2}=e^{-i\pi k}=(-1)^k$: la rampa è così ripida che, riportata modulo $2\pi$ sui punti della griglia, salta tra $0$ e $\pi$ (si vede un'alternanza invece di una retta). Il blocco seguente verifica la regola per ritardi diversi: per $m_0$ campioni di ritardo il fattore è $e^{-i2\pi km_0/N}$ (una rampa vera, di pendenza dipendente da $m_0$); per $m_0=N/2$ diventa $(-1)^k$.

# Il perche': ritardare un segnale di m0 campioni moltiplica la fft per e^(-i 2 pi k m0 / N)
k = np.arange(N)
x_centrato_in_zero = np.fft.ifftshift(x)                # lo 0 in prima posizione
for m0 in [1, 7, N // 2]:
    x_rit = np.roll(x_centrato_in_zero, m0)             # segnale ritardato di m0 campioni
    rapporto = np.fft.fft(x_rit) / np.fft.fft(x_centrato_in_zero)
    atteso = np.exp(-1j * 2 * np.pi * k * m0 / N)
    ok = np.allclose(rapporto[np.abs(np.fft.fft(x_centrato_in_zero)) > 1e-6], atteso[np.abs(np.fft.fft(x_centrato_in_zero)) > 1e-6])
    print(f"ritardo di {m0:>3} campioni: rapporto = e^(-i 2 pi k {m0}/N)?  {ok}")
print("e per m0 = N/2 = 500 il fattore e^(-i pi k) = (-1)^k:", np.allclose(atteso[:6], (-1.0) ** np.arange(6)))

ritardo di   1 campioni: rapporto = e^(-i 2 pi k 1/N)?  True
ritardo di   7 campioni: rapporto = e^(-i 2 pi k 7/N)?  True
ritardo di 500 campioni: rapporto = e^(-i 2 pi k 500/N)?  True
e per m0 = N/2 = 500 il fattore e^(-i pi k) = (-1)^k: True

Parte 8 - correzione: portare l'istante 0 al primo posto

# PARTE 8 - correzione: portare il campione n = 0 al primo posto, scambiando le due meta' del vettore
Fp = 1 / T
F = Fp / N
f = np.arange(N) * F
x1 = np.concatenate((x[N // 2:], x[:N // 2]))      # MATLAB [x(N/2+1:N), x(1:N/2)]
print("equivale a ifftshift (e a fftshift, perche' N e' pari):", np.array_equal(x1, np.fft.ifftshift(x)))
X1 = T * np.fft.fft(x1)

plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.plot(f, np.abs(X1), color='r'); plt.xlabel('f [Hz]'); plt.ylabel('|X(f)|'); plt.grid(True)
plt.subplot(1, 2, 2)
plt.plot(f, np.imag(X1), color='b')                 # e' quello che traccia il file MATLAB (ma l'etichetta dice arg)
plt.xlabel('f [Hz]'); plt.ylabel('arg(X(f))'); plt.grid(True); plt.ylim(-1.5 * np.pi, 1.5 * np.pi)
plt.show()

print(f"parte immaginaria massima di X1: {np.max(np.abs(X1.imag)):.1e}  (il segnale e' pari: la trasformata e' reale)")
print(f"differenza massima con la formula analitica: {np.max(np.abs(X1 - X_vera)):.1e}")
print("fase corretta (np.angle(X1)) per k = 0..7:", np.round(np.angle(X1[:8]), 4), "...  e' 0 dove X1 > 0 e pi dove X1 < 0")
print("X1[0..7] =", np.round(X1[:8].real, 5))
print(f"zeri della trasformata: ogni Fp/N1 = {Fp / N1:.2f} Hz (prima di tutti: {Fp / N1:.2f})")
print("X1 vicino al primo zero: k = 90, 91, 92 ->", np.round(X1[[90, 91, 92]].real, 5))

Risultato:

equivale a ifftshift (e a fftshift, perche' N e' pari): True
parte immaginaria massima di X1: 2.2e-18  (il segnale e' pari: la trasformata e' reale)
differenza massima con la formula analitica: 1.1e-15
fase corretta (np.angle(X1)) per k = 0..7: [ 0. -0. -0.  0.  0. -0.  0. -0.] ...  e' 0 dove X1 > 0 e pi dove X1 < 0
X1[0..7] = [0.022   0.022   0.02198 0.02196 0.02193 0.02189 0.02184 0.02179]
zeri della trasformata: ogni Fp/N1 = 90.91 Hz (prima di tutti: 90.91)
X1 vicino al primo zero: k = 90, 91, 92 -> [ 2.3e-04 -2.0e-05 -2.6e-04]

Cosa si fa. Per il vettore della fft il tempo va da $0$ a $(N-1)T$; i tempi negativi $-\frac N2T,\dots,-T$ sono "equivalenti modulo $NT$" a $\frac N2T,\dots,(N-1)T$ (stesso discorso della periodicità del segnale per la DFT). Allora si ricostruisce il vettore mettendo prima i campioni da $t=0$ in poi (le posizioni da $N/2$ a $N-1$ di x) e dopo i campioni dei tempi negativi (le posizioni da $0$ a $N/2-1$): x1 = np.concatenate((x[N//2:], x[:N//2])), in MATLAB [x(N/2+1:N), x(1:N/2)]. È lo scambio delle due metà, cioè np.fft.ifftshift(x) (per $N$ pari coincide anche con fftshift).

Cosa cambia. Ora l'elemento 0 è l'istante $0$, il ritardo è nullo e $X_1=T\cdot\text{fft}(x_1)$ è reale (parte immaginaria $\sim10^{-18}$), uguale alla formula teorica entro $10^{-15}$ e senza l'alternanza di segno. La sua fase è $0$ dove $X>0$ ($k\le90$, cioè sotto il primo zero a $90{,}9$ Hz) e $\pi$ dove $X<0$ (il segno cambia a ogni zero: il rect ha una trasformata che oscilla tra lobi positivi e negativi, e questo dà $0$ o $\pm\pi$). Per un segnale pari e reale la fase è sempre $0$ o $\pi$.

Errore del file MATLAB (parte 8). Il secondo grafico ha ylabel('arg(X(f))') e ylim([-1.5*pi 1.5*pi]), ma traccia imag(X1): si vede una riga a zero (la parte immaginaria, nulla perché la trasformata è reale), non la fase. Per vedere la fase bisogna tracciare angle(X1) (in Python np.angle(X1)), come nella parte 7. Nel codice sopra si è lasciato np.imag(X1) proprio per riprodurre il grafico del laboratorio e farne notare la discordanza con l'etichetta.

Controllo

• Il fattore $T$: $X(0)=T\sum x=0{,}008$ (parti 1-5) e $0{,}022$ (parte 6-8), uguali alle aree dei segnali.
• Griglia: $N=10\Rightarrow F=100$ Hz (10 punti esatti), $N=1000\Rightarrow F=1$ Hz (1000 punti esatti).
• Modulo identico prima e dopo lo scambio delle metà; fase corretta $0/\pi$ solo dopo lo scambio.