Esercizio - filtri FIR e IIR e risposta in frequenza in Python (laboratorio 5)

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

Testo (laboratorio 5, parti 3-6, corso Teoria dei Segnali, UniPD). Filtri a tempo discreto implementati con l'equazione alle differenze ($\texttt{lfilter}$ di scipy.signal, in Matlab filter):

• Parte 3. Segnale $x(nT)=\sin(2\pi f_1nT)+0{,}5\sin(2\pi f_2nT)$ con $T=0{,}5$ s, $N=200$ campioni, $f_1=0{,}1$ Hz, $f_2=0{,}01$ Hz. Filtro FIR a media mobile su $L=20$ campioni: $y(nT)=\frac1L\sum_{k=0}^{L-1}x((n-k)T)$.
• Parte 4. Risposta in frequenza del FIR con la FFT e lo zero-padding ($N_{fft}=512$).
• Condizioni iniziali. Uscita con e senza il passato dell'ingresso (filtic).
• Parte 5. Filtro IIR di Butterworth passa-alto, ordine $4$, frequenza di taglio $f_c=0{,}2$ Hz, passo $T_s=0{,}02$ s, ingresso $\sin(2\pi f_1nT_s)+\sin(2\pi f_2nT_s)$ con $f_1=0{,}02$ Hz, $f_2=1{,}6$ Hz, $N=200$.
• Parte 6. Risposta in frequenza dell'IIR calcolata dalla sua risposta impulsiva.

Teoria usata: 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 →, 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 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 zeta e sistemi a tempo discretoUn sistema LTI discreto è descritto da un'equazione alle differenze Σ a_i y(n-i) = Σ b_i x(n-i) con condizioni iniziali. L'uscita è evoluzione libera (dalle condizioni iniziali) più risposta forzata gx. Si calcola in tre modi: soluzione dell'equazione, risposta impulsiva con segnale fittizio, trasformata zeta. FIR: memoria finita, sempre stabile; IIR: stabile se tutti i poli hanno modulo minore di 1. In frequenza G(f) = Σ b_i e^{-i2πfi} / Σ a_i e^{-i2πfi}.Trasformata zeta e sistemi a tempo discreto →, 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 →, 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 →.

(1) Convenzione del corso: coefficienti e risposta impulsiva

lfilter(b, a, x) implementa $a_0y[n]=\sum_kb_kx[n-k]-\sum_{k\ge1}a_ky[n-k]$ con $a_0=1$. Nel corso la convoluzione su $\mathbb Z(T)$ è una somma pesata con $T$ (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 →): $y(nT)=\sum_kT\,g(kT)\,x(nT-kT)$. Per un FIR i due modi coincidono se $$g(kT)=\frac{b_k}{T},\qquad H(f)=\sum_kT\,g(kT)\,e^{-i2\pi fkT}=\sum_kb_k\,e^{-i2\pi fkT},\qquad H \text{ periodica di periodo }F_p=\frac1T.$$ Quindi la risposta impulsiva "del corso" è $b/T$ (l'impulso ideale su $\mathbb Z(T)$ vale $\frac1T$ in $0$: dandolo in ingresso a lfilter si legge $g$), e ogni volta che si trasforma $g$ con la FFT serve il fattore $T$: $H=T\cdot\mathtt{fft}(g)$. Il guadagno in continua $H(0)=\sum_kb_k$ è così $1$ per la media mobile.

(2) Media mobile (parte 3)

Risposta in frequenza. Somma geometrica: $$H(f)=\frac1L\sum_{k=0}^{L-1}e^{-i2\pi fkT}=e^{-i\pi f(L-1)T}\,\frac{\sin(\pi fLT)}{L\,\sin(\pi fT)}.$$ Il modulo è $\left|\frac{\sin(\pi fLT)}{L\sin(\pi fT)}\right|$, con zeri a $f=m/(LT)$, cioè multipli di $0{,}1$ Hz per $L=20$, $T=0{,}5$ s (la finestra dura $LT=10$ s). La fase è lineare: ritardo $\frac{(L-1)T}{2}=4{,}75$ s.

Effetto sul segnale. $f_1=0{,}1$ Hz è proprio il primo zero: la finestra dura esattamente un periodo di $f_1$ e la media su un periodo intero è zero, quindi la componente a $f_1$ sparisce ($|H(f_1)|\approx4\times10^{-17}$). La componente a $f_2=0{,}01$ Hz passa con guadagno $|H(f_2)|=\frac{\sin(0{,}1\pi)}{20\sin(0{,}005\pi)}=0{,}9837$ e viene ritardata di $4{,}75$ s: $$y(nT)\approx0{,}5\cdot0{,}9837\,\sin\big(2\pi f_2(nT-4{,}75)\big),\qquad n\ge L-1.$$ (Nel file del corso il commento dice "lenta + rapida", ma è $f_1=0{,}1$ a essere la più veloce.)

Transitorio. lfilter parte con lo stato nullo, cioè come se l'ingresso fosse zero per $n<0$. Per i primi $L-1=19$ campioni la media contiene anche questi zeri finti, quindi l'uscita è diversa dal regime (scarto massimo $0{,}358$); da $n=19$ in poi la finestra è piena e l'uscita coincide con la formula a $4\times10^{-16}$.

(3) Risposta in frequenza con zero-padding (parte 4)

La DFT di $h$ su $L=20$ punti darebbe solo $20$ valori di $H$. fft(h, Nfft) aggiunge zeri in coda fino a $N_{fft}=512$: la trasformata (che è sempre la stessa $H(f)$ del tempo discreto) viene valutata su una griglia più fitta, con passo $df=\frac1{N_{fft}T}=0{,}0039$ Hz. Lo zero-padding non aggiunge informazione, interpola; è lo stesso principio della ricostruzione di una sequenza finita (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 →). Poi fftshift porta $f=0$ al centro, con asse $f_k=k\,df$, $k=-N_{fft}/2,\dots,N_{fft}/2-1$, che copre un periodo, da $-1$ a $1$ Hz.

Attenzione al fattore: nel file del corso la parte 4 usa $h=b$ (cioè $1/L$) e $H=T\cdot\mathtt{fft}(h)$, che dà $H(0)=T=0{,}5$ e non $1$; con la convenzione del punto (1), $h=b/T$, si trova $H(0)=1$.

(4) Condizioni iniziali con lfiltic (parte 3, secondo blocco)

Per evitare il transitorio basta fornire a lfilter lo stato iniziale zi che corrisponde ai $L-1=19$ ingressi passati $x(-T),x(-2T),\dots,x(-19T)$ (qui noti perché l'ingresso è una somma di sinusoidi). lfiltic(b, a, y_pass, x_pass) (Matlab filtic) costruisce zi a partire da questi valori, con il più recente per primo; per un FIR non c'è uscita passata ([]). Con zi l'uscita è subito a regime: $y[0]=-0{,}1446$ (media di $x$ su $n=-19,\dots,0$) invece di $0$, e coincide con la formula in tutti i campioni a $4\times10^{-16}$.

(5) Filtro IIR passa-alto (parti 5-6)

butter(4, Wn, 'highpass') progetta un Butterworth digitale con $W_n=f_c/(F_s/2)=0{,}2/25=0{,}008$ (frequenza normalizzata rispetto a metà della frequenza di campionamento, $F_s=1/T_s=50$ Hz). Il modulo analogico è $|H(f)|=\left[1+(f_c/f)^{2\cdot4}\right]^{-1/2}$, con $|H(f_c)|=\frac1{\sqrt2}=0{,}7071$ ($-3$ dB). Coefficienti ottenuti: $$b=[0{,}96769,\ -3{,}87078,\ 5{,}80617,\ -3{,}87078,\ 0{,}96769],\qquad a=[1,\ -3{,}93433,\ 5{,}80513,\ -3{,}80723,\ 0{,}93643].$$ È un sistema con poli: la risposta impulsiva è infinita. È stabile perché i poli stanno dentro il cerchio unitario (Trasformata zeta e sistemi a tempo discretoUn sistema LTI discreto è descritto da un'equazione alle differenze Σ a_i y(n-i) = Σ b_i x(n-i) con condizioni iniziali. L'uscita è evoluzione libera (dalle condizioni iniziali) più risposta forzata g*x. Si calcola in tre modi: soluzione dell'equazione, risposta impulsiva con segnale fittizio, trasformata zeta. FIR: memoria finita, sempre stabile; IIR: stabile se tutti i poli hanno modulo minore di 1. In frequenza G(f) = Σ b_i e^{-i2πfi} / Σ a_i e^{-i2πfi}.Trasformata zeta e sistemi a tempo discreto →): moduli $0{,}9904$ (due volte) e $0{,}9771$ (due volte). Sono molto vicini a $1$, e questo, con una frequenza di taglio così bassa rispetto a $F_s$, rende lenti sia la risposta in frequenza sia il transitorio.

Effetto sul segnale. La componente a $f_1=0{,}02$ Hz (periodo $50$ s: nei $4$ s osservati è soltanto una rampa lentissima) è attenuata di un fattore $|H(f_1)|=1{,}0\times10^{-4}$ e sparisce; la componente a $f_2=1{,}6$ Hz passa con guadagno $0{,}99999997$ e una anticipazione di fase di $0{,}326$ rad (tempo $0{,}326/(2\pi\cdot1{,}6)=0{,}032$ s). In tutto l'intervallo osservato ($4$ s $=6{,}4$ periodi di $f_2$) si vede la sinusoide a $1{,}6$ Hz, ma con transitorio iniziale: all'inizio lo stato è nullo e $y[0]=b_0s[0]=0$, mentre il regime vorrebbe $\sin(0{,}326)=0{,}32$. Lo scarto dal regime è $0{,}320$ in $n=0$, $0{,}090$ in $n=25$, $0{,}0033$ in $n=50$ e al più $0{,}017$ dopo $n=100$ (la costante di tempo dei poli è $1/(1-0{,}9904)\approx104$ campioni $\approx2$ s).

(6) Risposta in frequenza dalla risposta impulsiva (parte 6)

Si genera la risposta impulsiva dandole in ingresso l'impulso ideale su $\mathbb Z(T_s)$, che vale $\frac1{T_s}$ in $0$ e $0$ altrove: impulso = (1/Ts) * [1, 0, 0, ...], h = lfilter(b, a, impulso). Poi H = Ts * fft(h, Nfft) con $N_{fft}=10000$ (passo $df=\frac1{N_{fft}T_s}=0{,}005$ Hz, periodo $50$ Hz) e fftshift per centrare l'asse; è la convenzione del punto (1) già coerente col file.

Trappola: la risposta impulsiva è infinita. Il file del corso tronca $h$ a $h_{len}=100$ campioni ($2$ s), ma la coda è ancora importante. Confronto con la risposta esatta calcolata con freqz sulla stessa griglia:

I primi $100$ campioni contengono il $99{,}8\%$ dell'energia di $h$, ma la coda residua è un andamento lento che, trasformato, pesa proprio sulle basse frequenze (quelle che il filtro dovrebbe tagliare). Il troncamento rialza $H$ in continua e abbassa la banda passante vicino a $f_c$. Servono almeno $h_{len}\approx2000$ campioni, cioè $40$ s: molte costanti di tempo dei poli. Il modulo ottenuto, trasformando la risposta (con $h_{len}$ sufficiente) è una funzione pari di $f$, periodica di periodo $50$ Hz, uguale a $0{,}7071$ in $\pm0{,}2$ Hz e a $1$ per $|f|\gg f_c$ (in particolare $|H(1{,}6)|=1$).

import numpy as np
from numpy.fft import fft, fftshift
from scipy import signal
import matplotlib.pyplot as plt

# ---------- parte 3: media mobile (FIR) con lfilter ----------
T, N, L = 0.5, 200, 20                             # passo 0.5 s, 200 campioni, finestra di 20 campioni = 10 s
t = np.arange(N) * T
f1, f2 = 0.1, 0.01
x = np.sin(2 * np.pi * f1 * t) + 0.5 * np.sin(2 * np.pi * f2 * t)
b, a = np.ones(L) / L, 1.0                         # FIR: y[n] = (1/L) sum_{k<L} x[n-k], denominatore 1
y = signal.lfilter(b, a, x)                        # MATLAB: filter(b, a, x)

H2 = abs(np.sin(np.pi * f2 * L * T) / (L * np.sin(np.pi * f2 * T)))      # |H(f2)| teorico
ritardo = (L - 1) * T / 2                          # ritardo di gruppo del FIR a fase lineare: 4.75 s
y_reg = 0.5 * H2 * np.sin(2 * np.pi * f2 * (t - ritardo))                 # la sola componente lenta
print(f"|H(f2)| = {H2:.4f}, ritardo = {ritardo} s")
print("errore a regime (n >= L-1):", np.abs(y[L - 1:] - y_reg[L - 1:]).max(), "; nel transitorio:", np.abs(y[:L - 1] - y_reg[:L - 1]).max())

# ---------- parte 4: risposta in frequenza con zero-padding ----------
h = b / T                                          # risposta impulsiva discreta del corso: g(kT) = b_k / T
Nfft = 512
H = fftshift(T * fft(h, Nfft))                     # fft(h, Nfft) allunga con zeri; fattore T = peso di Haar
df = 1 / (Nfft * T)                                # 0.0039 Hz; periodo della risposta: Fp = 1/T = 2 Hz
f = np.arange(-Nfft // 2, Nfft // 2) * df
H_teo = lambda ff: np.abs(np.sin(np.pi * ff * L * T) / (L * np.sin(np.pi * ff * T) + 1e-300))
print("H(0) =", abs(H[Nfft // 2]), " (se si usasse h = b, senza 1/T: ", abs(T * fft(b, Nfft)[0]), ")")
for ff in (0.01, 0.05, 0.1, 0.15, 0.2):
    i = np.argmin(np.abs(f - ff))
    print(f"  |H({f[i]:.4f})| = {abs(H[i]):.4f}  teorico {H_teo(f[i]):.4f}")
print("zeri a multipli di 1/(L T) =", 1 / (L * T), "Hz; massimo del primo lobo secondario:",
      abs(H)[(f > 0.1) & (f < 0.2)].max())

# ---------- condizioni iniziali con lfiltic ----------
k = np.arange(1, L)                                # il FIR ha memoria di L-1 = 19 ingressi passati x(-T), x(-2T), ...
x_pass = np.sin(2 * np.pi * f1 * (-k * T)) + 0.5 * np.sin(2 * np.pi * f2 * (-k * T))
zi = signal.lfiltic(b, [1.0], [], x_pass)          # MATLAB: filtic(b, a, [], x_pass)
y_zi, _ = signal.lfilter(b, [1.0], x, zi=zi)       # MATLAB: filter(b, a, x, zi)
print("senza zi: y[0] =", y[0], "; con zi: y[0] =", y_zi[0], "; errore a regime di y_zi (tutti gli n):",
      np.abs(y_zi - y_reg).max())

# ---------- parte 5: filtro IIR di Butterworth passa-alto ----------
Ts, N = 0.02, 200                                  # Fs = 50 Hz, durata 4 s
t = np.arange(N) * Ts
f1, f2 = 0.02, 1.6                                 # componente molto lenta + componente a 1.6 Hz
s = np.sin(2 * np.pi * f1 * t) + np.sin(2 * np.pi * f2 * t)
fc, Fs = 0.2, 1 / Ts
Wn = fc / (Fs / 2)                                 # frequenza di taglio normalizzata a Fs/2: 0.008
bb, aa = signal.butter(4, Wn, "highpass")          # MATLAB: butter(4, Wn, 'high')
y_iir = signal.lfilter(bb, aa, s)
print("Wn =", Wn, "\n b =", np.round(bb, 5), "\n a =", np.round(aa, 5), "\n |poli| =", np.round(np.abs(np.roots(aa)), 5))
_, Hz = signal.freqz(bb, aa, worN=[f1, fc, f2], fs=Fs)
print("|H| in f1, fc, f2:", np.abs(Hz), "; fase in f2:", np.angle(Hz[2]))
y_reg = np.abs(Hz[0]) * np.sin(2 * np.pi * f1 * t + np.angle(Hz[0])) + np.abs(Hz[2]) * np.sin(2 * np.pi * f2 * t + np.angle(Hz[2]))
e = np.abs(y_iir - y_reg)
print("scarto dal regime: n=0", e[0], "; n=25", e[25], "; n=50", e[50], "; max per n>=100:", e[100:].max())

# ---------- parte 6: risposta in frequenza dalla risposta impulsiva ----------
Nfft = 10000
df = 1 / (Nfft * Ts)                               # 0.005 Hz
f = np.arange(-Nfft // 2, Nfft // 2) * df
_, H_ex = signal.freqz(bb, aa, worN=f, fs=Fs)      # riferimento esatto
for h_len in (100, 500, 2000):
    impulso = (1 / Ts) * np.r_[1.0, np.zeros(h_len - 1)]    # impulso ideale su Z(Ts): 1/Ts in 0
    h = signal.lfilter(bb, aa, impulso)                     # risposta impulsiva (troncata a h_len campioni)
    H = fftshift(Ts * fft(h, Nfft))                         # zero-padding fino a Nfft, fattore Ts
    i = np.argmin(np.abs(f - 0.2))
    print(f"h_len = {h_len:4d}: |H(0.2)| = {abs(H[i]):.4f} (esatto {abs(H_ex[i]):.4f}); |H(0)| = {abs(H[Nfft // 2]):.2e}; "
          f"errore max |f|<5: {np.abs(H - H_ex)[np.abs(f) < 5].max():.2e}")

fig, ax = plt.subplots(2, 1)
ax[0].stem(t, s); ax[0].stem(t, y_iir, "r")
ax[1].plot(f, np.abs(H)); ax[1].set_xlim(-5, 5)
plt.show()

Cosa si vede nelle figure. Media mobile: l'ingresso è una sinusoide da $0{,}1$ Hz (periodo $10$ s) sovrapposta a una più lenta da $0{,}01$ Hz; l'uscita è una sinusoide lenta di ampiezza $\approx0{,}49$ che parte da $0$ con una rampa nei primi $9{,}5$ s (transitorio). Nel grafico di $\lvert H\rvert$ si vede il lobo principale largo $\pm0{,}1$ Hz, lobi secondari bassi (il primo vale $0{,}219$) e la ripetizione con periodo $2$ Hz. Butterworth: l'uscita è la sinusoide da $1{,}6$ Hz quasi senza distorsione dopo $\approx1$ s; il modulo della risposta in frequenza è nullo vicino a $0$, sale con pendenza ripida attorno a $0{,}2$ Hz e vale $1$ per $|f|$ oltre circa $0{,}4$ Hz.

Errori tipici.

• Dimenticare il fattore $T$ (o $T_s$) quando si trasforma la risposta impulsiva: $H=T\cdot\mathtt{fft}(g)$ con $g=b/T$; senza di esso il guadagno in continua è sbagliato.
• Troncare la risposta impulsiva di un IIR: non è finita, e un taglio troppo corto falsa proprio le basse frequenze.
• Confondere la frequenza di taglio $f_c$ con quella normalizzata $W_n=f_c/(F_s/2)$ in butter.

Vedi anche: Esercizio - filtraggio di un impulso rettangolare con convoluzione e trasformate in Python (laboratorio 5), 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 →.