%matplotlib inline
import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt
from matplotlib import patches

from util import set_style
set_style()

# częstotliwość próbkowania
fs = 48000

Filtry o nieskończonej odpowiedzi impulsowej (IIR)¶

W przeciwieństwie do filtrów FIR, które są wyłącznie cyfrowym algorytmem, filtry o nieskończonej odpowiedzi impulsowej (NOI), określane angielskim skrótem IIR (Infinite Impulse Filter), są cyfrowym odpowiednikiem analogowych filtrów, budowanych w oparciu o rezystory, kondensatory i wzmacniacze operacyjne. Metody projektowania filtrów analogowych są szeroko opracowane i metody cyfrowego projektowania filtrów IIR w dużej mierze na nich bazują.

Cechą odróżniającą filtry IIR od FIR jest to, że są one rekursywne, tzn. posiadają sprzężenie zwrotne. W filtracji biorą udział nie tylko próbki sygnału, ale także poprzednie wyniki filtracji (wartości wyjściowe). Filtr rzędu $N$ używa do filtracji: bieżącej próbki sygnału, $N$ poprzednich próbek sygnału oraz $N$ poprzednich wyników filtracji. Równanie różnicowe filtru IIR wygląda następująco:

$$y(n) = b_0 x(n) + b_1 x(n-1) + b_2 x(n-2) ... + b_N x(n-N) - a_1 y(n-1) - a_2 y(n-2) - ... - a_N y(n-N)$$

Transmitancja filtru IIR jest następująca:

$$H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + ... + b_N z^{-N}}{1 + a_1 z^{-1} + a_2 z^{-2} + ... + a_N z^{-N}}$$

Współczynniki filtru $a$ (w mianowniku) nazywa się rekursywnymi, współczynniki $b$ (w liczniku) nazywa się nierekursywnymi.

Projektowanie filtrów IIR ¶

Aby zaprojektować filtr IIR w Pythonie, użyjemy modułu scipy.signal, a z niego funkcji iirfilter. Sposób projektowania jest nieco inny niż w przypadku filtrów FIR.

Pierwszy argumentem jest rząd filtru.
Drugim jest częstotliwość graniczna lub lista częstotliwości. Uwaga: nie podajemy ich w Hz, ale jako unormowane wartości od 0 do 1, gdzie 1 oznacza częstotliwość Nyquista.
Parametr rp jest wymagany dla filtru Czebyszewa I i eliptycznego, oznacza on maksymalny poziom zafalowań w paśmie przepustowym, w dB.
Parametr rs jest wymagany dla filtru Czebyszewa II i eliptycznego, oznacza minimalne tłumienie w paśmie zaporowym, w dB.
Parametr btype określa typ charakterystyki (w apostrofach), możemy podać: 'lowpass', 'highpass', 'bandpass', 'bandstop'.
Parametr ftype określa rodzaj filtru IIR (opisany poniżej), domyślnie: 'butter'.
Parametr output określa jaką formę wyniku chcemy uzyskać: współczynniki ('ba'), zera i bieguny ('zpk') czy sekcje drugiego rzędu ('sos').

Poniżej projektujemy dolnoprzepustowy filtr Butterwortha rzędu 8. W tybie ba funkcja zwraca dwie wartości: listę współczynników b i listę współczynników a. Obliczamy i wykreślamy odpowiedź impulsową filtru, przepuszczając przez niego impuls jednostkowy. Odpowiedź impulsowa jest nieskończona, pokazujemy pierwszych 150 próbek tej odpowiedzi. W wywołaniu funkcji lfilter podajemy najpierw współczynniki licznika, a następnie współczynniki mianownika. Do wygenerowania impulsu jednostkowego służy funkcja scipy.signal.unit_impulse, parametrem funkcji jest długość sygnału.

dp_b, dp_a = sig.iirfilter(8, 1000 / (fs / 2), btype='lowpass')
odp_imp = sig.lfilter(dp_b, dp_a, sig.unit_impulse(150))
plt.plot(odp_imp)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Odpowiedź impulsowa filtru IIR DP Butterwortha L=8')
plt.show()

Powyższy wykres pokazuje tylko fragment odpowiedzi impulsowej, która teoretycznie jest nieskończona. W stabilnych filtrach IIR, odpowiedź impulsowa ulega wygaszeniu. Widać, że dla indeksów próbek dużo większych niż rząd filtru, odpowiedź ipulsowa nadal jest niezerowa.

Parametr ftype określa rodzaj projektowanego filtru IIR. Cyfrowe filtry IIR są obliczane na podstawie projektu filtru analogowego. Istnieje kilka realizacji analogowego filtru, dlatego też mamy różne typy filtrów IIR, różniące się kształtem charakterystyki. Poniższy wykres ilustruje charakterystyki typów filtrów, które mamy do dyspozycji w funkcji iirfilter, dla filtru dolnoprzepustowego rzędu 8.

plt.axvline(1000, color='#c0c0c0')
b_butter, a_butter = sig.iirfilter(8, 1000 / (fs / 2), ftype='butter', btype='lowpass')
w, h_butter = sig.freqz(b_butter, a_butter)
f = w * fs / (2 * np.pi)
plt.plot(f, 20 * np.log10(np.abs(h_butter)), label='butter')
b_cheby1, a_cheby1 = sig.iirfilter(8, 1000 / (fs / 2), rp=1, ftype='cheby1', btype='lowpass')
_, h_cheby1 = sig.freqz(b_cheby1, a_cheby1)
plt.plot(f, 20 * np.log10(np.abs(h_cheby1)), label='cheby1')
b_cheby2, a_cheby2 = sig.iirfilter(8, 1000 / (fs / 2), rs=60, ftype='cheby2', btype='lowpass')
_, h_cheby2 = sig.freqz(b_cheby2, a_cheby2)
plt.plot(f, 20 * np.log10(np.abs(h_cheby2)), label='cheby2')
b_ellip, a_ellip = sig.iirfilter(8, 1000 / (fs / 2), rp=1, rs=60, ftype='ellip', btype='lowpass')
_, h_ellip = sig.freqz(b_ellip, a_ellip)
plt.plot(f, 20 * np.log10(np.abs(h_ellip)), label='ellip')
b_bessel, a_bessel = sig.iirfilter(8, 1000 / (fs / 2), ftype='bessel', btype='lowpass')
_, h_bessel = sig.freqz(b_bessel, a_bessel)
plt.plot(f, 20 * np.log10(np.abs(h_bessel)), label='bessel')
plt.legend()
plt.xlim(0, 5000)
plt.ylim(-120, 10)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('wzmocnienie [dB]')
plt.title('Porównanie typów charakterystyk filtrów IIR')
plt.show()

Filtr Butterwortha (ftype='butter') posiada maksymalnie płaską charakterystykę w paśmie przepustowym, nie zniekształca więc amplitud składowych w tym paśmie. Pasmo przejściowe jest natomiast bardzo szerokie. Dla filtru rzędu $N$, charakterystyka w paśmie zaporowym opada o $N \cdot 6$ dB na oktawę, gdzie oktawa oznacza podwojenie częstotliwości. Przykładowo, dla filtru DP rzędu 8 o częstotliwości granicznej 1 kHz, teoretyczna różnica wzmocnienia filtru pomiędzy 2 kHz a 4 kHz wynosi -48 dB. Aby uzyskać dobre tłuienie w paśmie zaporowym, musimy stosować filtry dużego rzędu.

Filtr Czebyszewa typu I (ftype='cheby1') dopuszcza zafalowania charakterystyki w paśmie przepustowym. Maksymalną wielkość tych zafalowań w dB podajemy jako parametr rp. Filtr zniekształca więc w pewnym stopniu amplitudy składowych sygnału w paśmie przepustowym. W zamian poprawia się tłumienie w paśmie zaporowym - szybkość opadania charaterystyki jest większa niż dla filtru Butterwortha. Większy rząd filtru poprawia tłumienie w pasmie zaporowym, ale zwiększa częstotliwość zafalowań w paśmie przepustowym.

Filtr Czebyszewa typu II (ftype='cheby2') nakłada ograniczenie na tłumienie w paśmie zaporowym — parametr rs. Tłumienie w paśmie zaporowym jest więc mniejsze niż dla poprzednich filtrów, otrzymujemy zafalowania w paśmie zaporowym, a ponadto pasmo przejściowe "zabiera" część pasma przepustowego. Charakterystyka w paśmie przepustowym jest w przybliżeniu płaska.

Filtr eliptyczny (ftype='ellip'), nazywany też filtrem Cauera albo filtrem Zołotarjewa, ma najwęższe pasmo przejściowe spośród wszyskich typów. Ceną za to jest jednak występowanie zafalowań w paśmie przepustowym oraz zaporowym, jak i zmniejszone tłumienie w paśmie zaporowym. Projekt filtru wymaga zatem podania obu parametrów: rp i rs. Ponieważ filtr eliptyczny wyrównuje wielkość zafalowań w paśmie przepustowym i zaporowym, jest też nazywany filtrem o równomiernych zafalowaniach (equiripple).

Filtr Bessela (ftype='bessel') ma "najgorszą" charakterystykę amplitudową spośród wszystkich typów, podobną do filtru Butterwortha, ale z jeszcze szerszym pasmem przejściowym. W technice analogowej filtry te są stosowane ze względu na charakterystykę fazową.

Aby porównać wpływ rzędu filtru na charakterystykę amplitudową, poniżej wykreślono charakterystykę filtru eliptycznego rzędu 4 i 8, oraz ten sam wykres ze zbliżeniem na pasmo przepustowe. Filtr rzędu 8 ma znacznie węższe pasmo przejściowe niż filtr rzędu 4, ale ma większą częstotliwość zafalowań w obu pasmach.

plt.axvline(1000, color='#c0c0c0')
b_ellip_n4, a_ellip_n4 = sig.iirfilter(4, 1000 / (fs / 2), rp=1, rs=60, ftype='ellip', btype='lowpass')
_, h_ellip_n4 = sig.freqz(b_ellip_n4, a_ellip_n4)
plt.plot(f, 20 * np.log10(np.abs(h_ellip_n4)), label='N=4')
plt.plot(f, 20 * np.log10(np.abs(h_ellip)), label='N=8')
plt.legend()
plt.xlim(0, 10000)
plt.ylim(-120, 10)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('wzmocnienie [dB]')
plt.title('Charakterystyki filtru eliptycznego DP')
plt.show()

plt.figure()
plt.plot(f, 20 * np.log10(np.abs(h_ellip_n4)), label='N=4')
plt.plot(f, 20 * np.log10(np.abs(h_ellip)), label='N=8')
plt.legend()
plt.xlim(0, 1200)
plt.ylim(-2, 1)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('wzmocnienie [dB]')
plt.title('Charakterystyki filtru eliptycznego DP (powiększenie)')
plt.show()

Którego typu filtru użyć? Wybór zależy od wymagań stawianych filtrowi. Na czym zależy nam najbardziej? Na płaskości charakterystyki w paśmie zaporowym, na konkretnej wartości tłumienia w paśmie zaporowym, na jak najwęższym paśmie przejściowym? W zależności od tego, dobierzemy właściwy typ. Ponadto musimy jeszcze rozpatrzyć charakterystykę fazową każdego typu.

Projektowanie środkowoprzepustowych i środkowozaporowych filtrów IIR ¶

Filtry środkowoprzepustowe i środkowozaporowe projektuje się tak jak filtr dolnoprzepustowy, z tym że należy podać listę dwóch częstotliwości granicznych, dolną i górną. Trzeba jednak zwrócić uwagę na pewien niuans. Popatrzmy na obliczone współczynniki dla filtru dolnoprzepustowego i środkowoprzepustowego, w obu przypadkach podając przy projektowaniu filtru rząd równy 8.

print('Filtr DP, N=8: liczba współczynników: {:2} wsp. b, {:2} wsp. a'.format(len(b_ellip), len(a_ellip)))
b_ellip_sp, a_ellip_sp = sig.iirfilter(8, (1000 / (fs / 2), 5000 / (fs / 2)), rp=1, rs=60, ftype='ellip', btype='bandpass')
print('Filtr SP, N=8: liczba współczynników: {:2} wsp. b, {:2} wsp. a'.format(len(b_ellip_sp), len(a_ellip_sp)))

Filtr DP, N=8: liczba współczynników:  9 wsp. b,  9 wsp. a
Filtr SP, N=8: liczba współczynników: 17 wsp. b, 17 wsp. a

Dla filtru dolnoprzepustowego, wszystko się zgadza, dostaliśmy filtr rzędu 8. Wyniki projektowania filtru środkowoprzepustowego mogą być jednak zaskoczeniem: dostaliśmy filtr rzędu 16, czyli inny niż ten, o który prosiliśmy. Aby zrozumieć ten efekt, trzeba wiedzieć w jaki sposób projektowane są filtry analogowe. Otóż projektowany jest filtr dolnoprzepustowy, a następnie z niego można uzyskać pozostałe typy charakterystyk, poprzez przekształcenie transmitancji. Dla filtrów SP i SZ, przekształcenie to powoduje dwukrotne zwiększenie rzędu filtru. Projektowanie filtrów cyfrowych opiera się na metodach projektowania filtrów analogowych, dlatego w filtrach IIR zachodzi ten sam efekt. Należy o tym pamietać, aby uniknąć zaskoczenia i podejrzewać funkcję projektującą filtr o błąd w algorytmie.

Charakterystyka fazowa i opóźnienia grupowego filtrów IIR ¶

Filtry o skończonej odpowiedzi impulsowej (FIR) projektowane do typowych zastosowań przepustowo-zaporowych charakteryzują się liniową charakterystyką fazową i stałym opóźnieniem grupowym. Dzięki temu opóźnienie wprowadzane przez filtr jest jednakowe dla wszystkich składowych częstotliwościowych i sygnał nie jest zniekształcany fazowo wskutek filtracji. Filtry IIR nie posiadają tej cechy. Popatrzmy na charakterystyki fazowe poszczególnych typów filtrów.

_, ax = plt.subplots(5, sharex=True, figsize=(8, 8))
ax[0].plot(f, np.degrees(np.angle(h_butter)), label='butter')
ax[1].plot(f, np.degrees(np.angle(h_cheby1)), label='cheby1')
ax[2].plot(f, np.degrees(np.angle(h_cheby2)), label='cheby2')
ax[3].plot(f, np.degrees(np.angle(h_ellip)), label='ellip')
ax[4].plot(f, np.degrees(np.angle(h_bessel)), label='bessel')
for a in ax:
    a.set_xlim(0, 5000)
    a.set_ylim(-180, 180)
    a.set_ylabel('faza')
    a.legend()
ax[-1].set_xlabel('częstotliwość {Hz}')
ax[0].set_title('Charakterystyki fazowe filtrów IIR różnego typu')
plt.tight_layout()
plt.show()

Widać, że charakterystyka fazowa nie jest odcinkami liniowa. Popatrzmy teraz na charakterystyki opóźnienia grupowego.

_, ax = plt.subplots(5, sharex=True, figsize=(8, 8))
np.warnings.filterwarnings('ignore')
ax[0].plot(f, sig.group_delay((b_butter, a_butter))[1], label='butter')
ax[1].plot(f, sig.group_delay((b_cheby1, a_cheby1))[1], label='cheby1')
ax[2].plot(f, sig.group_delay((b_cheby2, a_cheby2))[1], label='cheby2')
ax[2].set_ylim(bottom=-1)
ax[3].plot(f, sig.group_delay((b_ellip, a_ellip))[1], label='ellip')
ax[4].plot(f, sig.group_delay((b_bessel, a_bessel))[1], label='bessel')
for a in ax:
    a.set_xlim(0, 5000)
    #a.set_ylim(-180, 180)
    a.set_ylabel('op. grupowe')
    a.legend(loc='upper right')
ax[-1].set_xlabel('częstotliwość {Hz}')
ax[0].set_title('Charakterystyki opóźnienia grupowego filtrów IIR różnego typu')
plt.tight_layout()
plt.show()

Widać wyraźnie, że żaden z filtów nie zapewnia stałego opóźnienia grupowego. Ma to istotne znaczenie w praktyce, ponieważ w przypadku przetwarzania sygnałów szerokopasmowych (np. muzyki) przez filtr IIR, do sygnału zostaną wprowadzone zniekształcenia fazowe. Spośród wszstkich typów filtrów IIR, filtr Bessela jako jedyny posiada monotoniczną charakterystykę opóźnienia grupowego i takie jest uzasadnienie dla stosowania tego typu.

Zera i bieguny, stabilność filtru IIR ¶

Kolejnym problemem w projektowaniu filtrów IIR jest ich stabilność. Układem stabilnym nazywamy taki, w którym przy sygnale wejściowym przyjmującym wartości ograniczone do pewnego zakresu, sygnał wyjściowy również przyjmuje wartości z ograniczonego zakresu. Jest to tzw. stabilność BIBO (bounded input, bounded output). Pamiętajmy, że filtry IIR są układem ze sprzężeniem zwrotnym. Aby taki układ był stabilny, wzmocnienie w pętli sprzężenia zwrotnego musi być mniejsze od jedności. Przykładem stabilnego filtru IIR jest filtr dolnoprzepustowy pokazany na początku tego rozdziału. Patrząc na wykres odpowiedzi impulsowej tego filtru widzimy, że ma ona kształt gasnących oscylacji. To wygaszanie wynika ze wzmocnienia pętli sprzężenia zwrotnego mniejszego od jedności — po każdym "przejściu" przez filtr, wartości sygnału powracające na wejście stają się coraz mniejsze.

Co się stanie w przypadku gdy wzmocnienie w pętli sprzężenia zwrotnego będzie większe od jedności? Wartości z wyjścia układu zostają wzmocnione w pętli sprzężenia i dodane do wartości sygnału wejściowego. Po kolejnym przejściu przez filtr, zostają one ponownie wzmocnione i dodane do sygnału wejściowego. W ten sposób odpowiedź impulsowa zamiast wygasać, przyjmuje postać narastających oscylacji. Układ nie spełnia warunku BIBO i jest niestabilny. Oczywiście niestabilne filtry IIR nie mogą być stosowane w praktyce, ponieważ nie działają zgodnie z wymaganiami projektowymi.

Przypomnijmy wzór na transmitancję filtru IIR:

$$H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + ... + b_N z^{-N}}{1 + a_1 z^{-1} + a_2 z^{-2} + ... + a_N z^{-N}}$$

W liczniku i mianowniku transmitancji występują wielomiany stopnia $N$. Można również zapisać transmitancję w następującej formie:

$$H(z) = g \frac{\left(1-q_1 z^{-1}\right)\left(1-q_2 z^{-1}\right) ... \left(1-q_N z^{-1}\right)}{\left(1-p_1 z^{-1}\right)\left(1-p_2 z^{-1}\right) ... \left(1-p_N z^{-1}\right)}$$

Współczynniki $q$ są miejscami zerowymi licznika i są one nazywane zerami transmitancji. Współczynniki $p$ są miejscami zerowymi mianownika transmitancji i są one nazywane biegunami transmitancji (pole). Współczynnik $g$ to stałe wzmocnienie (gain) filtru.

Aby obliczyć zera i bieguny filtru, możemy w wywołaniu funkcji scipy.signal.iirfilter podać parametr output='zpk'. Wtedy wynikiem działania funkcji będa kolejno: tablica zer, tablica biegunów oraz wartość wzmocnienia. Jeżeli mamy już obliczone współczynniki (b, a), możemy uzyskać te same dane za pomocą funkcji scipy.signal.tf2zpk.

Zera i bieguny filtru wykreśla się na płaszczyźnie zespolonej: oś pozioma to część rzeczywista, oś pionowa to część urojona. Okrąg o promieniu równym 1 i środku w punkcie zerowym nazywa się okręgiem jednostkowym (unit circle). Zera zaznacza się na wykresie symbolem kółka, a bieguny – krzyżyka. Poniżej wykreślamy rozkład zer i biegunów dla filtru dolnoprzepustowego Butterwortha rzędu 8, o częstotliwości granicznej 1 kHz.

z, p, k = sig.tf2zpk(dp_b, dp_a)
_, ax = plt.subplots(subplot_kw={'aspect': 'equal'})
ax.add_patch(patches.Circle((0, 0), 1, fill=False, color="#808080"))
ax.axhline(0, color='#c0c0c0')
ax.axvline(0, color='#c0c0c0')
ax.plot(z.real, z.imag, 'o', color='#000080', ms=5)
ax.plot(p.real, p.imag, 'x', color='#800000', ms=5)
ax.set_xlim(-1.5, 1.5)
ax.set_ylim(-1.5, 1.5)
ax.set_xlabel('Re')
ax.set_ylabel('Im')
plt.title('Wykres zer i biegunów filtru DP Butterwortha 1 kHz, L=8')
plt.show()

Warunek stabliności jest następujący: moduł wszystkich biegunów musi być mniejszy od 1. Mówiąc bardziej obrazowo, wszystkie bieguny muszą leżeć wewnątrz okręgu jednostkowego.

Jeżeli przynajmniej jeden biegun transmitancji jest większy od 1, układ jest niestabilny. Jeżeli istnieją bieguny o module bliskim 1, układ taki nazywa się układem na granicy stabilności (marginally stable). Nie jest on ani stabilny, ani niestabilny, jego odpowiedź impulsowa ani nie wygasza się, ani nie narasta.

Mając obliczone bieguny, można sprawdzić stabilność układu w następujący sposób.

stabilny = np.all(np.abs(p) < 1)
print('Stabilny: {}, maksymalny moduł biegunów: {:.3f}'.format(stabilny, np.max(np.abs(p))))

Stabilny: True, maksymalny moduł biegunów: 0.975

Dla ilustracji problemu spróbujmy "popsuć" ten filtr czyniąc go niestabilnym i wykreślić jego odpowiedź impulsową.

p_niestabilny = np.copy(p)
p_niestabilny *= 1.1
b_niestabilny, a_niestabilny = sig.zpk2tf(z, p_niestabilny, k)
print('Stabilny: {}, maksymalny moduł biegunów: {:.3f}'.format(
    np.all(np.abs(p_niestabilny) < 1), np.max(np.abs(p_niestabilny))))
plt.plot(sig.lfilter(b_niestabilny, a_niestabilny, sig.unit_impulse(150)))
plt.title('Odpowiedź ipulsowa niestabilnego filtru IIR')
plt.show()

Stabilny: False, maksymalny moduł biegunów: 1.072

Widać, że odpowiedź impulsowa nie jest ograniczona i narasta coraz bardziej ze wzrostem indeksu próbki.

Trzeba szczególnie uważać przy implementacji filtrów w przypadku, gdy współczynniki filtru są poddawane kwantyzacji, np. na stałoprzecinkowych procesorach DSP. Jeżeli bieguny filtru leżą blisko okręgu jednostkowego, może zdarzyć się że po kwantyzacji "wyjdą" poza okrąg i filtr stanie się niestabilny.

Uwaga na koniec. Filtry o skończonej odpowiedzi impulsowej (FIR) mają mianownik transmitancji zawsze równy jedności. Wskutek tego, wszystkie bieguny filtru są położone w punkcie zerowym i w rezultacie filtry FIR zawsze są stabilne.

Sekcje drugiego rzędu ¶

Jest jeszcze jeden problem związany z rekursywnością filtrów IIR. Współczynniki filtru zapisywane są ze skończoną precyzją (kwantyzacja), przez co wyniki filtracji są obarczone błędem nazywanym szumem kwantyzacji (quantization noise). W filtrze rekursywnym, zaszumione wyniki filtracji są podawane pownownie na wejście filtru, przez co błędy kwantyzacji ulegają kumulacji i są coraz większe po każdym "przejściu przez filtr". Może to spowodować nie tylko powstanie niedokładnych wyników, ale również utratę stabilności filtru. Efekt ten jest zatem tym bardziej dotkliwy, im większy jest rząd filtru. Dla filtrów IIR rzędu 1 lub 2, efekt ten jest mniej istotny. W praktyce potrzebujemy jednak filtrów wyższego rzędu.

Przypomnijmy wzór na transmitancję filtru IIR:

$$H(z) = g \frac{\left(1-q_1 z^{-1}\right)\left(1-q_2 z^{-1}\right) ... \left(1-q_N z^{-1}\right)}{\left(1-p_1 z^{-1}\right)\left(1-p_2 z^{-1}\right) ... \left(1-p_N z^{-1}\right)}$$

Ten sam wzór możemy też zapisać w następującej równoważnej formie:

$$H(z) = g_1 \frac{\left(1-q_1 z^{-1}\right)\left(1-q_2 z^{-1}\right)}{\left(1-p_1 z^{-1}\right)\left(1-p_2 z^{-1}\right)} \cdot g_2 \frac{\left(1-q_3 z^{-1}\right)\left(1-q_4 z^{-1}\right)}{\left(1-p_3 z^{-1}\right)\left(1-p_4 z^{-1}\right)} \cdot ... \cdot g_{N/2} \frac{\left(1-q_{N-1} z^{-1}\right)\left(1-q_N z^{-1}\right)}{\left(1-p_{N-1} z^{-1}\right)\left(1-p_N z^{-1}\right)}$$

gdzie iloczyn współczynników $g_i$ jest równy $g$. Taki zapis transmitancji reprezentuje kaskadowe połączenie sekcji drugiego rzędu (SOS, second order section), nazywanych również sekcjami dwukadratowymi (biquad). Jak wspomniano powyżej, dla filtrów drugiego rzędu wpływ błędów kwantyzacji jest mniejszy niż dla wyższych rzędów, zatem np. kaskada czterech filtrów drugiego rzędu jest równoważna filtrowi ósmego rzędu, ale cechuje się mniejszym wpływem błędów kwantyzacji. Z tego powodu zaleca się stosowanie tek struktury, szczególnie w przypadkach, w których występują duże błędy kwantyzacji, np. w implementacji stałoprzecinkowej. Każdy filtr IIR parzystego rzędu można przedstawić jako kaskadę sekcji drugiego rzędu. Dla nieparzystych rzędów filtru trzeba zastosować dodatkowo jedną sekcję pierwszego rzędu.

Aby zaprojektować filtr IIR złożony z sekcji drugiego rzędu należy w wywołaniu funkcji scipy.signal.iirfilter podać parametr output='sos'. Dla filtru dolnoprzepustowego Butterwortha rzędu 8 o czestotliwości granicznej 1 kHz otrzymujemy następujący wynik:

dp_sos = sig.iirfilter(8, 1000 / (fs / 2), btype='lowpass', ftype='butter', output='sos')
print(dp_sos)

[[  2.43444902e-10   4.86889804e-10   2.43444902e-10   1.00000000e+00
   -1.75785265e+00   7.73021088e-01]
 [  1.00000000e+00   2.00000000e+00   1.00000000e+00   1.00000000e+00
   -1.78875835e+00   8.04193476e-01]
 [  1.00000000e+00   2.00000000e+00   1.00000000e+00   1.00000000e+00
   -1.84881984e+00   8.64773233e-01]
 [  1.00000000e+00   2.00000000e+00   1.00000000e+00   1.00000000e+00
   -1.93365048e+00   9.50335873e-01]]

Wynikiem działania funkcji iirfilter jest w tym przypadku macierz o wymiarach $N \times 6$. Każdy rząd macierzy to jedna sekcja drugiego rzędu, współczynniki dla każdej sekcji są podawane w kolejności: $[b_0, b_1, b_2, 1, a_1, a_2]$.

Dla filtru zaprojektowanego w formie SOS musimy użyć specjalnych funkcji. W celu przefiltrowania sygnału musimy użyć funkcji scipy.signal.sosfilt, podając macierz SOS jako pierwszy parametr i filtrowany sygnał jako drugi parametr.

szum = np.random.random(2048) * 2 - 1
przefiltrowany = sig.sosfilt(dp_sos, szum)
h_szum = np.fft.rfft(przefiltrowany * np.hamming(2048))
f_szum = np.fft.rfftfreq(2048, 1 / fs)
plt.plot(f_szum, 20 * np.log10(np.abs(h_szum / 1024)))
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Wynik filtracji szumu filtrem DP za pomocą sosfilt')
plt.xlim(0, 5000)
plt.show()

Do obliczenia odpowiedzi częstotliwościowej służy funkcja scipy.signal.sosfreqz. Warunki początkowe dla filtru możemy obliczyć funkcją scipy.signal.sosfilt_zi. Możemy również zamienić formę SOS na formę pełną (b, a) za pomocą funkcji scipy.signal.sos2tf. Konwersja w odwrotnym kierunku jest możliwa za pomocą funkcji scipy.signal.tf2sos, ale nie zaleca się tego robić, ze względu na dokładność.

Poniższy przykład (zaczęrpnięty z dokumentacji modułu scipy.signal) pokazuje porównanie charakterystyk filtrów zaprojektowanych metodą klasyczną i jako sekcje dwukwadratowe. Jest to filtr środkowoprzepustowy, eliptyczny, rzędu 15, o względnych częstotliwościach granicznych 0,2 i 0,4, zafalowaniach 0,5 dB w paśmie przepustowym i tłumieniu -60 dB w paśmie zaporowym.

# filtr projektowany jako (b, a)
b_pp, a_pp = sig.iirfilter(15, (0.2, 0.4), rp=0.5, rs=60, ftype='ellip', btype='bandpass')
w, h_ba = sig.freqz(b_pp, a_pp)
f = w * fs / (2 * np.pi)

# filtr projektowany jako SOS
sos_pp = sig.iirfilter(15, (0.2, 0.4), rp=0.5, rs=60, ftype='ellip', btype='bandpass', output='sos')
_, h_sos = sig.sosfreqz(sos_pp)

plt.plot(f, 20 * np.log10(np.abs(h_ba)), label='ba')
plt.plot(f, 20 * np.log10(np.abs(h_sos)), label='sos')
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Charakterystyki filtrów SP N=15')
plt.ylim(bottom=-90)
plt.legend()
plt.show()

Na powyższym wykresie widać wyraźnie błędy powstające w paśmie przepustowym dla filtru projektowanego metodą ba. Ten sam filtr projektowany metodą sos jest wolny od tych błędów. Ale to nie wszystko. Popatrzmy na odpowiedzi impulsowe obu filtrów.

impuls = sig.unit_impulse(150)
impuls_ba = sig.lfilter(b_pp, a_pp, impuls)
impuls_sos = sig.sosfilt(sos_pp, impuls)
plt.plot(impuls_ba, label='ab')
plt.plot(impuls_sos, label='sos')
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Odpowiedzi impulsowe filtrów SP L=15')
plt.ylim(-0.25, 0.25)
plt.legend()
plt.show()

Co się stało? Filtr rzędu 15 zaprojektowany jako pojedyncza sekcja jest niestabilny! Ten sam filtr zaprojektowany jako sekcje drugiego rzędu zachowuje stabilność, chociaż jest zaprojektowany na granicy stabilności.

_, p_ba, _ = sig.tf2zpk(b_pp, a_pp)
print('Filtr ba  - stabilny: {}, maksymalny moduł biegunów: {:.3f}'.format(
    np.all(np.abs(p_ba) < 1), np.max(np.abs(p_ba))))
_, p_sos, _ = sig.sos2zpk(sos_pp)
print('Filtr sos - stabilny: {},  maksymalny moduł biegunów: {:.3f}'.format(
    np.all(np.abs(p_sos) < 1), np.max(np.abs(p_sos))))

Filtr ba  - stabilny: False, maksymalny moduł biegunów: 1.048
Filtr sos - stabilny: True,  maksymalny moduł biegunów: 1.000

Na skutek błędów numerycznych powstających przy filtracji, które są dodawane do sygnału wejściowego wielokrotnie wskutek rekursji, bieguny zostały "wypchnięte" poza okrąg jednostkowy. Filtr stał się więc niestabilny.

Z powyższego przykładu wynika wniosek, że nie powinniśmy stosować filtrów IIR dużych rzędów w formie pojedynczej sekcji, ponieważ narażamy się na ryzyko niestabilności fitru i na zniekształcenia charakterystyki częstotliwościowej. Forma sekcji drugiego rzędu jest znacznie mniej wrażliwa na niedokładności i dlatego powinniśmy ją stosować.

Filtracja zerofazowa ¶

Na wcześniejszych przykładach pokazano, że filtry IIR projektowane w opisany wyżej sposób mają nieliniową charakterystykę fazową i zmienne opóźnienie grupowe. Jest to istotna wada przy przetwarzaniu sygnałów szerokopasmowych. Istnieje jednak możliwość zaprojektowania filtrów IIR nie zniekształcających fazy. Będą to jednak filtry nieprzyczynowe. Nie będzie można ich zastosować w przypadku, gdy próbki sygnału są przetwarzane "na żywo". Ale jeżeli dysponujemy całym sygnałem, możemy zastosować opisywaną metodę.

Zasada filtracji zerofazowej (zero-phase filtering) jest następująca. Najpierw filtrujemy sygnał normalnie. Następnie odwracamy wynik filtracji w czasie i ponownie przepuszczamy przez filtr. Jest to tzw. filtr "w przód-wstecz" (forward-backward filter) Wynikowy sygnał nie będzie zniekształcony fazowo względem oryginalnego sygnału.

Dlaczego to działa? Przy pierwszym przetwarzaniu, do sygnału dodawane są zniekształcenia fazowe spowodowane przez nieliniową charakterystykę filtru. Następnie odwracamy sygnał w czasie — operacja ta powoduje odwrócenie fazy oryginalnego sygnału względem zera (zmianę znaku). Po ponownym przejściu sygnału przez filtr, dodane zniekształcenia mają przeciwny znak do tych dodanych do sygnału w pierwszej fazie. W rezultacie mamy filtr, którego charakerystyka amplitudowa jest równa kwadratowi charakterystyki oryginalnego filtru, natomiast charakterystyka fazowa jest zerowa – stąd nazwa filtru.

W module scipy.signal mamy do dyspozycji funkcje o nazwach scipy.signal.filtfilt dla filtrów jednosekcyjnych oraz scipy.signal.sosfiltfilt dla filtrów złożonych z sekcji dwukwadratowych. Dziwne nazwy miały prawdopodobnie zaznaczać fakt dwukrotnej filtracji sygnału.

Zobaczmy wynik filtracji dla normalnego eliptycznego filtru dolnoprzepustowego rzędu 8, o częstotliwości granicznej 1 kHz. Filtrowanym sygnałem jest suma sinusów o częstotliwościach 200, 500 i 800 Hz, a więc leżących w paśmie przepustowym filtru. Wykres pokazuje charakterystykę fazową sygnału przed i po filtracji.

n = np.arange(2048)
trzysin = (np.sin(2 * np.pi * n * 200 / fs) +
           np.sin(2 * np.pi * n * 500 / fs) +
           np.sin(2 * np.pi * n * 800 / fs))
b, a = sig.iirfilter(8, 1000 / (fs / 2), rp=0.5, rs=60, btype='lowpass', ftype='bessel')
trzysin_filt = sig.lfilter(b, a, trzysin)
f_fft = np.fft.rfftfreq(2048, 1 / fs)
plt.plot(f_fft, np.degrees(np.angle(np.fft.rfft(trzysin))), label='oryginalny')
plt.plot(f_fft, np.degrees(np.angle(np.fft.rfft(trzysin_filt))), label='przefiltrowany')
plt.xlim(0, 2000)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('faza [st]')
plt.legend(loc='lower right')
plt.title('Charakterystyki fazowe sygnału przed i po filtracji IIR')
plt.show()

Wyrażnie widać że filtr zniekształcił fazę sygnału. Teraz wykonamy filtrację zerofazową tego samego sygnału.

trzysin_filtfilt = sig.filtfilt(b, a, trzysin)
plt.plot(f_fft, np.degrees(np.angle(np.fft.rfft(trzysin))), label='oryginalny')
plt.plot(f_fft, np.degrees(np.angle(np.fft.rfft(trzysin_filtfilt))), label='przefiltrowany')
plt.xlim(0, 2000)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('faza [st]')
plt.legend(loc='lower right')
plt.title('Charakterystyki fazowe sygnału przed i po filtracji zerofazowej IIR')
plt.show()

Pewne różnice charakterystyk fazowych są zauważalne, ale obie charakterystyki są znacznie bardziej zbliżone do siebie niż w poprzednim przypadku.

Funkcje filtfilt i sosfiltfilt posiadają dodatkowe parametry, które umożliwiają dostrojenie działania filtru. Funkcje te stosują padding, czyli dodają zera przed i po filtrowanym sygnale aby zredukować błędy wynikające z niezerowych warunków początkowych. Parametr padtype umożliwia ustalenie metody paddingu, a parametr padsize – liczby dodawanych elementów. Alternatywną metodą jest metoda Gustaffsona, włączana parametrem method='gust', wtedy parametr irlen określa długość odpowiedzi impulsowej filtru. Szczegóły można znaleźć w dokumentacji.

Zaleca się stosowanie filtracji zerofazowej zawsze w przypadkach gdy mamy do dyspozycji kompletny sygnał i chcemy go przetworzyć za pomocą filtru IIR. Redukujemy w ten sposób zniekształcenia fazowe wprowadzane przez filtr. Możemy również zastosować metodę zerofazową do filtru FIR. Filtry te nie zniekształcają fazy, ale filtracja zerofazowa powoduje że filtr nie wprowadza opóźnienia, sygnał oryginalny i przefiltrowany pozostają zsynchronizowane.

Alternatywna metoda projektowania filtrów IIR ¶

W module scipy.signal jest także dostępna funkcja iirdesign. Różni się ona od iirfilter tym, że umożliwia wyspecyfikowanie szerokości pasma przejściowego filtru, a rząd filtru jest dobierany do tej specyfikacji. Pierwsze dwa parametry, wp i ws, wyznaczają krańce pasma przepustowego i zaporowego. Dla filtrów pasmowo-przepustowych i pasmowo-zaporowych należy podać po dwie wartości graniczne. Podobnie jak dla iirfilter, wartości częstotliwości są unormowane od 0 do 1, gdzie 1 jest równe częstotliwości Nyquista. Parametry gpass i gstop wyznaczają odpowiednie: maksymalny pozio zafalowań w paśmie przepustowym i minimalne tłumienie w paśmie zaporowym. Parametry ftype i output mają takie samo znaczenie jak poprzednio.

Zaprojektujmy filtr pasmowo-przepustowy, który ma pasmo przepustowe od 2 do 4 kHz z pasmem przejściowym po 500 Hz z każdej strony, z poziomem zafalowań w paśmie przepustowym 0,5 dB i minimalnym tłumieniem 60 dB.

fs2 = fs / 2
sos_pp2 = sig.iirdesign((2000 / fs2, 4000 / fs2), (1500 / fs2, 4500 / fs2), 0.5, 60, ftype='ellip', output='sos')
print('Rząd filtru:', sos_pp2.shape[0] * 2)
w, h_pp2 = sig.sosfreqz(sos_pp2)
f = w * fs / (2 * np.pi)
plt.plot(f, 20 * np.log10(np.abs(h_pp2)))
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.xlim(0, 10000)
plt.ylim(bottom=-90)
plt.title('Charakterystyka filtru SP zaprojektowana funkcją iirdesign')
plt.show()

Rząd filtru: 14

Filtry IIR a FIR ¶

Przedstawiliśmy filtry o skończonej i nieskończonej odpowiedzi impulsowej. Podsumujmy najważniejsze cechy obu typów.

Cecha	FIR	IIR
Zawsze są stabilne	Tak	Nie
Mogą być liniowofazowe	Tak	Nie
Wpływ szumu kwantyzacji	Mały	Duży
Wymagany rząd filtru	Większy	Mniejszy

Ostatnia z wymienionych cech wyjaśnia dlaczego w ogóle używamy filtrów IIR, pomimo licznych wad tego typu w porównaniu z filtrami FIR. Otóż filtry IIR zazwyczaj wymagają wykonania znacznie mniejszej ilości obliczeń niż filtr FIR o zbliżonej charakterystyce. Wynika to z tego, że aby uzyskać pożądany kształt charakterystyki, w przypadku filtru IIR wystarczy znacznie mniejszy rząd. Popatrzmy na porównanie dwóch filtrów o zbliżonym kształcie charakterystyki: filtr FIR projektowany oknem Kaisera oraz filtr eliptyczny IIR.

# FIR
N_fir, beta = sig.kaiserord(50, 500 / fs2)
print('Rząd filtru FIR:', N_fir - 1)
dp_fir = sig.firwin(N_fir, 1000, window=('kaiser', beta), fs=fs)
w, h_fir = sig.freqz(dp_fir)
plt.plot(w * fs / (2 *np.pi), 20 * np.log10(np.abs(h_fir)), label='FIR')

# IIR
dp_iir = sig.iirdesign(1000 / fs2, 1500 / fs2, 0.5, 60, ftype='ellip', output='sos')
print('Rząd filtru IIR:', dp_iir.shape[0] * 2)
w, h_iir = sig.sosfreqz(dp_iir)
plt.plot(w * fs / (2 *np.pi), 20 * np.log10(np.abs(h_iir)), label='IIR')

plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Porównanie charakterystyk filtrów FIR i IIR')
plt.xlim(0, 5000)
plt.ylim(bottom=-90)
plt.legend()
plt.show()

Rząd filtru FIR: 282
Rząd filtru IIR: 6

Różnica rzędu obu filtrów jest bardzo wyrażna. Ten aspekt ma bardzo duże znaczenie przy implementacji filtrów np.na procesorze sygnałowym, gdzie liczy się czas obliczeń - filtry IIR pozwalają uzyskać ten sam efekt znacznie mniejszym kosztem.

Podsumowanie ¶

Filtry IIR (o nieskończonej odpowiedzi impulsowej) są cyfrowym odpowiednikiem analogowych filtrów.
Filtr IIR jest rekursywny. W filtracji biorą udział także poprzednie wyniki filtracji.
Transmitancja filtru jest opisana ułamkiem, ze współczynnikami b w liczniku i a w mianowniku.
Projektowanie filtrów IIR opiera się na metodach projektowania filtrów analogowych i stosowaniu przekształceń.
Przy projektowaniu filtru należy wybrać typ filtru. Wybór zależy od wymagań, typy filtru różnią się szerokością pasma przejściowego, wielkością zafalowań w paśmie przepustowym i wielkością tłumienia w paśmie zaporowym.
Transmitancję filtru IIR można zapisać w formie zer, biegunów i wzmocnienia.
Filtr IIR nadający się do praktycznego zastosowania musi być stabilny – moduł wszystkich biegunów musi być mniejszy od 1.
Filtry IIR nie mają liniowej charakterystyki fazowej, opóźnienie grupowe zależy od częstotliwości. Filtry IIR mogą zniekształcać fazowo sygnały szerokopasmowe.
Filtracja zerofazowa nie wprowadza zniekształceń fazowych, ale filtr jest nieprzyczynowy – musimy mieć pełny sygnał do dyspozycji.
Błędy numeryczne kumulują się przy filtracji IIR, zniekształcają one charakterystykę filtru i mogą spowodować jego niestabilność.
Implementacja filtrów IIR za pomocą sekcji drugiego rzędu zmniejsza wpływ błędów na filtr i jest zalecana, szczególnie dla większych rzędów filtru.
Filtry IIR potrzebują znacznie mniejszego rzędu niż filtry FIR aby uzyskać podobną skuteczność filtracji. Są więc bardziej wydajne, wymagają mniejszej ilości obliczeń.