Filtry

AS/ Funkcja przejścia i filtry

Cyfrowe filtry liniowe niezmiennicze w czasie [math]x[n] \longrightarrow \boxed{LTI} \longrightarrow y[n][/math] opisuje liniowe równanie różnicowe o stałych współczynnikach

[math] \displaystyle y[n] = \sum_{k=1}^K a_k y[n-k] + \sum_{l=0}^L b_l x[n-l] [/math]

gdzie [math]a_i[/math] i [math]b_i[/math] to współczynniki, [math]x[n][/math] to sygnał wejściowy, a [math]y[n][/math] — wyjście; w ogólniejszej postaci można je zapisać jako

[math] \displaystyle \sum_{k=0}^K a_k y[n-k] = \sum_{l=0}^L b_l x[n-l] [/math]

Funkcja przejścia (transfer function)

Stosując transformatę [math]Z[/math] możemy powyższe równanie z dziedziny czasu przenieść do dziedziny częstości:

[math] \displaystyle \mathcal{Z}\left\{\sum_{k=0}^K a_k y[n-k] \right\} = \mathcal{Z}\left\{ \sum_{l=0}^L b_l x[n-l] \right\} [/math]

[math] \displaystyle \sum_{k=0}^K a_k \mathcal{Z}\left\{ y[n-k]\right\} = \sum_{l=0}^L b_l \mathcal{Z} \left\{x[n-l]\right\} [/math]

[math] \displaystyle \sum_{k=0}^K a_k z^{-k} Y(z) = \sum_{l=0}^L b_l z^{-l} X(z) [/math]

[math] \displaystyle Y(z) \sum_{k=0}^K a_k z^{-k} = X(z) \sum_{l=0}^L b_l z^{-l} [/math]

Dla systemu przyczynowego dostajemy:

[math] \displaystyle \frac{Y(z)}{X(z)} \equiv H(z) = \frac{\sum_{l=0}^L b_l z^{-l}}{\sum_{k=0}^K a_k z^{-k}} [/math]

Filtry

Funkcja przejścia [math]H(z)[/math] pozwala spójnie przedstawić działanie filtra LTI na sygnał [math]x[/math] w przestrzeni transformaty [math]\mathcal{Z}[/math]:

[math] \displaystyle Y(z)=H(z)X(z)=\frac{b_0 + b_1 z^{-1}+\dots +b_{L} z^{-L}}{a_0+a_1 z^{-1}+\dots +a_{K} z^{-K}}X(z) [/math]

gdzie [math]X(z)[/math] to transformata [math]\mathcal{Z}[/math] filtrowanego sygnału (wejścia [math]x[/math]), [math]Y(z)[/math] — wyjścia, a [math]H(z)[/math] to wprowadzona powyżej funkcja przejścia filtra. Filtrowanie w przestrzeni transformaty [math]\mathcal{Z}[/math] odpowiada przemnożeniu transformaty sygnału przez transformatę funkcji przejścia filtru. Zobaczmy, co się stanie z widmem sygnału po przefiltrowaniu. Weźmy w tym celu [math]z = e^{i\omega}[/math]. Wówczas transmitancja jest funkcją częstości [math]\omega[/math].

[math] \displaystyle Y(e^{i\omega})=H(e^{i\omega})X(e^{i\omega})=\frac{b_0+b_1 e^{-i\omega}+\dots +b_L e^{-i L\omega}}{a_0+a_1 e^{-i\omega}+\dots +a_K e^{-i K\omega}}X(e^{i\omega}) [/math]

Działanie funkcji przejścia na składową o częstości [math]\omega_k[/math] odpowiada wymnożeniu jej przez liczbę zespoloną [math]H(e^{i\omega_k})[/math], którą, jak każdą liczbę zespoloną, można przedstawić we współrzędnych biegunowych jako [math]|A_k| e^{i \phi_k}[/math]:

[math]Y(\omega_k) = |A_k| e^{i \phi_k} |X_k| e^{i \theta_k} = |A_k| |X_k| e^{i ( \phi_k +\theta_k)} [/math]

W wyniku filtrowania sinusoidalna składowa sygnału o danej częstości może zmienić amplitudę i fazę, ale nie zmienia częstości.

Finite Impulse Response (FIR) — filtr o skończonej odpowiedzi impulsowej

Jeśli w równaniu

[math] \displaystyle \sum_{k=0}^K a_k y[n-k] = \sum_{l=0}^L b_l x[n-l] [/math]

położymy [math]a_i = 0[/math] poza [math]a_0=1[/math], dostaniemy

[math] \displaystyle y[n] = \sum_{l=0}^L b_l x[n-l] [/math]

W funkcji przejścia mianownik będzie stały i dostaniemy

[math] \displaystyle Y[z]=H[z]X[z]=\left(b_0+b_1 z^{-1}+\dots +b_L z^{-L}\right) X[z] [/math]

a w dziedzienie czasu, z pełnego równania

[math] \displaystyle y[n] = b_0 x[n]+ b_1 x[n-1] + \dots + b_L x[n-L] \,\,\, - \,\,\, a_1 y[n-1] - \dots - a_K y[n-K] [/math]

pozostaje

[math] \displaystyle y(n) = b_0 *x[n] + b_1 *x[n-1] + \dots + b_L *x[n-L] [/math]

Jeśli współczynniki [math]b_i[/math] zapiszemy jako sekwencję [math]b[i][/math], działanie takiego filtra LTI na sygnał [math]x[n][/math] w dziedzinie czasu będze wyglądała tak:

[math] \displaystyle y(n) = b[0]*x[n] + b[1]*x[n-1] + \dots + b[L]*x[n-L] = b[n]*x[n] [/math]

Taki filtr nazywamy filtrem o skończonej odpowiedzi impulsowej (Finite Impulse Response, FIR), bo odpowiedź na impuls kończy się po [math]L[/math] próbkach. Inna nazwa to średnia biegnąca (Moving Average, MA).

Infinite Impulse Response (IIR) — filtr o nieskończonej odpowiedzi impulsowej

Jeśli w równaniu

[math] \displaystyle \sum_{k=0}^K a_k y[n-k] = \sum_{l=0}^L b_l x[n-l] [/math]

położymy [math]b_i = 0[/math] poza [math]b_0=1[/math], dostaniemy

[math] \displaystyle \sum_{k=0}^K a_k y[n-k] = x[n] [/math]

w funkcji przejścia licznik będzie stały

[math] \displaystyle Y[z]=H[z]X[z]=\frac{1}{a_0+a_1 z^{-1}+\dots +a_{K} z^{-K}}X[z] [/math]

Operacja splotu działa tu na sekwencji wyjściowej:

[math] y[n] = x[n] - a[1]*y[n-1] - \dots - a[n_a]*y[n-K] [/math]

dlatego taki filtr nazwać można filtrem rekursywnym, autoregresyjnym (AR), lub filtrem o nieskończonej odpowiedzi impulsowej (Infinite Impulse Response IIR), bo potencjalnie raz wzbudzony (sekwencją jednostkową) może dowolnie długo produkować niezerowe wyjście. Faza filtrowanego sygnału zaburzana jest nieliniowo (nonlinear phase filter).

Liniowe i nieliniowe opóźnienie fazy, opóźnienie grupowe

Jak ustaliliśmy powyżej, filtry liniowe odpowiadają mnożeniu poszczególnych składowych sygnału wejściowego [math]x[n][/math] przez liczbę zespoloną, mającą charakterystyczną dla danego filtra amplitudę i fazę, które dla każdej częstości [math]\omega[/math] można wyliczyć z funkcji przenoszenia [math]H(\omega)[/math].

Różna dla różnych częstości zmiana amplitudy jest kwintesencją i celem całego procesu filtrowania, w którym z sygnału usuwamy (w praktyce raczej osłabiamy) niechciane częstości. Jednak wywołana tym procesem zmiana faz dla różnych częstości nie jest już tak oczywista. Weźmy sygnał złożony z dwóch oscylacji, o częstościach [math]\omega[/math] i [math]2 *\omega[/math]:

[math]x(t) = \sin(\omega t) + \sin(2 * \omega t)[/math]

Suma tych dwóch składowych (niebieskie krzywe) na poniższym rysunku przedstawiona jest kolorem czerwonym:

Przesuńmy teraz fazę obydwu sygnałów składowych o tę samą wartość, na przykład [math]\pi/2[/math], i ponownie dodajmy do siebie. Odpowiada to sytuacji, kiedy po przejściu przez filtr sygnału [math]x(t) = \sin(\omega t) + \sin(2 * \omega t)[/math] dostajemy [math]y(t) = \sin(\omega t + \pi/2) + \sin(2 * \omega t + \pi/2)[/math]

Po przejściu przez system, opóźniający fazę [math]\phi[/math] jednolicie (o [math]\Delta\phi=\pi/{2}[/math]) dla wszystkich częstości, sygnał wyjściowy wygląda zupełnie inaczej, niż sygnał wejściowy. Przed ponownym "złożeniem" każdy z sygnałów składowych został przesunięty o inny odcinek czasu, który dla częstości [math]\omega[/math] wynosi [math]\phi/{\omega}[/math]:

[math]\displaystyle \sin(\omega t + \phi) = sin\left(\,\,\omega (t + \phi/{\omega}) \,\, \right) = \sin\left(\,\, \omega (t + \Delta t)\,\, \right) [/math]

Formalnie, opóźnienie fazy [math]\phi[/math] dla danej częstości [math]\omega[/math] określamy mianem opóźnienia fazowego, i definiujemy jako [math]-\phi/\omega[/math].

Przesunięcie składowej sygnału w czasie zwiemy opóźnieniem grupowym, i definiujemy jako [math]-\frac{d \phi(\omega)}{d \omega}[/math].

Żeby otrzymać sygnał, którego wszystkie składowe będą przesunięte o ten sam odcinek czasu [math]\Delta t[/math], czyli będą miały stałe przesunięcie grupowe, trzeba fazę przesuwać liniowo w zależności od częstości: [math]\Delta\phi(\omega) = \omega\Delta t[/math], wtedy opóźnienie grupowe [math]\frac{d \phi(\omega)}{d \omega} = \frac{d (\omega\Delta t)}{d \omega} = \Delta t[/math] będzie stałe.

W każdym innym przypadku sygnał na wyjściu nie będzie odtwarzał własności fazowych (czyli "kształtu", obserwowanego w dziedzinie czasu) sygnału oryginalnego, nawet w zakresie częstości nie zmienianych w procesie filtrowania — zgadzać się będzie tylko widmo w paśmie przenoszenia.

kod (uproszczony)

import matplotlib.pyplot as plt
import numpy as np

t = np.arange(0.0, 3, 0.01)
s1 = np.sin(2 * np.pi * t)
s2= np.sin(2 * np.pi * t + np.pi/2)
s3= np.sin(2*2 * np.pi * t)
s4= np.sin(2*2 * np.pi * t + np.pi/2)
s5= np.sin(2*2 * np.pi * t + 2*np.pi/2)

fig1, ax1 = plt.subplots(3, 1)
ax1[0].plot(t, s1)
ax1[1].plot(t, s3)
ax1[2].plot(t, s1+s3, 'r')
ax1[0].set_title('sin(wt)')
ax1[1].set_title('sin(2wt)')
ax1[2].set_title('sin(wt) + sin(2wt)')

fig2, ax2 = plt.subplots(3, 1)
ax2[0].plot(t, s2)
ax2[1].plot(t, s4)
ax2[2].plot(t, s2+s4, 'r')
ax2[0].set_title('sin(wt+pi/2)')
ax2[1].set_title('sin(2wt+pi/2)')
ax2[2].set_title('sin(wt+pi/2) + sin(2wt+pi/2)')

fig3, ax3 = plt.subplots(3, 1)
ax3[0].plot(t, s2)
ax3[1].plot(t, s5)
ax3[2].plot(t, s2+s5, 'r')
ax3[0].set_title('sin(w(t+pi/2))')
ax3[1].set_title('sin(2wt+2*pi/2)')
ax3[2].set_title('sin(wt+pi/2) + sin(2wt+2*pi/2)')

plt.show()

kod generujący obrazki

import matplotlib.pyplot as plt
import numpy as np

t = np.arange(0.0, 3, 0.01)
s1 = np.sin(2 * np.pi * t)
s2 = np.sin(2 * np.pi * t + np.pi / 2)
s3 = np.sin(2 * 2 * np.pi * t)
s4 = np.sin(2 * 2 * np.pi * t + np.pi / 2)
s5 = np.sin(2 * 2 * np.pi * t + 2 * np.pi / 2)

def fix_axes():
    for i in range(0, 3):
        ax[i].set_yticks([0], [])
        ax[i].set_xlim([0, 3])
        ax[i].grid(True)
        ax[i].set_frame_on(False)
        ax[i].tick_params(length=0)
        ax[i].set_ylim(-2, 2)

fig1, ax = plt.subplots(3, 1, dpi=200)
ax[0].plot(t, s1)
ax[1].plot(t, s3)
ax[2].plot(t, s1 + s3, 'r')
ax[0].set_title('sin( $\omega$t )', loc='left')
ax[1].set_title('sin( 2*$\omega$t )', loc='left')
ax[2].set_title('sin( $\omega$t ) + sin( 2*$\omega$t )', loc='left')
fix_axes()
for i in range(0, 3):
    ax[i].set_xticks(np.arange(.25, 3, 1), [])
plt.show()

fig2, ax = plt.subplots(3, 1, dpi=200)
ax[0].plot(t, s2)
ax[1].plot(t, s4)
ax[2].plot(t, s2 + s4, 'r')
ax[0].set_title('sin( $\omega$t + $\pi$/2 )', loc='left')
ax[1].set_title('sin( 2*$\omega$t + $\pi/2$ )', loc='left')
ax[2].set_title('sin( $\omega$t + $\pi$/2 ) + sin( 2*$\omega$t + $\pi$/2 )', loc='left')
fix_axes()
for i in range(0, 3):
    ax[i].set_xticks(np.arange(0, 3, 1), [])
plt.show()

fig3, ax = plt.subplots(3, 1, dpi=200)
ax[0].plot(t, s2)
ax[1].plot(t, s5)
ax[2].plot(t, s2 + s5, 'r')
ax[0].set_title('sin( $\omega$(t + $\pi$/2) )', loc='left')
ax[1].set_title('sin( 2*$\omega$t + 2*$\pi$/2 )', loc='left')
ax[2].set_title('sin( $\omega$t + $\pi$/2 ) + sin( 2*$\omega$t + 2*$\pi$/2 )', loc='left')
fix_axes()
for i in range(0, 3):
    ax[i].set_xticks(np.arange(0, 3, 1), [])
plt.show()