In [1]:
%matplotlib inline
import numpy as np
import scipy.signal as sig
import scipy.linalg as linalg
import matplotlib.pyplot as plt
In [2]:
from util import set_style
set_style()
Out[2]:
In [3]:
# częstotliwość próbkowania
fs = 48000

Specjalne zastosowania filtrów cyfrowych

W poprzednich rozdziałach zostały przedstawione przykłady, w których filtry cyfrowe były wykorzystywane do tłumienia lub przepuszczania wybranych zakresów częstotliwości. W tym rozdziale zostaną omówione bardziej specyficzne zastosowania wybranych typów filtrów cyfrowych.

Wygładzanie sygnałów

Sygnały pozyskiwane z różnego rodzaju czujników (sensorów) są niemal zawsze zniekształcone przez szum. Aby uzyskać sygnał użyteczny, należy pozbyć się szumu. Jest to złożone zagadnienie, szum często ma zwykle niestacjonarny charakter (zmienny w czasie), do jego usuwania wykorzystywane są specjalne algorytmy estymujące szum. Tutaj zostanie pokazane prostsze podejście, polegające na wygładzaniu sygnału (smoothing). Nie usuwa ono w pełni szumu, a jedynie redukuje jego wpływ na sygnał poprzez wygładzenie jego przebiegu.

Poniżej pokazano przykładowe dane uzyskane z rzeczywistego czujnika. Widoczne jest zaszumienie sygnału, chcielibyśmy uzyskać sygnał oczyszczony z szumu, tak aby reprezentował on rzeczywiste zmiany mierzonej wartości.

In [4]:
zaszumiony = np.loadtxt('pliki/noisy_data.txt')
plt.plot(zaszumiony)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Zaszumiony sygnał z czujnika')
plt.show()

Najprostszym filtrem wygładzającym sygnał jest filtr FIR o długości $L$, którego wszystkie współczynniki są równe $1/L$:

$$y(n) = \frac{1}{L} \left(x(n) + x(n-1) + ... + x (n-L+1)\right)$$

Taki filtr nazywa się filtrem średniej ruchomej (MA, moving average filter). Długość filtru określa liczbę próbek sygnału uwzględnianych przy uśrednianiu. Poniżej pokazano charakterystyki częstotliwości filtrów o długości 5, 11, i 31. W kazdym przypadku widoczne są zafalowania charakterystyki. W miarę zwiększania długości filtru, maleje szerokość listka głównego przy zerowej częstotliwości (reprezentującej rzeczywistą średnią sygnału), maleje wzmocnienie dla większych częstotliwości oraz zwiększa się częstotliwość zafalowań.

In [5]:
for L in (5, 11, 31):
    w, h = sig.freqz(np.ones(L) / L, 1)
    plt.plot(w * fs / (2 * np.pi), 20 * np.log10(np.abs(h)), label='L={}'.format(L))
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('wzmocnienie [dB]')
plt.title('Charakterystyki filtrów średniej ruchomej')
plt.ylim(bottom=-50)
plt.legend()
plt.show()

Na kolejnym wykresie pokazano powiększony fragment sygnału i przebiegi sygnału po filtracji. Dla filtru o długości 5 nastąpiło zmniejszenie szumu, ale jest on nadal obecny. Filtr o długości 11 daje dość dobre wygładzenie sygnału. Filtr o długości 31 uśrednia zbyt duży fragment sygnału, wkutek czego wygładzenie jest nadmierne i tracimy szczegóły przebiegu.

In [6]:
plt.plot(zaszumiony, c='#404040', label='oryginalny')
for L in (5, 11, 31):
    b = np.ones(L) / L
    odszumiony = sig.filtfilt(b, 1, zaszumiony)
    plt.plot(odszumiony, label='L={}'.format(L))
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Wygładzanie sygnału z czujnika')
plt.xlim(250, 320)
plt.ylim(28, 50)
plt.legend()
plt.show()

Stosuje się również ważone filtry średniej ruchomej, w których współczynniki nie są jednakowe, można np. nadać młodszym próbkom większą wagę niż starszym.

Innym filtrem stosowanym do wygładzania sygnałów jest filtr Sawickiego-Golaya. Podczas gdy poprzedna metoda jedynie uśrednia odcinki sygnału, filtr S-G dopasowuje wielomian zadanego rzędu do poszczególnych odcinków sygnału, stosując metodę najmniejszych kwadratów. Fitry S-G, w porównaniu z filtrami średniej ruchomej, pozwalają zachować więcej szczegółów sygnału (zmian amplitudy nie będących szumem), wymagają jednak znacząco dłuższych obliczeń.

Dla filtrów S-G musimy podać dwa parametry. Długość okna wyznacza liczbę współczynników i określa liczbę próbek sygnału uwzględnianych przy wygładzaniu, długość okna musi być liczbą nieparzystą. Drugim parametrem jest rząd wielomianu dopasowywanego do punktów, musi on być mniejszy niż długość okna. Wielomiany rzędu trzeciego zazwyczaj dają tu dobre wyniki. Wykres charakterystyki filtru S-G jest dość podobny do filtru średniej ruchomej. Współczynniki filtru zostały obliczone za pomocą funkcji scipy.signal.savgol_coeffs.

In [7]:
for L in (5, 11, 31):
    b = sig.savgol_coeffs(L, 3)
    w, h = sig.freqz(b, 1)
    plt.plot(w * fs / (2 * np.pi), 20 * np.log10(np.abs(h)), label='L={}'.format(L))
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('wzmocnienie [dB]')
plt.title('Charakterystyki filtrów Sawickiego-Golaya')
plt.ylim(bottom=-50)
plt.legend()
plt.show()

Poniżej przedstawiono wynik wygładzania sygnału filtrem S-G o różnej długości. Dla długości równej 5, szum jest usunięty w niewielkim stopniu. Dla długości 11 sygnał zostaje wygładzony, ale nadal podąża za kształtem zaszumionego sygnału. Filtr o długości 31 daje pożądany stopień wygładzenia sygnału.

In [8]:
plt.plot(zaszumiony, c='#404040', label='oryginalny')
for L in (5, 11, 31):
    odszumiony = sig.savgol_filter(zaszumiony, L, 3)
    plt.plot(odszumiony, label='L={}'.format(L))
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Wygładzanie sygnału filtrem Sawickiego-Golaya')
plt.xlim(250, 320)
plt.ylim(28, 50)
plt.legend()
plt.show()

Kolejny wykres pokazuje porówanie działania filtru średniej ruchomej o długości 11 z filtrem S-G o długości 31. Przebiegi są dość podobne, ale filtr S-G stara się zachować szczegóły przebiegu, podczas gdy filtr średniej ruchomej daje bardzo gładki przebieg.

In [55]:
plt.plot(zaszumiony, c='#404040', label='oryginalny')
odszumiony_sr = sig.filtfilt(np.ones(11) / 11, 1, zaszumiony)
plt.plot(odszumiony_sr, label='ŚR L=11')
odszumiony_sg = sig.savgol_filter(zaszumiony, 31, 3)
plt.plot(odszumiony_sg, label='SG L=31')
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Porównanie filtrów wygładzających')
plt.xlim(250, 320)
plt.ylim(28, 50)
plt.legend(loc='upper right')
plt.show()

Ostatni wykres pokazuje cały sygnał w wersji oryginalnej oraz po wygładzeniu za pomocą opisywanych filtrów. W tym przypadku efekty są niemal identyczne, więc prostszy filtr średniej ruchomej sprawdził się dobrze w wygładzaniu sygnału. Gdyby jednak oryginalny sygnał miał większą zmienność (więcej składowych o wysokich częstotliwościach), filtr S-G mógłby dać lepszy wynik.

In [10]:
fig, ax = plt.subplots(3, sharex=True)
ax[0].plot(zaszumiony, label='zaszumiony')
ax[1].plot(odszumiony_sr, label='ŚR L=11')
ax[2].plot(odszumiony_sg, label='SG L=31')
for a in ax:
    a.set_ylabel('amplituda')
    a.legend()
ax[-1].set_xlabel('nr próbki')
fig.suptitle('Sygnał zaszumiony i odszumione')
plt.show()

Filtr różniczkujący

Zadaniem filtru różniczkowego (differentiator) jst obliczanie pochodnej sygnału. W obwodach analogowych, rolę układów różniczkujących pełnią układy zbudowane z rezystorów i cewek (RL). Cyfrowe filtry różniczkujące mogą mieć zastosowanie w przypadkach, gdw których pochodna sygnału jest potrzebna do uzyskania żądanych informacji. Przykładowo, mając sygnał opisujący zmiany położenia obiektu, możemy uzyskać informację o prędkości poruszania się tego obiektu. W analizie sygnałów, obliczenie pochodnej sygnału może być przydatne przy obliczaniu ekstremów sygnału.

Najprostszym cyfrowym filtrem różniczkującym jest filtr FIR pierwszego rzędu, który oblicza różnicę między bieżącą a poprzednią próbką:

$$y(n) = x(n) - x(n-1)$$

Taki efekt możemy uzyskać za pomocą filtru, albo za pomocą funkcji numpy.diff, która liczy różnicę "w przód".

Innym przykładem filtru różniczkującego jest filtr drugiego rzędu, dany równaniem różnicowym:

$$y(n) = \frac{1}{2}\left(x(n) - x(n-2)\right)$$

Oba filtry , nazywane filtrami różnicowymi, nie są jednak cyfrowym odpowiednikiem analogowych układów różniczkujących. Jeżeli filtr ma obliczać pochodną sygnału, to po podaniu na jego wejście sygnału $\sin(x)$ powinniśmy uzyskać na wyjściu sygnał $x \cos(x)$, czyli jego charakterystyka częstotliwościowa powinna narastać liniowo z częstotliwością, a charakterystyka fazowa powinna być stała i równa $\pi/2$. Opisane wyżej filtry mają charakterystyki znacznie odbiegające od oryginalnej, co pokazano na poniższym wykresie. Filtr działa zgodnie z założeniami tylko w paśmie do ok. 5 kHz. Będzie więc różniczkował składowe sygnału o niskich częstotliwościach, natomiast wyższe częstotliwości zostaną stłumione. Może to być pożądanym efektem, jeżeli nie chcemy uwzględniać szumu przy obliczaniu pochodnej.

In [11]:
w_d1, h_d1 = sig.freqz([1, -1], 1)
plt.plot(w_d1 * fs / (2 * np.pi), np.linspace(0, np.pi, len(w_d1)), '--', c='#808080', label='idealna char.')
plt.plot(w_d1 * fs / (2 * np.pi), np.abs(h_d1), label='filtr 1. rzędu')
w_d2, h_d2 = sig.freqz([0.5, 0, -0.5], 1)
plt.plot(w_d2 * fs / (2 * np.pi), np.abs(h_d2), label='filtr 2. rzędu')
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('wzmocnienie')
plt.title('Charakterystyki prostych filtrów różniczkujących')
plt.legend()
plt.show()

Możliwe jest także stosowanie filtrów różnicowych wyższego rzędu.

Nie da się zrealizować cyfrowego filtru mającego idealną charakterystykę układu różniczkującego, można ją jednak przybliżyć stosując metody optymalizacyjne w projektowaniu filtru, np. metodę Parksa-McClellana. Do realizacji takiego filtru musimy wykorzystać filtr FIR typu III lub IV, czyli o asymetrycznej odpowiedzi impulsowej. Do zaprojektowania filtru możemy użyć funkcji scipy.signal.remez, podając parametr type='differentiator'. Musimy założyć zakres częstotliwości, w którym charakterystyka częstotliwościowa ma być liniowa. Załóżmy zakres przepustowy od 0 do 23 kHz i pasmo przejściowe prawie do częstotliwości Nyquista. Wzmocnienie w paśmie przepustowym ustawiamy na równe wzmocnieniu dla częstotliwości próbkowania, czyli $2\pi$, algorytm sam przeskaluje wzmocnienia dla poszczególnych częstotliwości.

In [12]:
plt.plot(w_d1 * fs / (2 * np.pi), np.linspace(0, np.pi, len(w_d1)), '--', c='#808080', label='idealna char.')
for L in (31, 41, 51, 61): 
    b_d3 = sig.remez(L, (0, 23, 23.99, 24), (2 * np.pi, 0), type='differentiator', fs=48)
    w_d3, h_d3 = sig.freqz(b_d3, 1)
    plt.plot(w_d3 * fs / (2 * np.pi), np.abs(h_d3), label='L={}'.format(L))
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('wzmocnienie')
plt.title('Charakterystyki filtrów różniczkujących (remez)')
plt.legend()
plt.show()

Jak widać, zwiększanie długości filtru powoduje zbliżenie charakterystyki do idealnej. Dla długości równej 61 mamy już w miarę równą charakterystykę.

Poniżej pokazano wynik różniczkowania sygnału, który wcześniej został poddany wygładzaniu. Sygnały zostały wyrównane i przeskalowane tak, aby można było je porównać. Widoczne jest, że miejsca zerowe zróżniczkowanego przebiegu pokrywają się z minimami i maksimami oryginalnego sygnału. Analizowany sygnał jest wolnozmienny, a w paśmie niskich częstotliwości, prosty filtr pierwszego rzędu ma w przybliżeniu liniową charakterystykę, dlatego oba filtry dają niemal identyczny wynik. Dla sygnałów zawierających składowe o wysokich częstotliwościach, filtr pierwszego błędu mógłby wprowadzać znaczne blędy do wyniku.

In [13]:
rozniczka_rzad1 = sig.lfilter([1, -1], 1, odszumiony_sg)
rozniczka_remez = sig.lfilter(b_d3, 1, odszumiony_sg)
plt.plot(odszumiony_sg[1:] / 10, '#c0c0c0', label='oryginalny')
plt.plot(rozniczka_rzad1, label='filtr 1. rzędu')
plt.plot(rozniczka_remez[30:], label='remez L=61')
plt.axhline(0, lw=1, c='#404040')
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Wynik różniczkowania sygnału')
plt.legend()
plt.ylim(top=5)
plt.show()

Filtr całkujący

Filtr całkujący (integrator) oblicza całkę sygnału, jego działanie jest zatem odwrotne do filtru różniczkującego. W obwodach analogowych, układy różniczkujące składają się z rezystorów i kondensatorów (RC). W dziedzinie cyfrowej, całkowanie sygnału odpowiada sumowaniu dyskretnych wartości. Filtr całkujący musi mieć realizację IIR. Filtr całkujący pierwszego rzędu ma równanie różnicowe:

$$y(n) = y(n-1) + x(n)$$

Jest to tzw. filtr całkujący wstecz (backward rectangular integrator). Inna forma filtru całkującego:

$$y(n) = y(n-1) + x(n-1)$$

to tzw. filtr całkujący w przód (forward rectangular integrator). Trzecia forma filtru pierwszego rzędu to:

$$y(n) = y(n-1) + \frac{1}{2} \left(x(n) + x(n-1) \right)$$

Jest to tzw. filtr trapezoidalny (trapezoidal integrator).

Charakterystyki tych filtrów pokazano na poniższym wykresie. Filtr "w przód" i "wstecz" mają identyczne charakterystyki, leżące nieco powyżej idealnej charakterystyki $1/j\omega$. Charakterystyka filtru trapezoidalnego leży nieco poniżej idealnej. W celu obliczenia charakterystyk należy podać wprost wartości pulsacji (zmienna w), domyślnie funkcja freqz bierze punkty w zakresie od 0 do $\pi$, a charakterystyka filtru różniczkującego jest nieokreślona dla zerowej częstotliwości.

In [14]:
w = np.linspace(0.1 * np.pi, np.pi, 512)
w_r1, h_r1 = sig.freqz([1], [1, -1], w)
plt.plot(w_r1 * fs / (2 * np.pi), 1 / w, '--', c='#808080', label='idealna char.')
plt.plot(w_r1 * fs / (2 * np.pi), np.abs(h_r1), label='filtr r. wstecz')
w_r2, h_r2 = sig.freqz([0, 1], [1, -1], w)
plt.plot(w_r2 * fs / (2 * np.pi), np.abs(h_r2), label='filtr r. w przód')
w_r3, h_r3 = sig.freqz([0.5, 0.5], [1, -1], w)
plt.plot(w_r3 * fs / (2 * np.pi), np.abs(h_r3), label='filtr r. trapezoidalny')
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('wzmocnienie')
plt.title('Charakterystyki filtrów całkujących 1. rzędu')
plt.legend()
plt.show()

Wszystkie omówione powyżej filtry całkujące mają jedną wspólną cechę: są to filtry niestabilne. Łatwo zauważyć, że jeżeli na wejście filtru podamy sygnał o stałych wartościach, odpowiedź filtru będzie dążyła do nieskończoności. Jeżeli podamy sygnał przyjmujący tylko wartości dodatnie, obliczana suma będzie coraz większa.

Rozwiązaniem problemu jest wprowadzenie do równania różnicowego współczynnika $\alpha < 1$. Dla filtru całkującego wstecz, równanie różnicowe przyjmuje postać:

$$y(n) = x(n) + \alpha y(n-1)$$ $$0 < \alpha < 1$$

Układ o takim równaniu różnicowym nazywa się stratnym układem całkującym. Anglojęzyczna nazwa: leaky integrator, czyli "układ całkujący z przeciekiem", obrazowo opisuje efekt jego działania: część obliczonej wcześniej sumy "wycieka" z układu. Modyfikacja ta zapewnia, że moduł biegunów jest mniejszy od jedności i układ jest stabilny. Transmitancja takiego filtru ma postać:

$$H(z) = \frac{1}{1 - \alpha z^{-1}}$$

Wartości $\alpha$ stosowane w praktyce są zwykle bardzo bliskie jedności, najczęściej w zakresie 0,99-0,998. Poniżej pokazano charakterystyki częstototliwościowe stratnego układu całkującego.

In [15]:
w_r4, h_r4 = sig.freqz([1], [1, -0.998])
plt.plot(w_r4 * fs / (2 * np.pi), np.abs(h_r4), label='alpha=0.998')
w_r5, h_r5 = sig.freqz([1], [1, -0.9])
plt.plot(w_r5 * fs / (2 * np.pi), np.abs(h_r5), label='alpha=0.990')
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('wzmocnienie')
plt.title('Charakterystyka stratnego filtru całkującego "wstecz"')
plt.ylim(0, 10)
plt.legend()
plt.show()

Zmniejszanie parametru $\alpha$ powoduje większe tłumienie charakterystyki dla niskich częstotliwości. W podobny sposób można zmodyfikować charakterystykę filtru całkującego "w przód" oraz trapezoidalnego.

Wynik całkowania sygnału, który w poprzednim podrozdziale uzyskano przez różniczkowanie wygładzonych danych z czujnika, wygląda następująco.

In [16]:
calka = sig.lfilter([1], [1, -0.998], rozniczka_remez[30:])
plt.plot(odszumiony_sg, '#c0c0c0', label='oryginalny')
plt.plot(calka, label='po całkowaniu')
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Wynik całkowania sygnału')
plt.legend()
plt.show()

Widać, że po różniczkowaniu, a następnie całkowaniu, otrzymujemy sygnał odpowiadający oryginalnemu sygnałowi, aczkolwiek sygnały nie są identyczne na skutek utraty części informacji w trakcie przetwarzania. Przesunięcie sygnałów na osi pionowej wynika z tego, że składowa stała występująca w oryginalnym sygnale została wytłumiona przez filtr różniczkujący, co zostało opisane w kolejnym podrozdziale.

Filtr składowej stałej

Wiele sygnałów pozyskiwanych z czujników zawiera składową stałą, czyli są one przesunięte o pewną wartość na skali amplitudy. Składowa stała jest niepożądana i powinna być odfiltrowana z sygnału. Jeżeli mamy do dyspozycji kompletny sygnał, wystarczy obliczyć składową stałą i odjąć ją od sygnału. W praktycznych zastosowaniach przetwarzamy jednak sygnał ciągły. W widmie sygnału, składowa stała jest reprezentowana przez składową na zerowej częstotliwości. Filtr składowej stałej (DC blocker, DC oznacza składową stałą, direct component) jest filtrem górnoprzepustowym, który usuwa składową stałą z sygnału. Filtr nie powinien jednak tłumić składowych sygnału o niskich częstotliwościach.

Typowy filtr cyfrowy usuwający składową stałą powstaje poprzez połączenie filtru różniczkującego z filtrem całkującym. Filtr różniczkujący daje zerowe wzmocnienie składowej stałej, ale silnie tłumi niższe częstotliwości. Stratny filtr całkujący z kolei wzmacnia niskie częstotliwości. Wzmocnienie filtru różniczkującego zależy wprost od częstotliwości, natomiast filtru całkującego – od odwrotności częstotliwości, zatem te dwa filtry kompensują nawzajem swoje wzmocnienie, za wyjątkiem częstotliwości bliskich zeru, gdzie wzmocnienie wynikowego filtru zostaje zredukowane.

Transmitancja filtru składowej stałej pierwszego rzędu ma postać:

$$H(z) = \frac{1 - z^{-1}} {1 - \alpha z^{-1}}$$

Znaczenie parametru $\alpha$ jest takie samo jak dla filtru całkującego. Charakterystykę filtru dla $\alpha$=0,998 pokazano poniżej.

In [17]:
w_dc, h_dc = sig.freqz([1, -1], [1, -0.998], np.linspace(0, 2 * np.pi * 1000 / fs, 512))
plt.plot(w_dc * fs / (2 * np.pi), np.abs(h_dc))
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('wzmocnienie')
plt.title('Charakterystyka filtru składowej stałej')

plt.figure()
w_dc2, h_dc2 = sig.freqz([1, -1], [1, -0.998], np.linspace(0, 2 * np.pi * 200 / fs, 512))
plt.axhline(1, ls='--', lw=1, c='#808080')
plt.plot(w_dc2 * fs / (2 * np.pi), np.abs(h_dc2))
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('wzmocnienie')
plt.title('Charakterystyka filtru składowej stałej (niskie cz.)')

plt.show()

Widać, że dla składowych poniżej 150 Hz wprowadzane jest tłumienie.

Wynik usuwania składowej stałej z sygnału pochodzącego z czujnika pokazano na poniższym wykresie. Obliczenie średniej z próbek sygnału po filtracji pokazuje, że nie jest ona zerowa. Pamiętajmy, że stratny układ całkujący ma charakterystykę odbiegającej od idealnej. W celu dokładniejszego usunięcia składowej stałej należy stosować filtry wyższego rzędu.

In [18]:
print('Składowa stała przed filtracją = {}'.format(np.mean(zaszumiony)))
bez_dc = sig.lfilter([1, -1], [1, -0.998], zaszumiony)
print('Składowa stała po filtracji    = {}'.format(np.mean(bez_dc)))
plt.plot(zaszumiony, '#c0c0c0', label='oryginalny')
plt.plot(bez_dc, label='po filtracji')
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Wynik usuwania składowej stałej')
plt.legend()
plt.show()

Składowa stała przed filtracją = 1.824
Składowa stała po filtracji    = -0.43858145498372925

Filtr grzebieniowy

Filtr grzebieniowy (comb filter) jest to filtr o cyklicznie powtarzającej się charakterystyce amplitudowej, przypominającej kształtem grzebień. Filtr grzebieniowy typu FIR można uzyskać z filtru FIR pierwszego rzędu poprzez zwiększenie opóźnienia do wartości $L$, tak że równanie różnicowe przyjmuje postać:

$$y(n) = x(n) + \alpha x(n-L)$$

gdzie $|\alpha| \le 1$. Działanie filtru polega zatem na dodaniu do sygnału jego "echa", czyli opóźnionej i przeskalowanej kopii. Charakterystyki amplitudowe dla różnych wartości $L$ oraz $\alpha=1$ pokazano poniżej (współczynniki filtru są skalowane aby zachować amplitudę sygnału w zakresie [-1, 1]).

In [19]:
fig, ax = plt.subplots(3, sharex=True, figsize=(8, 8))
w, h1 = sig.freqz([0.5, 0.5], 1)
f = w * fs / (2 * np.pi)
ax[0].plot(f, np.abs(h1), label='L=1')
_, h3 = sig.freqz([0.5, 0, 0, 0.5], 1)
ax[1].plot(f, np.abs(h3), label='L=3')
_, h10 = sig.freqz([0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5], 1)
ax[2].plot(f, np.abs(h10), label='L=10')
ax[-1].set_xlabel('częstotliwość [Hz]')
for a in ax:
    a.set_ylabel('wzmocnienie')
    a.legend(loc='upper right')
plt.suptitle('Charakterystyki amplitudowe filtrów grzebieniowych')
plt.show()

Filtr grzebieniowy o opóźnieniu $L$ ma $L$ minimów w zakresie $\left[0, f_s \right)$. Częstotliwości minimów przy $\alpha>0$ są równe:

$$f_k = \frac{2k-1}{2L} f_s$$

natomiast dla $\alpha<0$:

$$f_k = \frac{k-1}{L} f_s$$

$$k = 1, 2, ..., L$$

Minima charakterystyki odpowiadają położeniu zer transmitancji filtru: $L$ zer jest rozłożonych w równych odległościach na okręgu o promieniu $|\alpha|$. Dla filtru o $L=10$ i $\alpha=1$, częstotliwości minimów są następujące:

In [20]:
z, p, k = sig.tf2zpk([0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5], 1)
print('alpha = 1')
print('Kąty fazowe:', np.sort(np.degrees(np.angle(z))))
print('Częstotliwości:', np.sort(fs * np.angle(z) / (2 * np.pi)))

alpha = 1
Kąty fazowe: [-162. -126.  -90.  -54.  -18.   18.   54.   90.  126.  162.]
Częstotliwości: [-21600. -16800. -12000.  -7200.  -2400.   2400.   7200.  12000.  16800.
  21600.]

Dla ujemnego $\alpha$, minima przesuwają się na osi częstotliwości:

In [21]:
z, p, k = sig.tf2zpk([0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.5], 1)
print('alpha = -1')
print('Kąty fazowe:', np.sort(np.degrees(np.angle(z))))
print('Częstotliwości:', np.sort(fs * np.angle(z) / (2 * np.pi)))

alpha = -1
Kąty fazowe: [-144. -108.  -72.  -36.    0.   36.   72.  108.  144.  180.]
Częstotliwości: [-19200. -14400.  -9600.  -4800.      0.   4800.   9600.  14400.  19200.
  24000.]
In [22]:
w, h10_plus = sig.freqz([0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5], 1)
f = w * fs / (2 * np.pi)
plt.plot(f, np.abs(h10_plus), label=r'$\alpha=1$')
_, h10_minus = sig.freqz([0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.5], 1)
plt.plot(f, np.abs(h10_minus), label=r'$\alpha=-1$')
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('wzmocnienie')
plt.title(r'Charakt. filtru grzebieniowego dla dodatniego i ujemnego $\alpha$')
plt.legend()
plt.show()

Filtr grzebieniowy może być wykorzystany do odfiltrowania składowych widmowych na częstotliwościach harmonicznych (będących wielokrotnościami częstotliwości podstawowej). Dla przykładu, wygenerujmy sygnał harmoniczny o częstotliwości podstawowej 1 kHz (wykorzystamy do tego celu sygnału piłokształtny i zignorujemy aliasing), po czym użyjmy filtru grzebieniowego do usunięcia składowych o częstotliwościach 2 kHz i ich wielokrotnościach. Przy częstotliwości próbkowania równej 48 kHz, $48/2=24$, więc potrzebny jest filtr o opóźnieniu $L=24$ i ujemnym współczynniku $\alpha$. Gdyby stosunek ten nie był liczbą całkowitą, konieczne byłoby użycie filtru o ułamkowym opóźnieniu (z interpolacją). Współczynnik $b_{24}$ został zmniejszony do 0,499 aby uniknąć wartości zero przy obliczaniu logarytmu.

In [23]:
pila = sig.sawtooth(2 * np.pi * np.arange(2048) * 1000 / fs)
w_pila = 20 * np.log10(np.abs(np.fft.rfft(pila * np.hamming(2048))) / 1024)
plt.plot(np.fft.rfftfreq(2048, 1 / fs), w_pila)
b_24 = np.zeros(25)
b_24[0] = 0.5
b_24[-1] = -0.499
w, h_24 = sig.freqz(b_24, 1)
plt.plot(w * fs / (2 * np.pi), 20 * np.log10(np.abs(h_24)))
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('wzmocnienie')
plt.xlim(0, 11000)
plt.ylim(bottom=-40)
plt.title('Sygnał harmoniczny i filtr grzebieniowy L=24')
plt.show()

Poniżej pokazano wynik filtracji sygnału.

In [24]:
pila_filtr = sig.lfilter(b_24, 1, pila)
w_pila_filtr = 20 * np.log10(np.abs(np.fft.rfft(pila_filtr * np.hamming(2048))) / 1024)
plt.plot(np.fft.rfftfreq(2048, 1 / fs), w_pila_filtr)
plt.xlim(0, 11000)
plt.ylim(bottom=-40)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('wzmocnienie')
plt.title('Widmo sygnału harmonicznego po filtracji grzebieniowej')
plt.figure()
plt.plot(pila_filtr[:500])
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał harmoniczny po filtracji grzebieniowej')
plt.show()

Filtr grzebieniowy usunął z sygnału parzyste wielokrotności częstotliwości podstawowej i pozostawił tylko nieparzyste wielokrotności. Zgodnie z teorią, tak uzyskane widmo reprezentuje sygnał prostokątny, co potwierdza wykres czasowy przefiltrowanego sygnału.

Filtry grzebieniowe stosuje się również w przetwarzaniu dźwięku, do uzyskania efektu brzmieniowego nazywanego opóźnieniem (delay) dla małych wartości opóźnień, gdy kopia "zlewa się" z oryginalnym dźwiękiem, oraz jednokrotnym echem dla większych czasów opóźnień (>50 ms), gdy można wyraźnie usłyszeć dwie osobne kopie dźwięku.

Rekursywny filtr grzebieniowy ma formę filtru IIR, można go uzyskać ze "zwykłego" filtru IIR poprzez zmianę wartości opóźnienia. Równanie różnicowe ma postać:

$$y(n) = x(n) + \alpha y(n-L)$$

gdzie $|\alpha|\le 1$, przy czym dla $\alpha=1$ filtr jest na granicy stabilności. Charakterystyka amplitudowa filtru dla $L=10$ i $\alpha=\pm 1$ wygląda następująco.

In [25]:
w, h10_plus_rek = sig.freqz([0.5], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5])
f = w * fs / (2 * np.pi)
plt.plot(f, np.abs(h10_plus_rek), label=r'$\alpha=1$')
_, h10_minus_rek = sig.freqz([0.5], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.5])
plt.plot(f, np.abs(h10_minus_rek), label=r'$\alpha=-1$')
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('wzmocnienie')
plt.title(r'Charakt. rekursywnego filtru grzebieniowego dla dodatniego i ujemnego $\alpha$')
plt.legend()
plt.show()

O ile filtr grzebieniowy FIR miał charakterystykę przypominającą złożenie wielu filtrów środkowozaporowych o wąskim paśmie zaporowym, tak tutaj mamy efekt złożenia wielu filtrów środkowoprzepustowych o wąskim paśmie przepustowym. "Grzebień" ma w tym przypadku ostrza skierowane w górę. W przypadku filtru rekursywnego, częstotliwości maksimów widmowych są wyznaczone przez położenie biegunów. Filtr o opóźnieniu $L$ ma $L$ biegunów rozłożonych równomiernie na okręgu o promieniu $|\alpha|$ (zależności są takie same, jak dla zer filtru grzebieniowego FIR).

In [26]:
z, p, k = sig.tf2zpk([0.5], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5])
print('alpha = 1')
print('Kąty fazowe:', np.sort(np.degrees(np.angle(p))))
print('Częstotliwości:', np.sort(fs * np.angle(p) / (2 * np.pi)))
z, p, k = sig.tf2zpk([0.5], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.5])
print('alpha = -1')
print('Kąty fazowe:', np.sort(np.degrees(np.angle(p))))
print('Częstotliwości:', np.sort(fs * np.angle(p) / (2 * np.pi)))

alpha = 1
Kąty fazowe: [-162. -126.  -90.  -54.  -18.   18.   54.   90.  126.  162.]
Częstotliwości: [-21600. -16800. -12000.  -7200.  -2400.   2400.   7200.  12000.  16800.
  21600.]
alpha = -1
Kąty fazowe: [-144. -108.  -72.  -36.    0.   36.   72.  108.  144.  180.]
Częstotliwości: [-19200. -14400.  -9600.  -4800.      0.   4800.   9600.  14400.  19200.
  24000.]

Interesująco wygląda odpowiedź impulsowa filtru grzebieniowego IIR. Poniżej przykład dla $L=10$ i $\alpha=0.8$.

In [27]:
odp_imp = sig.lfilter([1], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.8], sig.unit_impulse(150))
plt.stem(odp_imp)
plt.show()

Widać na wykresie, że kolejne kopie odpowiedzi impulsowej pojawiają się co $L$ próbek, każda kolejna kopia ma $\alpha$ razy mniejszą amplitudę niż poprzednia, a ponadto faza co drugiej kopii jest odwrócona. Dla $\alpha=1$ mamy filtr na granicy stabilności, kopie impulsu o jednakowej amplitudzie będą generowane w nieskończoność.

Rekursywny filtr grzebieniowy jest także stosowany w przetwarzaniu dźwięku, umożliwia on uzyskanie efektu wielokrotnego echa (multiple echo). W niektórych implementacjach, efekt ten jest określany jako pogłos (reverberation, reverb), ale nie jest to poprawna nazwa, gdyż właściwy pogłos, symulujący odbicia dźwięku w pomieszczeniu, nie ma charakteru kopii o stałych przesunięciach czasowych, liczba kopii wzrasta z czasem, a odstępy między nimi zmniejszają się.

Filtr wszechprzepustowy

Filtr wszechprzepustowy (all-pass filter) jest dość nietypowym filtrem, ponieważ jego charakterystyka amplitudowa jest równa jedności w całym zakresie częstotliwości. Zadaniem tego filtru jest natomiast modyfikacja fazy sygnału, bez wpływu na jego amplitudę. Jeżeli popatrzymy na zamieszczone wcześniej wykresy charakterystyk amplitudowych filtrów grzebieniowych – zwykłego i rekursywnego, zobaczymy że jedna jest odwrotnością drugiej, ich iloczyn jest równy jedności. W istocie, połączenie kaskadowe nierekursywnego i rekursywnego filtru grzebieniowego o tym samym opóźnieniu $L$ i o przeciwnych współczynnikach skalujących $\alpha$, daje filtr wszechprzepustowy. Rozpatrzmy filtr pierwszego rzędu ($L=1$). Połączenie obu filtrów daje filtr o transmitancji:

$$H(z) = \frac{-z_0 + z^{-1}} {1 - z_0 z^{-1}}$$

Wartość $z_0$ wyznacza położenie bieguna transmitancji, jest ona rzeczywista, natomiast położenie zera transmitancji określa wartość $1 / z_0$. Tak jest zawsze dla filtrów wszechprzepustowych: zera i bieguny tworzą pary, jedno z nich leży wewnątrz okręgu jednostkowego, drugie na zewnątrz, a moduł jednego jest odwrotnością drugiego. Poniżej pokazano przykład dla $z_0=0{,}5$, wykreślono charakterystykę amplitudową oraz fazową.

In [28]:
z0 = 0.5
b = [-z0, 1]
a = [1, -z0]
z, p, k = sig.tf2zpk(b, a)
print('Zero:', z)
print('Biegun:', p)
w, h = sig.freqz(b, a)
f = w * fs / (2 * np.pi)
plt.plot(f, np.abs(h))
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('wzmocnienie')
plt.title('Char. amplitudowa filtru wszechprzepustowego 1. rzędu')
plt.figure()
plt.plot(f, np.angle(h))
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('faza [rad]')
plt.title('Char. fazowa filtru wszechprzepustowego 1. rzędu')
plt.show()

Zero: [2.]
Biegun: [0.5]

Charakterystyka amplitudowa jest rzeczywiście równa 1 w całym paśmie, charakterystyka fazowa zależy od częstotliwości.

W praktycznych zastosowaniach, chcemy uzyskać filtr wszechprzepustowy, którego zera i bieguny są położone w ściśle określonych miejscach. Jeżeli $z_0$ jest wartością zespoloną, transmitancja filtru wszechprzepustowego ulega przekształceniu do następującej postaci:

$$H(z) = \frac{\left|z_0\right|^2 - 2\Re(z_0)z^{-1} + z^{-2}} {1 - 2\Re(z_0)z^{-1} + \left|z_0\right|^2 z^{-2}}$$

Jest to filtr IIR drugiego rzędu. Ma on parę zespolonych, sprzężonych biegunów $z_0$ i $z_0^*$ oraz parę zespolonych, sprzężonych zer $1/z_0^*$ i $1/z_0$. Zauważmy, że w liczniku i w mianowniku transmitancji występują dokładnie te same współczynniki, tylko w odwrotnej kolejności. Jest to regułą dla filtrów wszechprzepustowych.

Poniżej pokazano przykład obliczeń filtru wszechprzepustowego drugiego rzędu, którego bieguny i zera mają kąt fazowy ±45 stopni (częstotliwości $\pm f_s/4$).

In [29]:
z0 = 0.5 * np.exp(1j * np.pi / 4)
b = [np.abs(z0)**2, -2 * np.real(z0), 1]
a = b[::-1]
z, p, k = sig.tf2zpk(b, a)
print('Zera:', z, 'moduł:', np.abs(z[0]))
print('Bieguny:', p, 'moduł:', np.abs(p[0]))
w, h = sig.freqz(b, a)
f = w * fs / (2 * np.pi)
plt.plot(f, np.angle(h))
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('faza [rad]')
plt.title('Char. fazowa filtru wszechprzepustowego 2. rzędu')
plt.show()

Zera: [1.41421356+1.41421356j 1.41421356-1.41421356j] moduł: 2.0
Bieguny: [0.35355339+0.35355339j 0.35355339-0.35355339j] moduł: 0.5

Charakterystyka fazowa jest podobna do uzyskanej dla poprzedniego filtru wszechprzepustowego, jest jednak przesunięta. Dzięki temu można wprowadzić do sygnału pożądaną korekcję fazy.

Filtry wszechprzepustowe są często łączone kaskadowo z filtrami IIR. Wypadkowa charakterystyka fazowa układu filtrów jest sumą charakterystyk fazowych obu filtrów. Kształtując w odpowiedni sposób charakterystykę fazową filtru wszechprzepustowego, można zniwelować nieliniowość charakterystyki fazowej filtru IIR i w ten sposób uzyskać w przybliżeniu stałe opóźnienie grupowe układu filtrów w paśmie przepustowym. Charakterystyka amplitudowa filtru IIR nie ulega oczywiście zmianie. Filtry wszechprzepustowe mogą również pomóc w uzyskaniu stabilnego filtru, jeżeli położenie zer i biegunów zostanie tak dobrane, że koryguje ono zera i bieguny filtru IIR.

Filtry wszechprzepustowe są także stosowane do tworzenia efektów dźwiękowych. Jak wspomniano powyżej, filtr grzebieniowy może być użyty do uzyskania efektu opóźnienia dźwięku, jednak ze względu na "grzebieniową" charakterystykę amplitudową, widmo amplitudowe dźwięku jest zniekształcane. Filtry wszechprzepustowe nie mają tej wady, a również są w stanie uzyskać efekt odbicia dźwięku. Filtry wszechprzepustowe mogą również posłużyć do budowy układu realizującego rzeczywisty efekt pogłosu (np. układ Schroedera złożony z kilku filtrów grzebieniowych i wszechprzepustowych).

Filtr Hilberta i sygnał analityczny

Filtrem Hilberta nazywa się filtr, który dokonuje przesunięcia fazowego sygnału rzeczywistego o wartość $\pi$, czyli o 90 stopni. Sygnał analityczny (analytic signal) jest sygnałem zespolonym, którego część urojona jest równa części rzeczywistej przesuniętej o kąt $\pi$. Układ przekształcający sygnał rzeczywisty w sygnał analityczny nazywa się transormatorem Hilberta (Hilbert transformer).

Sygnał analityczny ma jedną ważną cechę. Pamiętamy, że widmo sygnału rzeczywistego jest symetryczne, ma ono dwie kopie na "dodatnich i ujemnych" częstotliwościach. Przekształcenie sygnału rzeczywistego w analityczny powoduje, że pozostaje tylko jedna część widma, druga kopia ("nadmiarowa") dla ujemnych częstotliwości (lub, patrząc inaczej, od częstotliwości Nyquista do częstotliwości próbkowania) zostaje wyzerowana. Funkcja scipy.signal.hilbert tworzy sygnał zespolony operując bezpośrednio na widmie sygnału. Oblicza ona FFT sygnału, zeruje wartości widma dla ujemnych częstotliwości, podwaja amplitudy pozostałych składowych (aby zachować energię sygnału) oraz wykonuje IFFT.

Jako przykład zastosowania sygnału analitycznego, pokażemy przypadek, w którym otrzymujmy sygnał będący sumą dwóch sinusoid, z których jedna ma częstotliwość 1000 Hz, a druga jest mniejsza o pewną wartość. Naszym zadaniem jest znalezienie różnicy częstotliwości obu składowych. Postać czasowa analizowanego sygnału rzeczywistego wygląda następująco.

In [31]:
syg1 = np.sin(2 * np.pi * np.arange(2048) * 1000 / fs)
syg2 = np.sin(2 * np.pi * np.arange(2048) * 800 / fs)
sygnal = syg1 - syg2
plt.plot(sygnal[:1000])
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Różnica dwóch sygnałów sin')
plt.show()

Jest to sygnał zmodulowany amplitudowo. Widmo sygnału zawiera dwa prążki na częstotliwościach: 800 Hz i 1200 Hz. Z takiego widma trudno jest uzyskać żądane informacje. Dlatego przekształcamy ten sygnał w formę analityczną. Widmo amplitudowe sygnału analitycznego wygląda następująco (uwaga: jest to sygnał zespolony, dlatego używamy funkcji fft, a nie rfft jak zazwyczaj).

In [32]:
analityczny = sig.hilbert(sygnal)
w_analityczny = 20 * np.log10(np.abs(np.fft.fft(analityczny * np.hamming(2048))) / 2048)
plt.plot(np.fft.fftfreq(2048, 1 / fs), w_analityczny)
plt.xlabel('częstotliwość [dB]')
plt.ylabel('amplituda [dB]')
plt.title('Widmo amplitudowe sygnału analitycznego')
plt.ylim(bottom=-70)
plt.show()

Widmo nadal zawiera dwa prążki, jak w sygnale rzeczywistym (nie widać ich dokładnie, ponieważ są położone blisko siebie), natomiast zwróćmy uwagę na lewą część widma, reprezentującą ujemne częstotliwości. W widmie sygnału rzeczywistego mielibyśmy również prążki na częstotliwościach -800 Hz i -1200 Hz. Tutaj "ujemna" część widma jest pusta.

Nadal nie widać zysku z zastosowania transformatora Hilberta. Ale widmo modułu sygnału analitycznego pokazuje ciekawe informacje. Poniższy wykres przedstawia moduł widma modułu sygnału analitycznego.

In [33]:
analityczny_m = np.abs(analityczny)
w_analityczny_m = 20 * np.log10(np.abs(np.fft.rfft(analityczny_m * np.hamming(2048))) / 1024)
plt.plot(np.fft.rfftfreq(2048, 1 / fs), w_analityczny_m)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Widmo amplitudowe modułu sygnału analitycznego')
plt.xlim(0, 5000)
plt.ylim(bottom=-70)
plt.show()

Widmo modułu sygnału analitycznego wygląda zupełnie inaczej niż widmo samego sygnału. Jeżeli zastosujemy dowolną metodę pomiaru częstotliwości prążków widma (np. autokorelacyjną), otrzymamy częstotliwość pierwszego prążka równą 200 Hz. Tyle dokładnie wynosiła wartość bezwzględna różnicy częstotliwości obu składowych sinusoidalnych oryginalnego sygnału. Nie da się jednak określić czy czy częstotliwość jednej składowej była mniejsza czy większa o 200 Hz od częstotliwości drugiej składowej.

Jako drugi przykład wykorzystania sygnału analitycznego rozpatrzymy sytuację, w której otrzymujemy dwa sygnały sinusoidalne (tym razem osobne) o jednakowej częstotliwości 200 Hz, których faza różni się o plus lub minus 90 stopni. Naszym zadaniem jest stwierdzenie, czy różnica fazy jest dodatnia czy ujemna. Problem ten można rozwiązać wieloma metodami, np. za pomocą FFT lub funkcji korelacji wzajemnej, tutaj zostanie wykorzystany transformator Hilberta. Najpierw przekształcamy oba sygnały w formę analityczną. Dalej są możliwe dwie metody. W pierwszej z nich oblicza się kąt fazowy stosunku obu sygnałów analitycznych.

In [34]:
syg1 = np.sin(2 * np.pi * np.arange(2048) * 200 / fs)
syg2 = np.sin(2 * np.pi * np.arange(2048) * 200 / fs + np.pi / 2)

an_syg1 = sig.hilbert(syg1)
an_syg2 = sig.hilbert(syg2)
r_fazy = np.degrees(np.angle(an_syg2 / an_syg1))
print('Mediana różnicy fazy:', np.median(r_fazy))
plt.plot(r_fazy)
plt.xlabel('nr próbki')
plt.ylabel('różnica fazy')
plt.title('Różnica fazy sygnałów')
plt.show()

Mediana różnicy fazy: 90.00712115202145

Charakterystyka jest dość nierówna (im dłuższy sygnał, tym gładsza charakterystyka), ale wartość mediany, jak i wartości w środku charakterystyki, wskazują na to, że drugi sygnał jest przesunięty o +90 stopni względem drugiego.

Sprawdźmy teraz przypadek, gdy różnica fazy wynosi -90 stopni. Użyjemy teraz drugiego podejścia: obliczamy rozwiniętą fazę każdego z sygnałów analitycznych i odejmujemy ją od siebie.

In [35]:
syg1 = np.sin(2 * np.pi * np.arange(2048) * 200 / fs)
syg2 = np.sin(2 * np.pi * np.arange(2048) * 200 / fs - np.pi / 2)

an_syg1 = sig.hilbert(syg1)
an_syg2 = sig.hilbert(syg2)
r_fazy = np.degrees(np.unwrap(np.angle(an_syg2)) - np.unwrap(np.angle(an_syg1)))
print('Mediana różnicy fazy:', np.median(r_fazy))
plt.plot(r_fazy)
plt.xlabel('nr próbki')
plt.ylabel('różnica fazy')
plt.title('Różnica fazy sygnałów')
plt.show()

Mediana różnicy fazy: 270.0071211520213

Obliczona różnica fazy wynosi 270 stopni, co uwzględniając cykliczność kąta fazowego, jest toższame z oczekiwaną wartością -90 stopni.

Filtry adaptacyjne

Redukcja szumów i zakłóceń w sygnale jest jednym z najczęstszych praktycznych zastosowań przetwarzania sygnałów. Szumy i zakłócenia są sygnałami losowymi (stochastycznymi), bardzo często również niestacjonarnymi (o zmiennych parametrach). Nie da się usunąć zakłóceń z sygnału za pomocą klasycznego filtru o stałej charakterystyce. Potrzebny jest filtr, który będzie dostosowywał swoją charakterystykę do zmiennego charakteru zakłóceń.

Filtr adaptacyjny (adptive filter) charakteryzuje się współczynnikami, których wartości nie są stałe, lecz dobierane przez algorytm. Cyfrowe filtry adaptacyjne są najczęściej filtrami o skończonej odpowiedzi impulsowej (FIR). Działanie filtru można opisać w następujący sposób. Sygnał $x$ jest poddawany filtracji, wynikiem jest sygnał $y$. Potrzebny jest również drugi sygnał odniesienia $d$, nazywany również sygnałem pożądanym (desired). Ważnym założeniem metody jest to, że sygnały $y$ i $d$ muszą być ze sobą skorelowane. Wynik filtracji jest odejmowany od sygnału odniesienia, powstaje w ten sposób sygnał błędu $e$. Zadaniem algorytmu jest minimalizacja sygnału błędu poprzez modyfikację współczynników filtru $h$.

Najczęściej stosowanym algorytmem adaptacji współczynników filtru jest algorytm minimalizacji błędu średniokowadratowego, LMS (least mean squares). Aktualizacja wartości współczynników filtru zachodzi zgodnie z zależnością:

$$h_k(n+1) = h_k(n) + \mu e(n) x(n-k) = h_k(n) + \mu \left(d(n) - y(n) \right) x(n-k)$$

Parametr $\mu$ nazywa się krokiem adaptacji. Zbyt duża wartość parametru nie pozwoli na znalezienie minimum sygnału błędu (nastąpi "przeskoczenie" minimum) i spowoduje zniekształcenie sygnału. Zbyt mała wartość parametru spowoduje zbyt wolną adaptację. Ponieważ algorytm LMS jest bardzo czuły na dobór wartości $\mu$, często stosuje się algorytm znormalizowany NLMS, w którym wartość dodawaną do współczynników filtru dzieli się przez sumę kwadratów wartości sygnału wejściowego.

Nie ma gotowej funkcji realizującej filtr LMS w module SciPy, można jednak łatwo ją napisać samodzielnie. Poniższa funkcja realizuje filtr LMS, lub NLMS gdy parametr norm jest ustawiony na True. Próbki sygnału $x$ są zapisywane w buforze kołowym, funkcja numpy.roll przesuwa cyklicznie zawartość tablicy w prawo. Na miejsce najstarszej wartości filtrowanego sygnału wpisywana jest nowa wartość. Filtracja jest wykonywana za pomocą funkcji numpy.dot. Funkcja zwraca sygnał z wyjścia filtru, sygnał błędu oraz końcowe współczynniki filtru adaptacyjnego.

In [36]:
def filtr_adaptacyjny(sygnal_x, sygnal_d, n, krok=0.01, norm=False):
    '''Implementacja filtru adaptacyjnego LMS.
    Parametry:
        sygnal_x: sygnał poddawany filtracji adaptacyjnej
        sygnal_d: sygnał referencyjny
        n: liczba współczynników filtru adaptacyjnego
        krok: krok aktualizacji współczynników
        norm: algorytm LMS (False) lub NLMS (True)
    Zwracane wartości:
        sygnal_y: wartości sygnału po filtracji adaptacyjnej
        sygnal_e: wartości sygnału błędu
        wsp_filtru: końcowe współczynniki filtru adaptacyjnego
    '''
    sygnal_y = []
    sygnal_e = []
    # początkowe wartości współczynników
    b = np.zeros(n)
    b[0] = 1
    bufor = np.zeros(n)
    for x, d in zip(sygnal_x, sygnal_d):
        # wpisanie nowej próbki do bufora kołowego
        bufor = np.roll(bufor, 1)
        bufor[0] = x
        # filtracja adaptacyjna
        y = np.dot(b, bufor)
        sygnal_y.append(y)
        # sygnał błędu
        e = d - y
        sygnal_e.append(e)
        # adaptacja filtru
        u = krok * e * bufor
        if norm:
            u /= (1 + np.dot(bufor, bufor))
        b += u
    return sygnal_y, sygnal_e, b

Niech sygnałem poddawanym filtracji będzie sygnał sinusoidalny o częstotliwości 1000 Hz, do którego zostało dodane zakłócenie w formie sygnału sinusoidalnego o częstotliwości 50 Hz. Zakładamy, że nie znamy sygnału zakłócającego. Postać czasowa zniekształconego sygnału jest następująca.

In [37]:
x = np.sin(2 * np.pi * np.arange(1000) * 1000 / fs)
n = np.sin(2 * np.pi * np.arange(1000) * 50 / fs)
xn = 0.9 * x + 0.1 * n
plt.plot(xn)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał z zakłóceniem')
plt.show()

Sygnał odniesienia musi być skorelowany z sygnałem, który chcemy uzyskać z wyjścia filtru adaptacyjnego, czyli z sygnałem bez zakłóceń. Niech sygnałem tym będzie kosinus o częstotliwości rownież 1000 Hz. Stosujemy filtr 30. rzędu, wartość kroku $\mu$ = 0,005. Poniższy wykres pokazuje sygnał y z wyjścia filtru oraz sygnał błędu. Widać, że filtr adaptacyjny jest w stanie dość szybko usunąć zakłócenie z sygnału.

In [38]:
d = np.cos(2 * np.pi * np.arange(1000) * 1000 / fs)
y, e, wsp = filtr_adaptacyjny(xn, d, 31, 0.005)
plt.plot(y, label='y(n)')
plt.plot(e, label='e(n)')
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał po filtracji adaptacyjnej')
plt.legend(loc='upper right')
plt.show()

W praktycznych sytuacjach nie mamy jednak dostępnego sygnału pożądanego, który można byłoby użyć do porównania z wynikiem filtracji. Możemy jednak często pobrać próbkę zakłóceń, np. z fragmentu, w którym nie występuje sygnał użyteczny. W takim przypadku zmienia się struktura algorytmu. Na wejście filtru adaptacyjnego podajemy próbkę sygnału zakłócającego, która musi być skorelowana z sygnałem zakłócającym obecnym w sygnale. Wynikiem filtracji jest estymata zakłócenia występującego w sygnale, sygnałem odniesienia d jest natomiast sygnał zakłócony. Minimalizacja sygnału błędu wpływa tylko na wynik filtracji (estymatę zakłóceń), nie minimalizuje natomiast sygnału użytecznego (nie jest on podawany na wejście filtru adaptacyjnego). Sygnałem, który chcemy uzyskać (bez zakłóceń) jest w tym przypadku sygnał $e$ (błędu).

W poniższym przykładzie, do sygnału sinusoidalnego o częstotliwości 500 Hz dodawany jest szum biały, przetworzony przez układ o nieznanej transmitancji (symulujemy go filtrem FIR o współczynnikach 0,1; 0,3; 0,2). Próbką zakłóceń jest szum biały, jest on podawany na wejście filtru adaptacyjnego. Oczekujemy, że na wyjściu filtru uzyskamy estymatę przefiltrowanego szumu, tak że po odjęciu od sygnału odniesienia, sygnałem błędu będzie sygnał wejściowy bez zakłóceń. Funkcja realizująca fitrację adaptacyjną jest ta sama co poprzednio, podaje się tylko inne sygnały jako argumenty. Dodatkowo, w tym przypadku został użyty algorytm NLMS. Wykorzystano filtr 10. rzędu (wyższe rzędy nie poprawiały skuteczności), współczynnik $\mu$ = 0,1.

In [39]:
np.random.seed(9000)
x = np.sin(2 * np.pi * np.arange(1000) * 500 / fs)
n = np.random.random(1000) - 0.5
nf = sig.lfilter([0.1, 0.3, 0.2], 1, n)
xn = x + nf
y, e, wsp = filtr_adaptacyjny(n, xn, 11, 0.1, True)
plt.plot(y, label='y(n)')
plt.plot(e, label='e(n)')
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał po filtracji adaptacyjnej')
plt.legend(loc='upper right')
plt.show()
print(wsp)

[ 0.0815132   0.25957684  0.1638007  -0.04153777 -0.04482144 -0.02375405
 -0.02544225 -0.03232618 -0.04666007 -0.04186204 -0.04978637]

Na wykresie widać działanie filtru adaptacyjnego, z sygnału $e$ zostały częściowo usunięte zakłócenia. Efekt nie jest jednak idealny, filtracja adaptacyjna w przypadku szumu jest problematyczna ze względu na niski stopień korelacji próbki szumu z rzeczywistym zakłóceniem. Obliczone współczynniki filtru są dość bliskie użytym w symulacji (0,1; 0,3; 0,2), widać więc, że filtr przybliża rzeczywistą charakterystykę zakłóceń.

Powyższe przykłady są bardzo proste, stanowią one jedynie ilustrację działania metody. Filtry adaptacyjne są jednak szeroko stosowane do redukcji szumu, usuwania echa telekomunikacyjnego i akustycznego, identyfikacji systemów i w wielu innych zastosowaniach.

Filtr predykcyjny

Predykcja liniowa (linear prediction) polega na obliczeniu kolejnej wartości sygnału jako liniowej kombinacji poprzednich wartości. Mówiąc prościej, wartość sygnału jest równa sumie poprzenich wartości sygnału przemnożonych przez współczynniki predykcji. Zatem filtr predykcyjny (linear predicition filter) może być zrealizowany jako filtr FIR, którego współczynnik $b_0$ jest zerowy:

$$\hat{x}(n) = b_1 x(n-1) +b_2 x(n-2) + ... + b_N x(n-N)$$

Jako przykład zostanie wykorzystany sygnał, który był użyty wcześniej do ilustracji operacji wygładzania sygnału. Sygnał po wygładzeniu ma następującą postać.

In [46]:
x = np.loadtxt('pliki/noisy_data.txt')
x = sig.savgol_filter(x, 31, 3)
plt.plot(x)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał do predykcji')
plt.show()

Do obliczenia współczynników filtru predykcyjnego używa się zwykle funkcji autokorelacji sygnału. Sposób obliczania tej funkcji pokazano w rozdziale Splot i korelacja. Za pomocą funkcji scipy.signal.correlate obliczana jest funkcja autokorelacji sygnału z odjętą składową stałą. Funkcja autokorelacji jest następnie normalizowana do jedności dla zerowego przesunięcia, poprzez podzielenie wartości funkcji przez sumę kwadratów funkcji, a następnie wybierana jest część funkcji od zerowego przesunięcia.

In [47]:
xnorm = x - np.mean(x)
autokor = sig.correlate(xnorm, xnorm)[len(x) - 1:] / np.sum(xnorm * xnorm)
plt.plot(autokor)
plt.xlabel('opóźnienie')
plt.ylabel('autokorelacja')
plt.title('Funkcja autokorelacji sygnału do predykcji')
plt.show()

Następnie zakłada się rząd predykcji $N$ i oblicza się współczynniki za pomocą następującej zależności:

$$\textbf{R} \cdot \textbf{b} = -\textbf{x}$$

gdzie:

  • $\textbf{r}$ jest wektorem wartości funkcji autokorelacji,
  • $\textbf{R}$ jest macierzą autokorelacji, w której $\textbf{R}_{i,j} = \textbf{r}_{|i-j|}$,
  • $\textbf{b}$ jest szukanym wektorem współczynników predykcji,
  • $\textbf{x}$ jest wektorem wartości funkcji autokorelacji z pominięciem "zerowej" wartości: $\textbf{x}_i = \textbf{r}_{i+1}$,
  • $i, j = 0, ..., N$.

Rozwiązanie powyższego równania macierzowego jest wykonywane metodą Levinsona-Durbina. Służy do tego funkcja solve_toeplitz z modułu scipy.linalg. Pierwszym argumentem jest wektor, z którego tworzona jest macierz $\textbf{R}$, nazywana macierzą Toeplitza, drugi argument to wektor $\textbf{x}$. W poniższym przykładzie, zakładamy rząd predykcji równy 30, zatem chcemy obliczyć 31 współczynników filtru. Funkcja zwraca obliczony wektor $-\textbf{b}$, zatem odwracamy jego znak. Zerowy współczynnik jest następnie zerowany, a pozostałe współczynniki są normalizowane, aby zachować jednostkowe wzmocnienie filtru.

In [49]:
L = 31
b = -linalg.solve_toeplitz(autokor[:L], autokor[1:L+1])
b[0] = 0
b /= np.abs(np.sum(b))
print(b)

[ 0.          0.14850395  0.14357786  0.12649404  0.09115896  0.08377893
  0.06876021  0.06687199  0.05807353  0.04640319  0.03889321  0.02657993
  0.02235171  0.02169929  0.01970002 -0.04058587  0.04097344  0.01579166
 -0.0091857   0.06046151 -0.01035889  0.00874652  0.0104003  -0.04305643
 -0.04354969 -0.01129355 -0.02985861  0.00881199  0.0622828   0.05136301
 -0.03378929]

Poniższy wykres pokazuje wynik filtracji sygnału przez filtr predykcyjny. Widać, że predykowany sygnał podąża za przebiegiem oryginalnego sygnału, jednak "spóźnia się" za nim o kilka próbek.

In [51]:
x_pred = sig.lfilter(b, 1, x)
plt.plot(x, label='oryginalny')
plt.plot(x_pred, label='po predykcji')
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał oryginalny i po predykcji')
plt.legend()
plt.show()

Różnica między oryginalnym sygnałem a sygnałem po predykcji nazywa się sygnałem błędu predykcji (prediction error) i oznacza się go jako $e(n)$. W praktyce nie oblicza się sygnału $\hat{x}(n)$, lecz od razu oblicza się błąd predykcji. Filtr predykcyjny do obliczania sygnału $e(n)$ ma współczynniki $1, -b_1, -b_2, ..., -b_N$, gdzie współczynniki $b_i$ są takie, jak obliczone wcześniej dla filtru predykcyjnego.

In [52]:
be = -b
be[0] = 1
e = sig.lfilter(be, 1, x)
plt.plot(e)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał błędu predykcji')
plt.show()

Wartości błędu predykcji są dla tego sygnału duże z powodu pokazanego wcześniej opóźnienia predykowanego sygnału względem oryginalnego. Gdyby wyrównać oba sygnały, błąd predykcji byłby znacznie mniejszy.

Mając sygnał błędu predykcji oraz współczynniki filtru predykcyjnego, można odtworzyć oryginalny sygnał. Filtr rekonstruujący sygnał, również nazywany filtrem predykcyjnym, jest odwrotnością filtru obliczającego błąd predykcji. Jest to filtr IIR zawierający wyłącznie bieguny transmitancji (all-pole filter).

In [54]:
x_rek = sig.lfilter([1], be, e)
plt.plot(x_rek)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał odtworzony z predykcji')
plt.show()
print('Błąd średniokwadratowy:', np.sqrt(np.sum((x - x_rek)**2)))

Błąd średniokwadratowy: 1.5550318294528701e-13

Nie wykreślono oryginalnego sygnału, ponieważ oba przebiegi pokrywają się na wykresie. Można więc odtworzyć sygnał mając współczynniki filtru predykcyjnego oraz wartości sygnału błędu.

Predykcja liniowa jest szeroko stosowana w kodowaniu oraz bezstratnej kompresji danych. Przy właściwym przeprowadzeniu predykcji, wartości sygnału błędu predykcji są małe, a dzięki temu można je wydajnie kodować. Wystarczy więc zapisać lub przesłać wartości błędu predykcji, a dekoder może odtworzyć sygnał znając współczynniki filtru predykcyjnego. Przykładowe zastosowania predykcji to kodowanie dźwięku metodą ADPCM oraz kompresja obrazu w standardzie PNG. Predykcja liniowa jest stosowana w liniowym kodowaniu predykcyjnym (LPC, linear predictive coding), które jest ważnym elementem systemów kodowania mowy w systemach telekomunikacyjnych, w tym telefonii GSM (np. kodeki CELP).

Podsumowanie

  • Wygładzanie sygnału polega na zredukowaniu szumu obecnego w sygnale. Można wykorzystać filtr średniej ruchomej (MA) lub filtr Sawickiego-Golaya. Większe rzędy filtrów mocniej tłumią szum, ale powodują większą utratę szczegółów sygnału.
  • Filtr różniczkujacy może być zrealizowany jako filtr różnicowy pierwszego lub drugiego rzędu. Jednak aby uzyskać poprawną charakterystykę częstotliwościową, należy użyć metody optymalnego projektowania filtru FIR o asymetrycznej odpowiedzi impulowej, np. metody Parksa-McClellana.
  • Układ całkujący sygnał może zostać zrealizowany jako filtr IIR pierwszego rzędu. Aby uniknąć niestabilności filtru, wprowadza się współczynnik straty, nieznacznie mniejszy od jedności.
  • Filtr tłumiący składową stałą sygnału można zrealizować łącząc kaskadowo filtr różniczkujący z filtrem całkującym.
  • Filtr grzebieniowy ma charakterystykę powtarzającą się cyklicznie. Może on być nierekursywny lub rekursywny. Filtr grzebieniowy może być wykorzystany jako filtr środkowozaporowy lub środkowoprzepustowy o wielu pasmach zaporowych/przepustowych rozłożonych równomiernie na skali częstotliwości.
  • Filtr wszechprzepustowy powstaje przez kaskadowe połączenie filtrów grzebieniowych: nierekursywnego i rekursywnego. Powoduje on przesunięcie fazowe sygnału, nie modyfikując jego charakterystyki amplitudowej. Jest przydatny w korekcji fazowej sygnałów, np. po filtracji IIR.
  • Filtr i transformator Hilberta tworzą sygnał analityczny - zespolony sygnał, którego widmo amplitudowe jest zerowe dla ujemnych częstotliwości. Sygnał analityczny jest użyteczny np. przy wykrywaniu przesunięć fazowych czy demodulacji sygnałów.
  • Filtry adaptacyjne to takie, których współczynniki są obliczane na podstawie różnicy sygnału odniesienia i sygnału przefiltrowanego, tak aby zminimalizować tę różnicę. Do adaptacji filtru stosuje się zwykle algorytm LMS i jego modyfikacje. Filtry adaptacyjne są użyteczne m.in. w redukcji zakłóceń w sygnale.
  • Predykcja liniowa polega na obliczeniu kolejnej wartości sygnału jako liniowej kombinacji poprzednich wartości. Filtr predykcyjny oblicza sygnał błędu predykcji. Znając współczynniki filtru predykcyjnego, można z sygnału błedu odzyskać oryginalny sygnał. Predykcja liniowa ma zastosowanie m.in w kodowaniu i bezstratnej kompresji dźwięku (zwłaszcza mowy) i obrazu.