Akademia Górniczo-Hutnicza

Katedra Robotyki i Mechatroniki

Identyfikacja i analiza sygnałów

Laboratorium 2

Wprowadzenie do przetwarzania sygnałów w dziedzinie

czasu

Przetwarzanie sygnału dla potrzeb badania dynamiki układów
mechanicznych

Problemy wstępnego przetwarzania sygnału w dziedzinie czasu

Podstawową czynnością w przetwarzaniu sygnałów jest przetwarzanie sygnału o naturze
analogowej do postaci cyfrowej. Proces ten polega na zamianie wielkości analogowej na ciąg
dyskretnych wartości w zadanych chwilach czasowych. Operacja ta nie wnosi nowych
informacji o sygnale, czasem może prowadzić do utraty pewnej ilości informacji.
Przetwarzanie analogowo-cyfrowe sygnałów składa się z dwóch podstawowych czynności:

−

kwantowanie sygnału,

−

dyskretyzacja sygnału.

0

1

2

3

4

5

6

7

8

n

∆

t

t

1

2

3

4

5

6

7

(0,4) (1,7) (2,6) (3,2)

(4,1) (5,4) (6,7) (7,5)...(n,7)

X

Rys. 1 Schematyczne przedstawienie procesu przetwarzania analogowo

−

cyfrowego.

Dyskretyzacja sygnału polega na wyborze na osi czasu dyskretnych chwil czasowych, w
których określona jest wartość amplitudy sygnału w procesie przetwarzania. Dyskretyzacja
sygnału odbywa się z okresem próbkowania zależnym od wymaganych parametrów analizy,
zakresu częstości w analizowanym sygnale tak, aby spełnione było twierdzenie Shanona o
próbkowaniu sygnałów. Twierdzenie to mówi, że: "Sygnał nie zawierający częstości
większych niż f

N

(częstość graniczna) może być jednoznacznie opisany za pomocą

dyskretnych próbek, otrzymanych przez próbkowanie analizowanego sygnału z częstością

f

f

S

N

≥

2

" W praktyce częstość f

S

jest 4 do 10 razy większa niż częstość f

N

.

Kolejną czynnością jest kwantyzacja polegająca na zamianie wartości analogowej na postać
cyfrową o wartości całkowitej związanej z ustalonym poziomem kwantyzacji.. Przetwarzanie
analogowo-cyfrowe sygnałów jest wykonywane za pomocą specjalizowanych układów
elektronicznych.

Przykłady związane z doborem częstotliwości próbkowania

Odpowiedni dobór częstotliwości próbkowania:

f1 = 1; % Częstotliwości przykładowych przebiegów sinusoidalnych
f2 = 4;

f3 = 6;

fs = 200; % Częstotliwość próbkowania sygnału
t=0:(1/fs):1; % Wektor czasu

% Definicje przykładowych przebiegów sinusoidalnych

x1 = sin(2*pi*f1*t);
x4 = -sin(2*pi*f2*t);

x6 = sin(2*pi*f3*t);
plot(t,x1,t,x4,t,x6)

Wybór zbyt niskiej częstotliwości próbkowania:

fs = 5; % Częstotliwość próbkowania sygnału
t=0:(1/fs):1; % Wektor czasu

% Definicje przykładowych przebiegów sinusoidalnych

x1 = sin(2*pi*f1*t);
x4 = -sin(2*pi*f2*t);

x6 = sin(2*pi*f3*t);
plot(t,x1,t,x4,t,x6)

Kolejny przykład złego doboru częstotliwości próbkowania zawarty jest poniżej.

fs = 20; % Częstotliwość próbkowania sygnału

t=0:(1/fs):10; % Wektor czasu

x = sin(2*pi*40*t);
plot(t,x)

Przedstawiona w powyższych przykładach niejednoznaczność przetwarzania analogowo
cyfrowego jest podstawowym problemem przetwarzania sygnałów. Wynika ona z faktu, że w
przypadku gdy w sygnale znajdują się składowe o częstości większej niż f

S

/ 2 wówczas nie

mogą one być odtworzone w sposób jednoznaczny, a mogą mieć wpływ na wartość amplitudy
sygnału. Zjawisko to jest nazywane aliazingiem.
Aby uniknąć zjawiska aliazingu należy przed przetwarzaniem usunąć z sygnału wszystkie
składowe o częstości większej niż f

S

/ 2 . W tym celu stosuje się na etapie kondycjonowania

sygnału filtrację analogową za pomocą filtru dolnoprzepustowego, o częstości przepuszczania

mniejszej niż

f

S

2

. Ze względu na to, że nie istnieją filtry idealne o pionowym nachyleniu

zbocza, częstość graniczną filtru przyjmuje się w praktyce ok. 40 % f

S

. Nachylenie zbocza

filtru zależy od jego rzędu i dla najczęściej stosowanych filtrów typu Buttworta 6 rzędu
wynosi 48 dB na oktawę. Do bardzo efektywnych filtrów należą filtry Cautera, dla których
nachylenie zbocza wynosi 96 dB/oktawę.

Zastosowanie filtracji cyfrowej

W tej części przedstawiono możliwości toolboxu Signal Processing Toolbox (SPT) z zakresu
analizy i projektowania filtrów cyfrowych.

Splot i filtrowanie
Matematyczną podstawą filtracji jest splot. Matlabowska funkcja

conv

dokonuje

standardowego jednowymiarowego splotu dwóch wektorów.

conv([1 1 1],[1 1 1])

•

ans =

•

1 2 3 2 1

Splotu macierzy prostokątnych dla dwuwymiarowego przetwarzania sygnałów można
natomiast dokonać za pomocą funkcji

conv2

.

Sekwencja danych y(n) otrzymywana na wyjściu filtru cyfrowego jest splotem ciągu danych
wejściowych x(n) przez splot z odpowiedzią impulsowego filtru h(n):

∑

∞

− ∞

=

−

=

∗

=

m

)

m

(

x

)

m

n

(

h

)

n

)(

x

h

(

)

n

(

y

( 1 )

i może być interpretowana jako ruchoma średnia ważona wejścia.
Jeżeli odpowiedź impulsowa filtru cyfrowego h(n) oraz sygnał wejściowy x(n) mają
skończone długości, można wtedy dokonać filtracji za pomocą funkcji

conv

. Przykładowo

utworzono sygnał wejściowy x(n) jako wektor x, sygnał h(n) jako wektor h, a następnie
dokonano ich splotu.

x=randn(5,1);

h=[1 1 1 1]/4;
y=conv(h,x);

Funkcja przejścia filtru
Definicja funkcji przejścia filtru oparta jest na własnościach przekształcenia Z. Transformata
Z Y(z) z cyfrowego wyjścia filtru y(n) jest powiązana z transformacją Z X(z) wejścia x(n)
następującym wzorem:

)

z

(

X

z

)

1

na

(

a

...

z

)

2

(

a

1

z

)

1

nb

(

b

...

z

)

2

(

b

)

1

(

b

)

z

(

X

)

z

(

H

)

z

(

Y

na

1

nb

1

⋅

⋅

+

+

+

⋅

+

⋅

+

+

+

⋅

+

=

⋅

=

−

−

−

−

( 2 )

gdzie H(z) jest funkcją przejścia filtru. Stałe b(i) oraz a(i) są współczynnikami filtru, a
wartość na i nb reprezentuje rząd filtru. Współczynniki filtru są zebrane w dwóch wektorach:
b - dla licznika i a - dla mianownika. Według przyjętej konwencji Matlab używa wektora
wierszowego dla zapamiętania współczynników.

Uwaga: Indeksowanie współczynników w wektorach rozpoczyna się od liczby 1, zamiast od
0. Jest to standardowy schemat stosowanego w Matlabie sposobu indeksowania wektorów.

Filtrowanie przy pomocy funkcji

filter

W celu przedstawienia sposobu działania funkcji

filter

cofnijmy się do zależności (2). Po

pomnożeniu równania (2) przez mianownik ułamka i dokonaniu odwrotnej transformaty Z
otrzymuje się następującą zależność:

)

nb

n

(

x

)

1

nb

(

b

)

1

n

(

x

)

2

(

b

...

)

n

(

x

)

1

(

b

)

na

n

(

y

)

1

na

(

a

...

)

1

n

(

y

)

2

(

a

)

n

(

y

−

⋅

+

+

+

−

⋅

+

+

⋅

=

−

⋅

+

+

+

−

⋅

+



( 3 )

Z równania (3) wynika, że wartość y(n) można wyznaczyć znając wartości wejściowe x(n),
x(n-1), ...,x(n-nb):

)

na

n

(

y

)

1

na

(

a

...

)

1

n

(

y

)

2

(

a

)

nb

n

(

x

)

1

nb

(

b

...

)

1

n

(

x

)

2

(

b

)

n

(

x

)

1

(

b

)

n

(

y

−

⋅

+

−

−

−

⋅

−

−

⋅

+

+

+

−

⋅

+

⋅

=

( 4 )

( )

)

k

n

(

y

)

1

k

(

a

)

m

n

(

x

)

1

m

(

b

n

y

na

0

k

nb

0

m

−

⋅

+

−

−

⋅

+

=

∑

∑

=

=

Równanie (4) jest standardową reprezentacją filtru cyfrowego w dziedzinie czasu i nosi
nazwę tzw. równania różnicowego. Obliczenia rozpoczyna się dla y(1), zakładając zerowe
warunki początkowe. Algorytm obliczeń przedstawiają równania (5):

⋮

)

1

(

y

)

3

(

a

)

2

(

y

)

2

(

a

)

1

(

x

)

3

(

b

)

2

(

x

)

2

(

b

)

3

(

x

)

1

(

b

)

3

(

y

)

1

(

y

)

2

(

a

)

1

(

x

)

2

(

b

)

2

(

x

)

1

(

b

)

2

(

y

)

1

(

x

)

1

(

b

)

1

(

y

−

−

+

+

=

−

+

=

=

( 5 )

W tej formie filtr cyfrowy jest wprowadzany przez funkcję

filter

. Dla przykładu, prosty

rekurencyjny filtr dolnoprzepustowy z pojedynczym biegunem (zero wielomianu mianownika
funkcji przejścia H(z) (2)) jest dany jako:

b=1;

a=[1 -0.9];

gdzie wektory a i b reprezentują współczynniki filtru w formie funkcji przejścia. Filtr
realizowany jest poprzez wyrażenie:

y=filter(b,a,x);

Funkcja

filter

oblicza tyle próbek wyjściowych, ile jest próbek wejściowych, znaczy to, że

długość wektora y jest taka sama, jak długość wejściowego wektora x. Jeśli pierwszy element
a(1) nie jest 1, funkcja

filter

dzieli przez niego współczynniki uprzednio wyprowadzonych

równań różnicowych postaci (4).

Odpowiedź impulsowa
Odpowiedź impulsową (charakterystykę czasową) filtru cyfrowego otrzymujemy podając na
wejście filtru następujący ciąg próbek:









≠

=

=

1

n

dla

0

1

n

dla

1

)

n

(

x

W Matlabie funkcję impulsową, można wygenerować na wiele sposobów. Najprostszy to
wykonanie poniższego polecenia:

imp=[1 zeros(1,49)];

Odpowiedź impulsową prostego filtru o współczynnikach

b=1;

a=[1 -0.9];

otrzymamy

po wykonaniu polecenia:

h=filter(b,a,imp);

Funkcja

impz

w przyborniku upraszcza powyższą operację wybierając liczbę generowanych

punktów oraz tworząc wykres słupkowy odpowiedzi filtru (używając wewnętrznie funkcji

stem

):

impz(b,a)

0

20

40

60

80

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Rys. 2 Odpowiedź impulsowa przykładowego filtru uzyskana za pomocą funkcji impz

Wykres pokazuje eksponencjalny rozkład 0.9

n

układu z pojedynczym biegunem.

Funkcja filter
Funkcja

filter

jest wprowadzana jako transpozycja struktury przedstawionej na rys. 3,

gdzie n-1 jest rzędem filtru.

Σ

Σ

Σ

Σ

z

n-1

(m)

z

-1

z

-1

z

-1

z

2

(m)

z

1

(m)

...

...

...

b(n)

b(1)

b(2)

b(3)

-a(n)

y(m)

-a(2)

-a(3)

x(m)

Rys. 3 Struktura blokowa filtru wprowadzanego w SPT

Jest to podstawowa forma filtru cyfrowego, posiadająca minimalną ilość elementów
opóźniających.
Przy m - tej próbce sygnału, funkcja

filter

rozwiązuje następujący układ równań

różnicowych:

)

m

(

y

)

n

(

a

)

m

(

x

)

n

(

b

)

m

(

z

)

m

(

y

)

1

n

(

a

)

1

m

(

z

)

m

(

x

)

1

n

(

b

)

m

(

z

)

m

(

y

)

2

(

a

)

1

m

(

z

)

m

(

x

)

2

(

b

)

m

(

z

)

1

m

(

z

)

m

(

x

)

1

(

b

)

m

(

y

1

n

1

n

2

n

2

1

1

−

=

−

−

−

+

−

=

=

−

−

+

=

−

+

=

−

−

−

⋮

⋮

( 6 )

W najbardziej podstawowej formie, funkcja

filter

wprowadza opóźnienia z

i

(1), i=1,...,n-1,

równe zero. Opóźnienia te można jednak nastawiać wprowadzając czwarty parametr
wejściowy w funkcji. Można także uzyskać dostęp do opóźnienia wprowadzanego przez
funkcję poprzez wykorzystanie drugiego parametru wyjściowego:

[y,zf]=filter(b,a,x,zi);

Dostęp do wyżej wymienionych parametrów jest użyteczny podczas filtrowania danych w
sekcjach, zwłaszcza gdy należy się liczyć z ograniczeniami pamięci. Przykładowe dane
zgromadzone zostały w dwóch segmentach po 5000 punktów każdy:

x1=randn(5000,1);
x2=randn(5000,1);

Załóżmy, że pierwsza sekwencja x1 odpowiada pierwszym 10 minutom danych, druga zaś x2
kolejnym 10 minutom. Cała sekwencja ma długość x = [x1;x2]. Jeśli nie posiadamy
wystarczającej ilości pamięci do przetworzenia całej sekwencji od razu, przefiltrujmy
podsekwencje jedna po drugiej. Do zapewnienia ciągłości filtrowanych danych, użyjmy
warunków końcowych z filtracji x1 jako warunków początkowych w wywołaniu filtracji x2:

[y1,zf]=filter(b,a,x1);
y2=filter(b,a,x2,zf);

W przyborniku dostępna jest również dodatkowa funkcja

filtic

, generująca zewnętrznie

warunki początkowe dla funkcji

filter

. Utwórzmy taką samą wartość opóźnienia jak w

przykładzie powyżej używając do tego celu funkcji

filtic

:

zf=filtic(b,a,flipud(y1),flipud(x1));

Ta funkcja jest użyteczna podczas filtracji krótkich sekwencji danych, ponieważ odpowiednie
warunki początkowe eliminują startowe stany nieustalone.

Podstawowe funkcje stosowane do filtracji
Oprócz funkcji

filter

w SPT znajdują się jeszcze dwie inne instrukcje, które dokonują

podstawowych operacji filtrowania. Są to

filtfilt

eliminująca zniekształcenia fazowe w

procesie filtrowania, oraz

fftfilt

, która dokonuje filtracji FIR (Finite Impulse Response) w

dziedzinie częstotliwościowej.

Ogólne zasady filtracji cyfrowej w Matlabie
Implementacja filtru FIR o liniowej fazie może być zrealizowane dla funkcji

filter

lub

conv

poprzez dodanie do sekwencji danych opóźnienia wyjścia o stałą liczbę próbek. Dla

filtrów IIR (Infinite Impulse Response) natomiast zniekształcenia fazowe są stosunkowo
duże. Funkcja

filtfilt

eliminuje te zniekształcenia, także w przypadku filtrów typu FIR.

Prześledźmy sposób działania funkcji

filtfilt

. Oznaczmy transformatę Z sekwencji x(n)

jako X(z) oraz transformatę Z sekwencji x(n) odwróconej w czasie jako X(1/z). Opisana
funkcja realizuje obliczenia według schematu przedstawionego na rys. 4.

Odwrotność

czasu

Odwrotność

czasu

H(z)

H(z)

X(z)

X(z)H(z)

X(1/z)H(1/z)

X(1/z)H(1/z)H(z)

Y(z)=X(z)H(z)H(1/z)

Rys. 4 Schemat działania funkcji filtfilt

Gdy



z



=1, co znaczy że z=e

j

ω

, wyjście redukuje się do:

2

j

j

)

e

(

H

)

e

(

X

ω

ω

⋅

. Wszystkie próbki

ciągu wyjściowego y(n) stanowią podwójnie przefiltrowany ciąg danych wejściowych x(n).
Dla przykładu przefiltrowano próbkę sygnału o długości 1 sekundy, próbkowanego z
częstotliwością 100 Hz, składającego się z dwóch sinusoidalnych składników, jednego o
częstotliwości 3 i drugiego o częstotliwości 40 Hz.

Fs=100;
t=0:1/Fs:1;

x=sin(2*pi*t*3)+0.25*sin(2*pi*t*40);

Utworzono również prosty uśredniający filtr FIR dziesiątego rzędu, a następnie
przefiltrowano sygnał

x

używając w tym celu obu funkcji

filtfilt

i

filter

.

b=ones(1,10)/10;
y=filtfilt(b,1,x);

yy=filter(b,1,x);
subplot(1,1,1);

plot(t,x,'black',t,y,'--black',t,yy,':black')

0

0 . 1

0 . 2

0 . 3

0 . 4

0 . 5

0 . 6

0 . 7

0 . 8

0 . 9

1

- 1 . 5

- 1

- 0 . 5

0

0 . 5

1

1 . 5

P r z e b i e g w e j c i o w y

ś

F i l t r a c j a f u n k c j f i l t f i l t

ą

F i l t r a c j a f u n k c j f i l t e r

ą

Rys. 5 Porównanie wyników filtracji sygnału dokonanej przez funkcje filtfilt i filter

Powyższy wykres pokazuje różnice pomiędzy zastosowaniem funkcji

filter

i

filtfilt

.

Obydwa filtry usunęły sinusoidę o częstotliwości 40 Hz z oryginalnego sygnału. Różnice
pomiędzy wykresami dotyczą zarówno ich amplitudy jak i fazy. Wynik działania funkcji

filtfilt

jest dokładnie w fazie z oryginalnym sygnałem o częstotliwości 3 Hz, podczas

gdy wynik funkcji

filter

jest przesunięty o około pięć próbek. W przypadku funkcji

filtfilt

także amplituda sygnału wyjściowego jest mniejsza, a to ze względu na kwadrat

modułu funkcji przenoszenia.
Funkcja

filtfilt

zmniejsza również wpływ stanów nieustalonych. Dla uzyskania

poprawnych wyników należy się upewnić, czy sekwencja przeznaczona do filtrowania,
posiada co najmniej trzykrotną długość rzędu filtru oraz czy kończy się zerami przy obydwu
brzegach.

Specyfikacja filtrów
Celem projektowania filtrów jest dokonanie zmian w zawartości składników
częstotliwościowych ciągu danych. Typowym wymaganiem może być np. usunięcie z
sygnału szumów o częstotliwościach powyżej 30 Hz z sekwencji danych próbkowanych z
częstotliwością 100 Hz. Ponadto w specyfikacji filtrów muszą być podane wymagania co do
szerokości pasma przepuszczania, tłumienia w paśmie zaporowym, dopuszczalnych
zniekształceń amplitudy lub fazy.
Dla wymagań, jak w przypadku podanym wyżej, wystarczający jest filtr typu Butterwortha
(IIR). Dla przykładu zaprojektowano filtr dolnoprzepustowy piątego rzędu, o częstotliwości
granicznej 30 Hz, dla danych wejściowych znajdujących się w wektorze

x

. Zaprojektowania

filtru i filtracji sygnału dokonano używając następujących poleceń:

[b,a]=butter(5,30/50);

y=filter(b,a,x);

Drugi argument wejściowy w funkcji

butter

wskazuje częstotliwość odcięcia filtru, która

jest normalizowana do połowy częstotliwości próbkowania (częstotliwość Nyquista).

Uwaga: Wszystkie funkcje Matlaba wykorzystujące filtry działają z normalizowaną
częstotliwością, więc nie jest wymagane podawanie częstotliwości próbkowania jako
dodatkowego argumentu. W przyborniku założono częstotliwości Nyquista jako jednostkową.
Normalizowana częstotliwość jest więc zawsze zawarta w przedziale 0

≤

f

≤

1. Dla układu z

częstotliwością próbkowania 1000 Hz, częstotliwość 300 Hz wynosi 300/500=0.6. W celu
przeliczenia znormalizowanej częstotliwości na częstość kątową (rad/s) mnożymy ją przez

π

.

By przeliczyć częstotliwość znormalizowaną na częstotliwość w [Hz], mnożymy ją przez
połowę częstotliwości próbkowania.

Najbardziej rygorystyczne wymagania stawiane filtrom zawierają takie parametry jak
dopuszczalne oscylacje w paśmie przepuszczania (Rp, w decybelach), tłumienie w paśmie
zaporowym (Rs, w decybelach), oraz stromość zbocza (Ws - Wp, w Hz)

0

Wp

Ws

10

10

-Rs/20

-Rp/20

A

m

pl

it

ud

a

fi

lt

ru

Częstotliwość

Rys. 6 Typowa charakterystyka filtru

Można zaprojektować filtry: Butterwortha, Chebyshewa typu I i typu II oraz filtr eliptyczny
spełniające te typy wymagań. Istnieją także funkcje w przyborniku estymujące minimalną
klasę filtru spełniającego żądane warunki.
Do wymagań zawierających bardziej rygorystyczne ograniczenia, jak np. liniowość fazy lub
stromość zbocza filtru, powinno stosować się filtry typu FIR lub typu IIR projektowane
bezpośrednio.

Projektowanie filtrów typu IIR
Zaletą filtrów IIR w stosunku do filtrów FIR jest to, że zapewniają większą stromość zbocza
przy mniejszym rzędzie filtru. Filtry IIR posiadają jednak nieliniową fazę, co jest zjawiskiem
negatywnym, jednak przetwarzanie danych w Matlabie jest realizowane głównie off line (tzn.
cała sekwencja dostępna jest przed rozpoczęciem operacji filtrowania). Sytuacja taka pozwala
na zastosowanie specjalnego rodzaju filtracji (funkcja

filtfilt

), w celu wyeliminowania

nieliniowych zniekształceń fazowych filtru IIR.

Klasyczne filtry IIR, a więc Butterwortha, Chebyshewa typu I i typu II, eliptyczny oraz filtr
Bessela, aproksymują w różny sposób idealny filtr „prostokątny”. Przybornik SPT zawiera
funkcje do tworzenia wszystkich tych typów klasycznych filtrów IIR w dziedzinach:
analogowej (ciągłej) i cyfrowej (dyskretnej) (z wyjątkiem filtru Bessela, który można
zrealizować jedynie w dziedzinie ciągłej). Przybornik umożliwia także skonfigurowanie filtru
jako dolnoprzepustowego, górnoprzepustowego, środkowoprzepustowego, lub

środkowozaporowego. Dla głównych typów filtrów można także znaleźć ich najmniejszy
rząd, spełniający zadane wymagania w części związanej z przepuszczaniem lub tłumieniem
oraz stromością zboczy.
Za pomocą funkcji

yulewalk

projektuje się filtry, których odpowiedź amplitudowa jest

aproksymacją zadanej funkcji. Pozwala ona więc na tworzenie filtrów pasmowych
niekoniecznie „prostokątnych”.

Tabela 1 Opisy konstrukcji filtrów IIR dostępnych w przyborniku.

Technika

Opis

Funkcja

Prototypo-
wanie
analogowe

Uzyskanie filtru
cyfrowego na podstawie
jego analogowego
dolnoprzepustowego
prototypu
unormowanego,
podanego w formie
transmitancji iloczynowej
modelu ciągłego
(Laplace’a), poprzez
transformację
częstotliwości i
dyskretyzację filtru.

Kompletne funkcje
projektujące:

bessel

,

butter

,

cheby1

,

cheby2

,

ellip

Klasa estymacji funkcji:

buttord

,

cheb1ord

,

cheb2ord

,

ellipord

Dolnoprzepustowe,
unormowane, analogowe
prototypy funkcji:

lp2lp

,

lp2hp

,

lp2bp

,

lp2bs

Funkcje dyskretyzujące filtry:

bilinear

,

impinvar

Projektowanie
bezpośrednie

Projektowanie filtrów
cyfrowych bezpośrednio
w dziedzinie dyskretnej
poprzez przybliżanie
charakterystyki
amplitudowej odcinkami
liniowymi.

youlewalk

Modelowanie
parametryczne

Projektowanie filtrów
cyfrowych, których
charakterystyka w
dziedzinie czasu lub
częstotliwości jest
przybliżeniem
charakterystyki
założonej.

Funkcje modelujące w
dziedzinie czasu:

lpc

,

prony

,

stmcb

Funkcje modelujące w
dziedzinie częstotliwości:

invfreqz

,

invfreqs

Do projektowania filtrów IIR można także zastosować metody modelowania
parametrycznego lub metody identyfikacji.

Filtry IIR projektowane poprzez prototypy analogowe
Podstawowa technika projektowania cyfrowych filtrów IIR wprowadzana w przyborniku
bazuje na przekształceniu klasycznego unormowanego analogowego filtru
dolnoprzepustowego na jego cyfrowy odpowiednik. Ten rozdział przedstawia metodę
projektowania oraz główne charakterystyki podstawowych typów filtrów.

Filtr dowolnego rzędu w konfiguracji dolnoprzepustowej, górnoprzepustowej,
środkowoprzepustowej lub środkowozaporowej, można łatwo utworzyć za pomocą funkcji
zestawionych w tabelce 2.
Standardowo każda z tych funkcji oblicza parametry filtru dolnoprzepustowego, należy
jedynie podać wymaganą częstotliwość odcięcia Wn w znormalizowanej częstotliwości
(częstotliwość Nyquista = 1 Hz). Dla filtru górnoprzepustowego należy dołączyć parametr
‘high’. Dla filtrów środkowozaporowych i środkowoprzepustowych należy zadać Wn jako
dwuelementowy wektor, zawierający częstotliwości krawędzi pasma przepuszczania. Dla
konfiguracji środkowozaporowej dodatkowo należy dołączyć parametr ‘stop’.

Tabela 2 Zestawienie funkcji do projektowania filtru IIR

Filtr

Funkcje projektujące

Butterwortha

[b,a]=butter(n,Wn,opcje)

[z,p,k]= butter(n,Wn,opcje)
[A,B,C,D]= butter(n,Wn,opcje)

Chebyshewa typ I

[b,a]=cheby1(n,Rp,Wn,opcje)
[z,p,k]= cheby1(n,Rp,Wn,opcje)

[A,B,C,D]= cheby1(n,Rp,Wn,opcje)

Chebyshewa typ II

[b,a]=cheby2(n,Rs,Wn,opcje)

[z,p,k]= cheby2(n,Rs,Wn,opcje)
[A,B,C,D]= cheby2(n,Rs,Wn,opcje)

Eliptyczny

[b,a]=ellip(n,Rp,Rs,Wn,opcje)
[z,p,k]= ellip (n,Rp,Rs,Wn,opcje)

[A,B,C,D]= ellip (n,Rp,Rs,Wn,opcje)

Bessela (jedynie
analogowy)

[b,a]=besself(n,Wn,opcje)

[z,p,k]= besself(n,Wn,opcje)
[A,B,C,D]= besself(n,Wn,opcje)

Przykłady tworzenia niektórych filtrów cyfrowych:

[b,a]=butter(5,0.4); - dolnoprzepustowy filtr Butterwortha

[b,a]=cheby1(4,1,[0.4 0.7]); - środkowoprzepustowy filtr Chebyshewa typu

I

[b,a]=cheby2(6,60,0.8,’high’); - górnoprzepustowy filtr Chebyshewa typu

II

[b,a]=ellip(4,1,60,[0.4,0.7],'stop');

- środkowozaporowy filtr eliptyczny

Do skonstruowania filtru analogowego (np. w celu jego symulacji) należy jako ostatni
dołączyć parametr ‘

s

’ i podać częstotliwość odcięcia w radianach na sekundę

[b,a]=butter(5,0.4,’s’); - analogowy filtr Butterwortha
Wszystkie przedstawione tu funkcje umożliwiają wyznaczenie parametrów równania filtru w
postaci funkcji przejścia, transmitancji iloczynowej lub przestrzeni stanu, zależnie od ilości
zadawanych argumentów wyjściowych.

Przybornik zawiera również funkcje obliczające minimalny rząd filtru, przy którym jeszcze
spełnia on zadane wymagania. Funkcje te zebrane zostały w tabelce 3. Są one bardzo
przydatne i często stosowane w połączeniu z funkcjami do konstrukcji filtrów.
Załóżmy, że poszukuje się filtru środkowoprzepustowego z pasmem przepuszczania od 1000
do 2000 Hz, pasmem zaporowym zaczynającym się przy 500 i 2500 Hz, o częstotliwość
próbkowania 10 000 Hz, z maksymalnie 1 dB błędem w paśmie przepuszczania i z co
najmniej 60 dB tłumieniem w paśmie zaporowym. Do konstrukcji użyto funkcji

buttord

.

Tabela 3 Zestawienie funkcji wyznaczających minimalny rząd filtrów.

Typ filtru

Klasa estymacji funkcji

Butterworth

[n,Wn]=buttord(Wp,Ws,Rp,Rs)

Chebyshew typ
I

[n,Wn]=cheb1ord(Wp,Ws,Rp,Rs)

Chebyshew typ
II

[n,Wn]=cheb2ord(Wp,Ws,Rp,Rs)

Elliptic

[n,Wn]=ellipord(Wp,Ws,Rp,Rs)

[n,Wn]=buttord([1000 2000]/5000,[500 2500]/5000,1,60)

•

n =

•

12

•

Wn =

•

0.1951 0.4080

[b a]=butter(n,Wn);

Filtr eliptyczny, spełniający te same wymagania, jest dany przez funkcję

ellip

.

[n,Wn]=ellipord([1000 2000]/5000,[500 2500]/5000,1,60)

•

n =

•

5

•

Wn =

•

0.2000 0.4000

[b a]=ellip(n,1,60,Wn);

Przybornik dostarcza pięć typów klasycznych filtrów IIR, z których każdy jest w sensie
pewnego kryterium optymalny. W rozdziale tym pokazano podstawowe postacie
analogowych prototypów dla tych filtrów oraz podano ogólną charakterystykę każdego z
nich.

Bezpośrednie projektowanie filtrów IIR
W przyborniku użyto terminu metody bezpośrednie do opisu metody projektowania filtrów
IIR, która znajduje filtr bazując na stawianych mu wymaganiach w dziedzinie dyskretnej.
Metoda bezpośrednia, odmiennie niż metoda analogowego prototypowania, nie jest
ograniczona standardowymi konfiguracjami filtrów (dolno-, środkowo- czy górno-
przepustową). Służy ona raczej do projektowania dowolnych filtrów o zadanej odpowiedzi
częstotliwościowej.

Funkcja

yulewalk

zrealizowana w SPT umożliwia projektowanie rekursywnego cyfrowego

filtru IIR poprzez dopasowanie odpowiedzi częstotliwościowej filtru do zadanego kształtu.
Nazwa funkcji

yulewalk

pochodzi od metody znajdowania przez nią współczynników filtru.

Funkcja ta znajduje odwrotną transformatę FFT dla idealnego zadanego widma mocy i
rozwiązuje zmodyfikowane równania Yule-Walkera na podstawie znajomości funkcji
autokorelacji. Składnia funkcji ma postać:

[b,a]= yulewalk(n,f,m);

Oblicza ona rzędowe wektory b i a składające się z n+1 współczynników licznika i
mianownika filtru IIR n-tego rzędu, którego charakterystyki amplitudowo-częstotliwościowe
aproksymują dane zawarte w wektorach f i m.
Wektor f jest wektorem punktów częstotliwościowych zawartych w granicach od 0 do 1,
gdzie 1 reprezentuje częstotliwość Nyquista. Wektor m jest wektorem zawierającym zadaną
amplitudę odpowiedzi dla punktów o częstotliwościach zawartych w f. Wektory f i m mogą

opisywać każdy składający się z liniowych odcinków kształt odpowiedzi częstotliwościowej,
włączając w to także odpowiedzi wielopasmowe (rys 7). Odpowiednik filtrów FIR tej funkcji
to

fir2

, która także projektuje filtr o odpowiedzi amplitudowej bazującej na dowolnym

kształcie składającym się z liniowych odcinków.
Należy zauważyć, iż funkcja

yulewalk

nie zawiera informacji o fazie, oraz brak jest

instrukcji pozwalającej na optymalizację rezultatów filtracji.

Dla przykładu przedstawiono projekt filtru pasmowego wykonanego za pomocą funkcji

yulewalk

. Wynikowa i przybliżana odpowiedź częstotliwościowa zostały przedstawione na

wykresie.

m=[0 0 1 1 0 0 1 1 0 0];
f=[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1];

[b,a]=yulewalk(10,f,m);
[h,w]=freqz(b,a,128);

plot(f,m,'black-',w/pi,abs(h),'black:')
legend('Zadana funkcja dla filtru',...

'Uzyskana charakterystryka amplitudowa filtru')

0

0 . 2

0 . 4

0 . 6

0 . 8

1

0

0 . 2

0 . 4

0 . 6

0 . 8

1

1 . 2

1 . 4

Z a d a n a f u n k c j a d l a f i l t r u

U z y s k a n a c h a r a k t e r y s t r y k a a m p l i t u d o w a f i l t r u

Rys. 7 Porównanie kształtu charakterystyki zadanej i uzyskanej przy pomocy funkcji youlewalk

Projektowanie filtrów typu FIR
Cyfrowe filtry ze skończonym czasem trwania odpowiedzi impulsowej (filtry FIR) mają w
porównaniu z filtrami dającymi nieskończoną odpowiedź impulsową (IIR) zarówno zalety jak
i wady.
Podstawowe zalety filtrów FIR to:

•

liniowa faza,

•

duża stabilność,

•

startowe warunki początkowe filtru posiadają skończony czas trwania.

Główną wadą filtrów FIR jest to, że często wymagają dużo wyższego rzędu niż filtry IIR do
osiągnięcia podobnych własności. Opóźnienie wprowadzane przez te filtry jest często dużo
większe niż dla podobnych filtrów typu IIR.

Liniowość fazy
Funkcje projektujące filtry FIR, konstruują filtry z liniową fazą. Współczynniki takiego filtru
posiadają parzystą lub nieparzystą symetrię. Zależnie od rodzaju tej symetrii oraz od tego, czy
rząd n filtru jest parzysty czy nieparzysty, liniowa faza filtru (zawarta w wektorze b o
długości n+1) posiada pewne ograniczenia co do jego charakterystyki amplitudowo-
częstotliwościowej.
Opóźnienia fazowe oraz opóźnienia grupowe filtrów FIR z liniową fazą są sobie równe i stałe
w całym paśmie przepuszczania. Dla n-tego rzędu takiego filtru całkowite opóźnienie
grupowe jest 2/n, a filtrowany sygnał jest opóźniony o 2/n (moduł jego transformaty Fouriera
jest skalowany przez odpowiedź amplitudową filtru). Te właściwości zabezpieczają kształt
przebiegu w paśmie przepuszczania, to znaczy nie wprowadzają dodatkowych zniekształceń
fazowych.
Funkcje

fir1

,

fir2

,

firls

, oraz

remez

umożliwiają projektowanie dolno-, górno-,

środkowoprzepustowych oraz środkowozaporowych filtrów FIR. Wszystkie one umożliwiają
dobór standardowych filtrów z liniową fazą typu I oraz II.
Uwaga: Ponieważ charakterystyka amplitudowo-częstotliwościowa filtrów typu II jest równa
zero przy częstotliwości Nyquista, za pomocą funkcji

fir1

nie może projektować filtrów

typu II w wypadku filtru górnoprzepustowego i środkowozaporowego. Dla parzystej wartości
n w tym przypadku,

fir1

dodaje jeden do rzędu filtru i tym samym tworzy filtr typu I.

Funkcje

firls

i

remez

umożliwiają projektowanie filtrów FIR z liniową fazą typu III i IV

pod warunkiem, że zostanie w ich wywołaniach zadany parametr 'hilbert' lub 'differentiator'.

Metoda okien
Rozważmy idealny, dolnoprzepustowy filtr cyfrowy z częstotliwością odcięcia

ω

0

[rad/s].

Moduł tego filtru wynosi jeden dla wszystkich częstotliwości z amplitudą mniejszą niż

ω

0

oraz zero dla częstotliwości z amplitudą pomiędzy

ω

0

i

π

. Jego odpowiedź impulsowa h(n)

wynosi:

)

n

(

c

sin

d

e

2

1

d

e

)

(

H

2

1

)

n

(

h

0

0

n

j

n

j

⋅

π

ω

⋅

π

ω

=

ω

⋅

π

⋅

=

ω

ω

⋅

π

⋅

=

∫

∫

ω

ω

−

ω

π

π

−

ω

( 7 )

gdzie n zmienia się od -

∞

do +

∞

.

Filtr ten nie jest realizowalny, ponieważ jego odpowiedź impulsowa jest nieskończona. W
celu utworzenia odpowiedzi impulsowej o skończonym czasie trwania należy zastosować
okno czasowe, ograniczające czas trwania h(n).
Przykładowy, filtr dolnoprzepustowy o długości 51 i częstotliwości odcięcia

ω

0.4

π

rad/sek

utworzony został za pomocą poniższej instrukcji:

b=0.4*sinc(0.4*(-25:25));

W tym przypadku do wyznaczenia odpowiedzi filtru zastosowano zwykłe okno prostokątne.
Według twierdzenia Parsevala spośród wszystkich filtrów o długości 51 ten z oknem
prostokątnym najlepiej przybliża idealny filtr dolnoprzepustowy.
Na wykresie (rys. 8) przedstawiono odpowiedź częstotliwościową tego filtru.

clf;

[H,w]=freqz(b,1,512,2);
plot(w,abs(H),'black'),

xlabel('Częstotliwość normalizowana do częst.
Nyquista'),

ylabel('Odpowiedź amplitudowa')

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

1.2

Częstotliwość znormalizowana do częstotliwości Nyquista

|H(j

ω

)|

Rys. 8 Charakterystyka częstotliwościowa przykładowego filtru dolnoprzepustowego

Na charakterystyce filtru, szczególnie w pobliżu krawędzi pasma, można zobaczyć oscylacje.
To zjawisko, nazywane efektem Gibbsa, nie znika dla wzrastającej długość filtru, ale
odpowiednie okno redukuje jego amplitudę. Mnożenie przez okno w dziedzinie czasu jest
równoważne splotowi lub wygładzaniu w dziedzinie częstotliwości.
Przykładowo dla filtru zaprojektowanego w poprzednim punkcie zastosowano okno
Hamminga o długości 51.

b=b.*hamming(51)';
[H,w]=freqz(b,1,512,2);

clf;
plot(w,abs(H),'black')

xlabel('Częstotliwość normalizowana do częst.
Nyquista'),

ylabel('Odpowiedź amplitudowa')

Jak można zaobserwować (rys. 9), oscylacje w znacznym stopniu zostały zredukowane.
Zostało to jednak okupione zmniejszeniem się nachylenia charakterystyki filtru oraz
wzrostem błędu kwadratowego.

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

1.2

Częstotliwość znormalizowana do częstotliwości Nyquista

|H(j

ω

)|

Rys. 9 Charakterystyka częstotliwościowa filtru wygładzona widmem okna Hamminga

Funkcje

fir1

i

fir2

w procesie projektowania filtrów FIR bazują na stosowaniu okien

czasowych. Mając zadany rząd oraz opis wymaganego filtru, funkcje te obliczają odwrotną
transformację Fouriera odpowiedzi zadanego filtru i poddają ją działaniu wybranego okna

czasowego. Obie standardowo używają okna Hamminga, można jednak stosować w nich okna
czasowe innego typu.

Projektowanie filtrów FIR dla pojedynczego pasma; funkcja

fir1

Funkcja

fir1

stosuje klasyczną metodę związaną z zastosowaniem okien do projektowania

cyfrowych filtrów FIR z liniową fazą. Przypomina ona projektowanie filtrów IIR w tym
sensie, że sformułowana jest do tworzenia filtrów w standardowych: dolno-, środkowo-,
górnoprzepustowych oraz środkowozaporowych konfiguracjach.
Sekwencja instrukcji:

n=50;
Wn=0.4;

b=fir1(n,Wn);

tworzy wektor

b

, zawierający współczynniki filtru n-tego rzędu, zaprojektowanego z

zastosowaniem okna Hamminga. Jest to dolnoprzepustowy filtr FIR z liniową fazą, o
częstotliwości odcięcia

Wn

.

Wn

jest liczbą z zakresu 0

÷

1, przy czym 1 odpowiada

częstotliwości Nyquista czyli połowie częstotliwości próbkowania. Występuje tu różnica w
stosunku do innych metod, gdzie

Wn

odpowiada punktowi o 6 dB spadku amplitudy.

Dla uzyskania górnoprzepustowego filtru należy dołączyć do listy parametrów ciąg znaków

‘high’

. Podczas tworzenia filtrów środkowo-przepustowych lub środkowo-zaporowych,

należy zadać

Wn

jako dwuelementowy wektor składający się z częstotliwości krawędzi pasma

przepuszczania, dodając parametr

‘stop’

dla konfiguracji zaporowej.

Instrukcja

fir1(n,Wn,okno)

stosowana do projektowania filtru zawiera okno podane w

wektorze kolumnowym

okno

. Wektor ten musi mieć długość n+1. Jeśli parametr nie zostanie

zadany, funkcja

fir1

standardowo zastosuje okno Hamminga.

Projektowanie filtrów FIR pasmowych; funkcja

fir2

Funkcja

fir2

umożliwia projektowanie filtrów FIR korzystając z metody okien, pozwalając

uzyskać dowolnie ukształtowaną charakterystykę częstotliwościową, przybliżaną za pomocą
odcinków liniowych. Tym różni się funkcja

fir2

w stosunku do funkcji

fir1

, która

umożliwia projektowanie filtrów jedynie w standardowych konfiguracjach.

n=50;
f=[0 0.4 0.5 1];

m=[1 1 0 0];
b=fir2(n,f,m);

Funkcja

fir2

daje w wyniku wektor wierszowy

b

, składający się z

n

+1 współczynników

filtru FIR

n

-tego rzedu, którego charakterystyka amplitudowo-fazowa jest dana poprzez

wektor

f

i

m

, gdzie

f

jest wektorem punktów częstotliwościowych w zakresie od 0 do 1 (1

reprezentuje częstotliwość Nyquista), zaś

m

jest wektorem składającym się z zadanych

odpowiedzi amplitudowych w punktach podanych w wektorze

f

. Odpowiednikiem tej funkcji

dla filtrów IIR jest funkcja

yulewalk

, tworząca filtry o odpowiedzi amplitudowej, bazującej

na dowolnym kształcie aproksymowanym odcinkami.

Projektowanie pasmowych filtrów FIR z zadanym pasmem przepuszczania
Funkcje

firls

i

remez

wykorzystują bardziej ogólne metody umożliwiające zadawanie

idealnego filtru, który ma być aproksymowany przez filtr cyfrowy.

Firls

i

remez

tworzą

transformatory Hilberta, układy różniczkujące oraz inne filtry posiadające nieparzystą
symetrię współczynników (typy III i IV filtrów z liniową fazą). Pozwalają one także na
zadawanie obszarów przejść, w których błędy nie są minimalizowane, podczas gdy w
pozostałej części charakterystyki filtru minimalizacja taka jest dokonywana.

Funkcja

firls

jest rozszerzeniem funkcji

fir1

i

fir2

o minimalizację całkowego błędu

kwadratowgo pomiędzy zadaną a rzeczywistą odpowiedzią częstotliwościową.
Funkcja

remez

wykorzystuje algorytm Pearks-McClellana, który przy wykorzystaniu

algorytmu wymiany Remeza oraz teorii aproksymacji Chebyshewa umożliwia projektowanie
filtrów, których rzeczywista odpowiedź częstotliwościowa jest w sposób optymalny
dopasowana do zadanej. Filtry te są optymalne w tym sensie, że minimalizują maksymalny
błąd pomiędzy zadaną i rzeczywistą odpowiedzią częstotliwościową (nazywane często
filtrami minimaksowymi). Filtry projektowane tą drogą wykazują zachowanie stałych
oscylacji w ich odpowiedziach częstotliwościowych. Algorytm Parks-McClellana
projektujący filtry FIR jest jednym z najbardziej znanych i szeroko stosowanych w praktyce
projektowania filtrów FIR.
Składnia dla funkcji

firls

i

remez

jest taka sama, jedyne różnice występują w

zastosowanych algorytmach minimalizacji. W kolejnym przykładzie pokazano, jak można
projektować filtry przy pomocy funkcji

firls

i

remez

z uwzględnieniem tych różnic.

Standardowy sposób działania funkcji

firls

i

remez

w projektowaniu filtrów typu I i II z

liniową fazą zależy od tego, czy zadany rząd filtru jest parzysty czy nieparzysty. Przykładowy
filtr dolnoprzepustowy z amplitudą 1 dla zakresu częstotliwości od 0 do 0.4 Hz oraz
amplitudą 0 dla przedziału od 0.5 do 1 Hz, zadano za pomocą następujących przedziałów.

n=20;
f=[0 0.4 0.5 1];

m=[1 1 0 0];
b=remez(n,f,m);

W zakresie częstotliwości od 0.4 do 0.5 Hz, funkcja

remez

nie dokonuje minimalizacji

błędów, ponieważ jest to obszar nachylenia zbocza filtru. Minimalizacja błędu w tym paśmie
jest również możliwa, ale należy liczyć się z dłuższym procesem przejściowym.
Dla porównania metody najmniejszych kwadratów projektowania filtrów z minimaksowym
projektowaniem filtrów, użyto funkcji

firls

do zbudowania podobnego filtru, jak w

przykładzie poprzednim:

bb=firls(n,f,m);

a następnie porównano ich charakterystyki częstotliwościowe.

[H,w]=freqz(b);

[HH,w]=freqz(bb);
plot(w/pi,abs(H),'black',w/pi,abs(HH),'black:')

grid
xlabel('Częstotliwość normalizowana do częst.

Nyquista'),
ylabel('Odpowiedź amplitudowa')

legend('Funkcja remez','Funkcja firls')

Na rysunku 10 można łatwo zauważyć, że filtr zaprojektowany za pomocą funkcji

remez

wykazuje stałe oscylacje tak w paśmie przepuszczania, jak i zaporowym. Filtr projektowany
funkcją

firls

posiada lepszy przebieg charakterystyki dla większej części obu pasm, ale

przy krawędziach (f=0.4 i f=0.5) jego charakterystyka ma przebieg dużo gorszy od
poprzedniego. Świadczy to o tym, że maksymalny błąd filtru

remez

ponad pasmem

przepuszczania i zaporowym jest mniejszy i w rzeczywistości ma on dla tej konfiguracji
krawędzi pasma i długości filtru wartość najmniejszą.

Funkcja remez

Funkcja firls

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

1.2

Częstotliwość znormalizowana do częstotliwości Nyquista

O

d

po

w

ie

d

ź

a

m

pl

itu

do

w

a

Rys. 10 Porównanie wyników działania funkcji remez i firls

W funkcjach

firls

i

remez

poszczególne pasma filtru traktowane są jako odcinki liniowe

reprezentujące przedziały częstotliwości. Bazując na funkcji o dowolnym kształcie,
składającym się z liniowych odcinków, można przedstawić filtr w dowolnej konfiguracji.
Funkcje

firls

i

remez

umożliwiają projektowanie filtrów dolno-, górno-,

środkowoprzepustowych oraz środkowozaporowych.
Przykłady parametrów dla wybranych filtrów zestawiono poniżej.

•

filtr środkowoprzepustowy:

f=[0 0.3 0.4 0.7 0.8 1];
m=[0 0 1 1 0 0];

(Wektory

f

i

m

definiuje się jako pięć pasm:

•

dwa zaporowe od 0.0 do 0.3 oraz od 0.8 do 1.0

•

pasmo przepuszczania od 0.4 do 0.7

•

dwa zbocza, od 0.3 do 0.4 i od 0.7 do 0.8 )

•

filtr górno- i środkowozaporowy:

f=[0 0.7 0.8 1];
m=[0 0 1 1];

f=[0 0.3 0.4 0.5 0.8 1];

m=[1 1 0 0 1 1];

•

filtr pasmowo-przepustowy:

f=[0 0.1 0.15 0.25 0.3 0.4 0.45 0.55 0.6 0.7 0.75 0.85 0.9 1];

m=[1 1 0 0 1 1 0 0 1 1 0 0 1 1];

Jeszcze innym przykładem zastosowania opisywanych funkcji jest projektowanie filtru,
posiadającego zbocze w formie liniowego połączenia pasma przewodzenia z zaporowym.
Zastosowanie takiego zabiegu może być pomocne w sterowaniu przeciekami odpowiedzi
amplitudowej przy słabo nachylonych zboczach.

f=[0 0.4 0.42 0.48 0.5 1];

m=[1 1 0.8 0.2 0 0];

Wektor wagi
Obie funkcje

firls

i

remez

pozwalają z różną wagą minimalizować błędy w pewnych

pasmach częstotliwości. Przykładowo, jeśli jako czwarty parametr w funkcji

remez

zostanie

podany wektor wagi

w=[1 10]

, zostanie zaprojektowany dolnoprzepustowy filtr z 10-krotnie

mniejszymi oscylacjami w paśmie zaporowym niż w paśmie przepuszczania.

n=20;

f=[0 0.4 0.5 1];
m=[1 1 0 0];

w=[1 10];
b=remez(n,f,m,w);

Wektor wagi jest zawsze połową długości wektorów

f

i

m

(musi być dokładnie jedna waga na

pasmo).

Projektowanie filtrów cyfrowych z wykorzystaniem narzędzi graficznych

W przyborniku SPT dostępne jest narzędzie udostępniające graficzny interfejs użytkownika
dla czynności opisywanych w poprzednich podpunktach konspektu. Narzędzie uruchamiane
jest poprzez wpisanie polecenia:

sptool;

w oknie poleceń Matlaba. Widok okna po uruchomieniu pokazany jest na rysunku 11.

Rys. 11 Widok okna narzędzia SPTool

Idea użycia narzędzia w celu wykonania filtracji sygnału jest następująca:

1. Sygnał obrabiany należy zaimportować do interfejsu (File ->Import...)
2. Utworzyć za pomocą narzędzi graficznych nowy projekt filtru (Przycisk New w

kolumnie Filters).

3. Zaprojektować filtr zgodnie z wymaganiami.
4. Zapisać utworzony projekt filtru.
5. Wykonać filtrację poprzez zaznaczenie w kolumnie Signals zaimportowanego sygnału

oraz w kolumnie Filters zaprojektowanego filtru. Operacja filtracji wykonywana jest
po naciśnięciu przycisku Apply.

Zadania do wykonania

1. Skomentować przykłady zaprezentowane w części konspektu dotyczącego doboru

częstotliwości próbkowania sygnału. Wyjaśnić krótko mechanizmy powstawania
niejednoznaczności w prezentowanych przypadkach.

2. Wygenerować kilkusekundowy sygnał będący złożeniem trzech przebiegów

sinusoidalnych o częstotliwościach: 10, 80 i 120 Hz. Amplitudy sygnałów powinny

znajdować się w proporcjach 1:3:1 oraz fazy powinny być przesunięte o ok. 20

o

.

Należy zwrócić szczególną uwagę na dobór próbkowania sygnału.

3. Zaprojektować filtry pozwalające na odfiltrowania kolejnych składowych ze sygnału.

Filtry powinny być wykonane każdorazowo w wersji IIR oraz FIR.

4. Wykonać porównanie parametrów zaprojektowanych filtrów zestawiając ze sobą filtry

FIR i IIR odpowiednio w konfiguracjach dolnoprzepustowych, pasmowych oraz
górnoprzepustowych.

5. Dla jednej wybranej konfiguracji filtru wykonać porównanie filtracji za pomocą

polecenia filter oraz filtfilt. Porównanie przedstawić w postaci graficznej poprzez
nałożenie na siebie:
•

sygnału oryginalnego – nieprzefiltrowanego,

•

sygnału odfiltrowanego poleceniem filter,

•

sygnału odfiltrowanego poleceniem filtfilt.

6. Na wygenerowany sygnał testowy uzyskany w punkcie 2 nałożyć sygnał szumu

losowego na poziomie 10% wartości amplitudy sygnału oraz skomentować wyniki
filtracji zaprojektowanymi wcześniej filtrami.