Esercizio - Antitrasformata zeta con Python (lezione 10)

Testo (dispensa "Focus on z anti-transform" della lezione 10 del corso Multimedia Signal Processing, UniPD, con il file MATLAB z_transform_exercises.m; qui rifatto in Python). Quattro esercizi: (1) antitrasformata di $X(z)=\frac1{1-0{,}5z^{-1}}$ e verifica numerica; (2) frazioni parziali di $X(z)=\frac1{(1-z^{-1})(1-0{,}5z^{-1})}$; (3) antitrasformata di un FIR con la divisione lunga; (4) diagramma poli-zeri e stabilità di $H(z)=\frac{1-z^{-1}}{1-0{,}9\cos(\pi/4)z^{-1}+0{,}81z^{-2}}$. In MATLAB si usano impz, residuez e zplane; in Python scipy.signal.lfilter, residuez e tf2zpk.

Teoria usata: Antitrasformata zetaPer tornare da X(z) a x[n] serve anche la ROC. Formula di inversione x[n] = (1/2πj)∮ X(z)z^{n-1}dz (somma dei residui). Nella pratica: ispezione (riconoscere coppie note, scegliendo la ROC), divisione lunga (sviluppo in serie di potenze), frazioni parziali per le razionali. Causale con poli semplici p_k: X = Σ A_k/(1-p_k z^-1) con A_k = [(1-p_k z^-1)X]_{z=p_k}, x[n] = Σ A_k p_k^n u[n]; se il grado del numeratore in z^-1 è ≥ di quello del denominatore c'è anche una parte polinomiale (impulsi). Poli doppi: n a^n u[n]; poli complessi coniugati: r^n cos(ω0 n + φ).Antitrasformata zeta → (ispezione, frazioni parziali, divisione lunga), Funzione di sistema, poli, zeri e stabilitàLa funzione di sistema H(z) è la trasformata zeta della risposta impulsiva: con ingresso z^n l'uscita è H(z) z^n. Per un FIR H(z) = Σ b_k z^{-k} è un polinomio con M zeri e M poli in z = 0; in generale H = B(z)/A(z) dall'equazione alle differenze. Sulla circonferenza unitaria H(e^{jω̂}) è la risposta in frequenza: |H| = prodotto delle distanze dagli zeri / prodotto delle distanze dai poli, quindi gli zeri bloccano frequenze e i poli le esaltano. Un LTI causale è BIBO stabile se e solo se tutti i poli hanno modulo < 1 (a meno di cancellazioni polo-zero); i FIR sono sempre stabili.Funzione di sistema, poli, zeri e stabilità → (poli, zeri, stabilità), Trasformata zeta - definizione e regione di convergenzaLa trasformata zeta bilatera X(z) = Σ x[n] z^{-n} associa a una sequenza una funzione della variabile complessa z, definita nella regione di convergenza (ROC), sempre una corona circolare |z| in (R1, R2). Segnale a durata finita: ROC tutto il piano (tranne eventualmente 0 e ∞); causale: |z| > R1 (teorema di Abel); anticausale: |z| < R2; bilatero: intersezione, se non vuota. La stessa espressione algebrica con ROC diverse è la trasformata di segnali diversi: la ROC fa parte della trasformata. Sulla circonferenza unitaria, se è nella ROC, X(e^{jθ}) è la trasformata di Fourier. Le ROC non contengono poli.Trasformata zeta - definizione e regione di convergenza →.

Traduzione degli strumenti

Una funzione razionale $B(z)/A(z)$ si rappresenta con i due vettori dei coefficienti in potenze di $z^{-1}$, in ordine crescente (b = [b0, b1, ...], a = [1, a1, ...]). La risposta impulsiva causale è la soluzione della ricorsione con ingresso $\delta[n]$, quindi coincide con l'antitrasformata:

import numpy as np
from scipy.signal import lfilter, residuez, tf2zpk

def impz(b, a, N):
    """Risposta impulsiva (antitrasformata causale) di B(z)/A(z)."""
    d = np.zeros(N); d[0] = 1
    return lfilter(b, a, d)

residuez(b, a) dà residui $r_k$, poli $p_k$ e la parte diretta $k$ (i coefficienti degli impulsi), come in MATLAB; tf2zpk(b, a) dà zeri e poli.

Esercizio 1 - un polo

$X(z)=\sum_{n\ge0}(0{,}5)^nz^{-n}=\frac1{1-0{,}5z^{-1}}$, $\lvert z\rvert>0{,}5$ (geometrica di ragione $0{,}5z^{-1}$). La coppia $a^nu[n]\leftrightarrow\frac1{1-az^{-1}}$ ($\lvert z\rvert>\lvert a\rvert$) con $a=0{,}5$ restituisce $x[n]=(0{,}5)^nu[n]$: $1;\ 0{,}5;\ 0{,}25;\ 0{,}125;\dots\to0$. La sequenza decade perché il polo è dentro il cerchio unitario ($\lvert a\rvert<1$).

n = np.arange(20)
h = impz([1], [1, -0.5], 20)
print(np.max(np.abs(h - 0.5**n)))    # 0.0

impz risolve $y[n]=0{,}5y[n-1]+\delta[n]$ campione per campione, ed è esattamente $0{,}5^n$ (errore $0$ in aritmetica double).

Esercizio 2 - frazioni parziali

$X=\dfrac1{(1-z^{-1})(1-0{,}5z^{-1})}=\dfrac1{1-1{,}5z^{-1}+0{,}5z^{-2}}$: poli semplici in $z=1$ e $z=0{,}5$, causale. $A=\big[(1-z^{-1})X\big]_{z=1}=\frac1{1-0{,}5}=2$ e $B=\big[(1-0{,}5z^{-1})X\big]_{z=0{,}5}=\frac1{1-2}=-1$. Quindi $$X(z)=\frac2{1-z^{-1}}-\frac1{1-0{,}5z^{-1}},\qquad x[n]=\big[2-(0{,}5)^n\big]u[n]$$ (tabella: $x[0]=1$, $x[1]=1{,}5$, $x[2]=1{,}75$, $x[3]=1{,}875$, $\dots\to2$: il polo in $1$ dà il valore costante a regime). Il polo in $z=1$ sta sul cerchio unitario: se questa $X$ fosse la funzione di sistema di un filtro, il filtro non sarebbe BIBO stabile ($x[n]\to2\neq0$, $\sum\lvert x[n]\rvert$ diverge).

b2, a2 = [1], [1, -1.5, 0.5]
r, p, k = residuez(b2, a2)        # r = [-1, 2], p = [0.5, 1], k = []
n = np.arange(25)
x_pf = np.real(r[0]*p[0]**n + r[1]*p[1]**n)       # sum_k r_k p_k^n
print(np.max(np.abs(impz(b2, a2, 25) - x_pf)))     # 0.0

Non c'è parte diretta $k$ perché il numeratore ($1$) ha grado minore del denominatore. La ricostruzione $\sum_kr_kp_k^n$ coincide con la risposta impulsiva (errore $0$) e con $2-0{,}5^n$.

Esercizio 3 - FIR e divisione lunga

$X(z)=1+2z^{-1}+3z^{-2}+2z^{-3}+z^{-4}$ è un polinomio in $z^{-1}$: confrontando con $\sum x[n]z^{-n}$ i coefficienti sono i campioni, $x=\{1,2,3,2,1\}$ per $n=0,\dots,4$ e $0$ altrove. La divisione lunga è banale (il denominatore è $1$). Per un IIR invece non termina: $\frac1{1-0{,}5z^{-1}}=1+0{,}5z^{-1}+0{,}25z^{-2}+\dots$ non si ferma mai, ed è per questo che per gli IIR si preferiscono le frazioni parziali.

print(impz([1, 2, 3, 2, 1], [1], 5))    # [1. 2. 3. 2. 1.]

Differenza tra FIR e IIR: il FIR ha risposta impulsiva finita (la divisione termina), poli solo in $z=0$ e stabilità sempre; l'IIR ha poli non nulli e la stabilità dipende dalla loro posizione.

Esercizio 4 - poli, zeri e stabilità

$H(z)=\dfrac{1-z^{-1}}{1-0{,}9\cos(\frac\pi4)\,z^{-1}+0{,}81z^{-2}}$.

Zeri. $1-z^{-1}=0\Rightarrow z=1$: zero in $\hat\omega=0$ (blocca la continua; il numeratore è la differenza prima).

Poli. $z^2-0{,}9\cos\frac\pi4\,z+0{,}81=0$, cioè $z^2-0{,}6364z+0{,}81=0$. Il prodotto delle radici è $0{,}81$ (modulo $0{,}9$ per ciascuna, essendo complesse coniugate) e la somma è $0{,}6364=2\cdot0{,}9\cos\omega_0$, quindi $\cos\omega_0=\frac{0{,}9\cos(\pi/4)}{2\cdot0{,}9}=0{,}3536$ e $\omega_0=1{,}2094$ rad $=69{,}3^\circ$. I poli sono $0{,}318\pm0{,}842j=0{,}9e^{\pm j1{,}2094}$.

Errore nella dispensa del corso. La dispensa scrive "$z_{1,2}=0{,}9e^{\pm j\pi/4}$" (angolo $45^\circ$). Ma $0{,}9e^{\pm j\pi/4}$ ha somma $2\cdot0{,}9\cos\frac\pi4=1{,}273$, mentre il polinomio ha somma $0{,}636$: il coefficiente $0{,}9\cos\frac\pi4$ vale $2r\cos\omega_0$ con $r=0{,}9$ solo se $\cos\omega_0=\frac12\cos\frac\pi4$. Il modulo $0{,}9$ è giusto e quindi anche la conclusione sulla stabilità; la frequenza di risonanza è però $\approx1{,}21$ rad/campione, non $\frac\pi4=0{,}785$ (con il calcolo numerico: massimo di $\lvert H\rvert$ in $\hat\omega\approx1{,}215$, vedi Funzione di sistema, poli, zeri e stabilitàLa funzione di sistema H(z) è la trasformata zeta della risposta impulsiva: con ingresso z^n l'uscita è H(z) z^n. Per un FIR H(z) = Σ b_k z^{-k} è un polinomio con M zeri e M poli in z = 0; in generale H = B(z)/A(z) dall'equazione alle differenze. Sulla circonferenza unitaria H(e^{jω̂}) è la risposta in frequenza: |H| = prodotto delle distanze dagli zeri / prodotto delle distanze dai poli, quindi gli zeri bloccano frequenze e i poli le esaltano. Un LTI causale è BIBO stabile se e solo se tutti i poli hanno modulo < 1 (a meno di cancellazioni polo-zero); i FIR sono sempre stabili.Funzione di sistema, poli, zeri e stabilità →).

Stabilità. Sistema causale con poli $\lvert p\rvert=0{,}9<1$: stabile (verificato: $\sum\lvert h\rvert\approx8{,}1$ su $40$ campioni, che converge).

Risposta impulsiva. Con poli $re^{\pm j\omega_0}$, $h[n]$ è una sinusoide smorzata: scrivendo $g[n]=r^n\frac{\sin((n+1)\omega_0)}{\sin\omega_0}u[n]$ per $\frac1{1-2r\cos\omega_0z^{-1}+r^2z^{-2}}$ e usando $H=(1-z^{-1})G$, si ha $h[n]=g[n]-g[n-1]$. Valori: $1;\ -0{,}364;\ -1{,}041;\ -0{,}368;\ 0{,}609;\ 0{,}686;\ -0{,}057;\ -0{,}592;\dots$ (verificato con lfilter), con oscillazione di periodo $\frac{2\pi}{1{,}209}\approx5{,}2$ campioni.

b4, a4 = [1, -1], [1, -0.9*np.cos(np.pi/4), 0.81]
z, p, kk = tf2zpk(b4, a4)
print(z, p, np.abs(p), np.degrees(np.angle(p)))
# [1.] [0.318+0.842j 0.318-0.842j] [0.9 0.9] [ 69.295 -69.295]

Lettura del diagramma. Numeratore $\to$ zeri $\to$ frequenze bloccate; denominatore $\to$ poli $\to$ frequenze esaltate; tutti i poli con $\lvert p\rvert<1$ $\to$ stabile.