Prenosne funkcije digitalnih sistema i osnovni digitalni filtri

Digitalna obrada signala, Vladimir Petrović

U ovom notebook-u ćemo uvesti alate za analizu digitalnih sistema. Definisaćemo prenosnu funkciju sistema, nacrtaćemo frekvencijsku karakteristiku sistema i naći odziv sistema na proizvoljnu pobudu. Na kraju ćemo projektovati najjednostavnije filtre sa beskonačnim impulsnim odzivom, analizirati uticaj rasporeda nula i polova prenosne funkcije na odziv i frekvencijsku karakteristiku filtra.

Funkcija prenosa digitalnog filtra

Ako je impulsni odziv filtra konačan, onda se odziv na proizvoljnu pobudu može izračunati konvolucijom. Prenosna funkcija tog sistema predstavlja Z-transformaciju impulsnog odziva:

$$ H(z) = \mathcal{Z} \left\lbrace h[n] \right\rbrace = \sum_{n=0}^{N - 1} h[n] z^{-n}. $$

Ovakve sisteme nazivamo sistemima sa konačnim impulsnim odzivom (Finite Impulse Response - FIR). Međutim, ako je impulsni odziv beskonačan, tj. ako je beskonačan broj odbiraka impulsnog odziva različit od nule, obično se takvi sistemi ne predstavljaju impulsnim odzivom, već diferencnom jednačinom. Na primer:

$$ y[n] - \frac{5}{2}y[n-1] + y[n-2] = 3x[n] - 7x[n-1] + 5x[n-2]. $$

Korišćenjem osobine Z-transformacije da je $\mathcal{Z} \left\lbrace x[n-m] \right\rbrace = y^{-m}X(z)$, možemo doći do prenosne funkcije ovako definisanog sistema:

$$ \mathcal{Z} \left\lbrace y[n] - \frac{5}{2}y[n-1] + y[n-2] \right\rbrace = \mathcal{Z} \left\lbrace 3x[n] - 7x[n-1] + 5x[n-2] \right\rbrace \\ Y(z) - \frac{5}{2}z^{-1}Y(z) + z^{-2}Y(z) = 3X(z) - 7X(z) + 5X(z)\\ \Downarrow \\ H(z) = \frac{Y(z)}{X(z)} = \frac{3 - 7z^{-1} + 5z^{-2}}{1 - \frac{5}{2}z^{-1} + z^{-2}} = \frac{3z^2 - 7z + 5}{z^2 - \frac{5}{2}z + 1} $$

Ovakve sisteme nazivamo sistemima sa beskonačnim impulsnim odzivom (Infinite Impulse Response - IIR). Korene polinoma u brojiocu nazivamo nulama sistema, dok korene polinoma u imeniocu nazivamo polovima sistema:

$$ H(z) = \frac{3z^2 - 7z + 5}{z^2 - \frac{5}{2}z + 1} = k \frac{(z - z_{o1})(z - z_{o2})}{(z - z_{p1})(z - z_{p2})}. $$

Nule funkcije prenosa su obeležene sa $z_{oi}$, a polovi sa $z_{oi}$. Konstanta $k$ predstavlja pojačanje.

Ako je kompleksna promenljiva $z$ jednaka nekoj od nula sistema, očigledno, prenosna funkcija ima vrednost $0$, pa odatle i ime. Ako je kompleksna premenljiva $z$ jednaka nekom od polova sistema, očigledno, vrednost prenosne funkcije teži beskonačnosti. Za funkciju prenosa kažemo da je N-tog reda, ako je polinom u imeniocu funkcije prenosa N-tog reda.

Nađimo nule, polove i konstantu pojačanja, a zatim i nacrtajmo raspored nula i polova funkcije prenosa date sledećim izrazom:

$$ Y(z) = \frac{z^2 + 2z + 1}{z^2 - 1,1421z + 0,41421}. $$

zplane je funkcija napisana u fajlu dosutils.py. Možete je proučiti. Ona crta jedinični krug i položaj nula i polova. U slučaju da je funkcija prenosa više nula ili polova u istoj tački, pored odgovarajućeg markera se ispisuje broj nula ili polova u tim tačkama (Napomena: Usled konačne preciznosti izračunavanja, može se desiti da se dve vrednosti u istoj tački na realnoj osi protumače kao $a+jb$ i $a-jb$ (ili $-a+jb$ i $a+jb$ na imaginarnoj osi), ako je $b$ (odnosno $a$, za imaginarnu osu) mnogo mala vrednost, pa u tim slučajevima funkcija ne napiše broj istih polova jer ih smatra različitim). Takođe, funkcija vraća nule, polove i pojačanje k, kao i objekat tipa AxesSubplot na kome su nacrtani svi elementi prikaza. Funkciji se može, kao treći argument, proslediti drugi objekat tipa AxesSubplot na koji će se izvršiti crtanje.

Koreni polinoma se lako nalaze korišćenjem Numpy funkcije roots.

Odziv sistema na pobudu se dobija Scipy funkcijom lfilter. Frekvencijska karakteristika filtra se može odrediti korišćenjem Scipy funkcije freqz. Nacrtajmo impulsni odziv gorenavedenog sistema i njevu frekvencijsku karakteristiku.

Primećujemo da frekvencijska karakteristika ovog filtra ima velike vrednosti za male učestanosti, a male vrednosti za velike učestanosti. To znači da je ovaj filtar filtar propusnik niskih učestanosti. Stoga, izgenerišimo reprezentativan signal i isfiltrirajmo ga. Reprezentativan signal bi mogao da bude zbir tri sinusoidalna signala, na primer:

$$ x[n] = \cos(2\pi F_1 n) + \cos(2\pi F_2 n) + \cos(2\pi F_3 n), $$

gde su $F_1 = 0,05$, $F_2 = 0,075$ i $F_1 = 0,3$.

Veza između Furijeove i Z transformacije

Ponovimo izraze za Furijeovu i Z transformaciju signala $x[n]$ konačnog trajanja:

$$ X(e^{j\Omega}) = \mathcal{F} \left\lbrace x[n] \right\rbrace = \sum_{n=0}^{N - 1} x[n] e^{-jn\Omega}, $$$$ X(z) = \mathcal{Z} \left\lbrace x[n] \right\rbrace = \sum_{n=0}^{N - 1} x[n] z^{-n}. $$

Primećujemo da je Furijeova transformacija samo specijalan slučaj Z transformacije kada kompleksna promenljiva $z$ može uzeti vrednosti samo na jedičnom krugu u kompleksnoj ravni, tj. $z = e^{j\Omega}$. Dakle kompleksna promenljiva $z$ se lako može povezati sa relativnim kružnim učestanostima signala. Tako, ako na primer imamo Z transformaciju nekog signala, vrlo lako, smenom, možemo doći do njegove frekvencijske karaktiristike. Evo primera:

$$ x[n] = u[n] - u[n - 6]\\ X(z) = \sum_{n=0}^{N - 1} x[n] z^{-n} = \sum_{n=0}^{5} z^{-n} = \frac{1 - z^{-6}}{1 - z^{-1}}\\ \Downarrow \\ X(e^{j\Omega}) = X(z)|_{z = e^{j\Omega}} = \frac{1 - e^{-j6\Omega}}{1 - e^{-j\Omega}} = e^{-j5\Omega/2} \frac{e^{j3\Omega} - e^{-j3\Omega}}{e^{j\Omega/2} - e^{-j\Omega/2}} = e^{-j5\Omega/2} \frac{\sin (3\Omega)}{\sin(\Omega/2)} $$

Kada imamo Furijeovu transformaciju, onda je jednostavno dobiti i diskretnu Furijeovu transformaciju, odabiranjem Furijeove transformacije:

$$ X[k] = X(e^{j2\pi k/N}) = X(z)|_{z = e^{j2\pi k/N}} = e^{-j5\pi k/N} \frac{\sin (6\pi k / N)}{\sin(\pi k /N)} $$

Obrnuto, ako je poznata diskretna Furijeova transformacija, iz njenih odbiraka se može odrediti Z transformacija na sledeći način:

$$ X(z) = \mathcal{Z} \left\lbrace x[n] \right\rbrace = \sum_{n=0}^{N - 1} x[n] z^{-n} = \sum_{n=0}^{N - 1} \left( \frac{1}{N} \sum_{k=0}^{N - 1} X[k] e^{j2\pi k n/N} \right) z^{-n} = \frac{1}{N} \sum_{k=0}^{N - 1} X[k] \sum_{n=0}^{N - 1} \left( e^{j2\pi k /N} z^{-1} \right)^n \\ \Downarrow \\ X(z) = \frac{1-z^{-N}}{N} \sum_{k=0}^{N - 1} \frac{X[k]}{1 - e^{j2\pi k /N} z^{-1} } $$

Da li DFT može da se koristi za iračunavanje odziva IIR sistema?

Ako imamo funkciju prenosa sistema $H(z)$, onda znamo da je Z transformacija odziva sistema na pobudu $x[n]$

$$ Y(z) = H(z)X(z) = \frac{B(z)}{A(z)}X(z) = \frac{b_0 + b_1 z^{-1} + ... + b_m z^{-m}}{a_0 + a_1 z^{-1} + ... + a_n z^{-n}} X(z) $$

S obzirom na vezu između Z tranformacije i diskretne Furijeove transformacije možemo napisati sledeći izraz. $$ Y(e^{j2\pi k/N}) = H(e^{j2\pi k/N})X(e^{j2\pi k/N}) = \frac{B(e^{j2\pi k/N})}{A(e^{j2\pi k/N})}X(e^{j2\pi k/N})\\ \Downarrow\\ Y[k] = \frac{B[k]}{A[k]} X[k] $$

Postavlja se pitanje da li je onda odziv sistema zapravo $y[n] = IDFT\left\lbrace Y[k] \right\rbrace$, tj. da li se odziv sistema može izračunati korišćenjem sledećeg izraza:

$$ y[n] = IFFT \left( \frac{FFT(b, N)}{FFT(a, N)} FFT(x, N) \right). $$

Najlakše je ovo pokazati na primeru. Najpre treba napomenuti da Scipy funkcija lfilter izračunava odziv egzaktno korišćenjem diferencne jednačine, pa taj odziv možemo smatrati tačnim. U narednom primeru treba probati različite vrednosti za dužinu ulaznog signala i uporediti rezultate.

Primećujemo da odziv nije isti i to najpre u početnom delu. Moramo napomenuti da je ovaj sistem sistem sa beskonačnim impulsnim odzivom i da ne postoji način da se impulsni odziv produži nulama do dužine rezultata konvolucije kao što je to bio slučaj kod sistema sa konačnim impulsnim odzivom. S obzirom na to da je DFT izračunat u konačnom broju tačaka, a da je odziv beskonačan, dobija se preklapanje u vremenskom domenu na početku odziva. Međutim, impulsni odziv stabilnih IIR sistema, iako beskonačnog trajanja, ima veoma male vrednosti za veliko $n$. Zbog toga je ovo preklapanje izraženo samo na početku, pa ipak DFT možemo iskorititi za estimaciju odziva. Ako želimo izbeći ovako značajno preklapanje, ulazni signal se mora pordužiti nulama bar za onoliko odbiraka koliko odbiraka impulsnog odziva ima značajne vrednosti.

Korišćenje DFT-a za izračunavanje odziva IIR sistema se relativno retko radi jer nema značajnih prednosti. IIR sistemi obično nisu prevelikog reda, pa je izračunavanje diferencne jednačine zapravo brži i tačniji način da se dobije odziv sistema. Naravno, sve ranije što smo učili o sistemima sa konačnim impulsnim odzivom velikog trajanja i dalje važi i DFT tu ima svoju primenu.

Uticaj položaja nula i polova na impulsni odziv i frekvencijsku karakteristiku

Stabilnost sistema

Jedan od najvažnijih postupaka prilikom projektovanja IIR sistema je analiza njihove stabilnosti. Za LTI sistem važi da je stabilan ako za njegov impulsni odziv važi sledeća nejednakost: $$ \sum_{n=-\infty}^{\infty}|h[n]| < \infty. $$

Funkcija prenosa se, kao što je ranije rečeno, može definisati kao Z transformacija impulsnog odziva: $$ H(z) = \sum_{n=-\infty}^{\infty} h[n]z^{-n}. $$

Posmatrajući prethodni izraz se vidi da je uslov stabilnosti zadovoljen ako jedinični krug u kompleksnoj ravni $|z| = 1$ pripada oblasti konvergencije funkcije $H(z)$. Drugim rečima, linearan vremenski nepromenljiv diskretni sistem je stabilan ako i samo ako oblast konvergencije funkcije prenosa sistema uključuje jedinični krug.

Od posebnog interesa su kauzalni sistemi kod kojih je oblast konvergencije izvan kruga na kome se nalazi najudaljeniji pol od kordinatnog početka: $|z| > |z_{pmax}|$. Zbog toga možemo da kažemo da je kauzalni linearni vremenski nepromenljivi diskretni sistem stabilan ako i samo ako svi polovi funkcije prenosa $H(z)$ leže unutar jediničnog kruga.

Funkcija prenosa 2. reda, bez nula

Neka je data funkcija prenosa drugog reda sa parom konjugovano kompleksnih polova i bez nula: $$ H(z) = \frac{1}{(z - z_{p1})(z - z_{p2})}, $$ gde su $z_{p1,2} = r e^{\pm j\theta}$. Posmatrajmo kako izgledaju impulsni odziv i frekvencijska karakteristika ako menjamo položaj ovih polova.

Primećujemo da ugao $\theta$ određuje učestanost oscilacija u impulsnom odzivu, ali i položaj najveće vrednosti u amplitudskoj karakteristici. Poluprečnik polova $r$ određuje stabilnost sistema. Ako je veći od 1, sistem je nestabilan, ako je manji od 1 sistem je stabilan, a ako je jednak 1, sistem je granično stabilan. Što su polovi bliži jediničnom krugu to oscilacije impulsnog odziva duže traju, ali to povlači oštrije prelaze u frekvencijskoj karakteristici.

Funkcija prenosa 2. reda - filtar propusnik niskih učestanosti

Imajući u vidu sve navedeno, možemo pokušati da projektujemo filtre direktno postavljanjem polova i nula na odgovarajuće mesto. Videli smo da je veza kompleknse promenljive $z$ i učestanosti $\Omega$ data sa $z = e^{j\Omega}$. To znači da ako želimo da potisnemo frekvencijske komponente na visokim učestanostima, funkcija prenosa treba da ima nule u $z = e^{j\Omega}|_{\Omega = \pi} = e^{j\pi} = -1$.

Neka je data funkcija prenosa drugog reda sa parom konjugovano kompleksnih polova i dve nule u $z = -1$: $$ H_{LP}(z) = k\frac{(1 + z^{-1})^2}{(1 - r e^{j\theta} z^{-1})(1 - r e^{-j\theta} z^{-1})} = k\frac{(1 + z^{-1})^2}{1- 2r \cos(\theta) z^{-1} + r^2 z^{-2}} $$ Posmatrajmo kako izgledaju impulsni odziv i frekvencijska karakteristika ako menjamo položaj polova.

Funkcija prenosa 2. reda - filtar propusnik visokih učestanosti

Slično kao i ranije, ako želimo da potisnemo frekvencijske komponente na niskim učestanostima, funkcija prenosa treba da ima nule u $z = e^{j\Omega}|_{\Omega = 0} = e^{j0} = 1$.

Tada je funkcija prenosa drugog reda sa parom konjugovano kompleksnih polova i dve nule u $z = 1$: $$ H_{HP}(z) = k\frac{(1 - z^{-1})^2}{(1 - r e^{j\theta} z^{-1})(1 - r e^{-j\theta} z^{-1})} = k\frac{(1 - z^{-1})^2}{1- 2r \cos(\theta) z^{-1} + r^2 z^{-2}} $$ Posmatrajmo kako izgledaju impulsni odziv i frekvencijska karakteristika ako menjamo položaj polova.

Funkcija prenosa 2. reda - filtar propusnik opsega učestanosti

Konačno, ako želimo da potisnemo frekvencijske komponente i na niskim i na visokim učestanostima, a zadržimo komponente na srednjim učestanostima, funkcija prenosa treba da ima nulu i u $z = -1$ i u $z = 1$.

Tada je funkcija prenosa drugog reda sa parom konjugovano kompleksnih polova i sa po jednom nulom u $z = -1$ i $z = 1$: $$ H_{BP}(z) = k\frac{(1 + z^{-1})(1 - z^{-1})}{(1 - r e^{j\theta} z^{-1})(1 - r e^{-j\theta} z^{-1})} = k\frac{(1 + z^{-1})(1 - z^{-1})}{1- 2r \cos(\theta) z^{-1} + r^2z^{-2}} $$ Posmatrajmo kako izgledaju impulsni odziv i frekvencijska karakteristika ako menjamo položaj polova.

Funkcija prenosa 2. reda - filtar nepropusnik opsega učestanosti

Dodatak na ovu temu je filtar nepropusnik opsega učestanosti. Nule ovog filtra treba da postavimo na učestanosti na kojima želimo da potisnemo frekvencijske komponente.

Tada je funkcija prenosa drugog reda sa parom konjugovano kompleksnih polova i parom konjugovano kompleksnih nula: $$ H_{BS}(z) = k\frac{(1 - e^{j\theta_z} z^{-1})(1 - e^{-j\theta_z} z^{-1})}{(1 - r e^{j\theta_p} z^{-1})(1 - r e^{-j\theta_p} z^{-1})} = k\frac{1- 2 \cos(\theta_z) z^{-1} + z^{-2}}{1- 2r \cos(\theta_p) z^{-1} + r^2z^{-2}} $$ Posmatrajmo kako izgledaju impulsni odziv i frekvencijska karakteristika ako menjamo položaj nula i polova.

Literatura