INSTYTUT AUTOMATYKI I

ROBOTYKI

CYFROWE PRZETWARZANIE SYGNAŁÓW

Sprawozdanie z wykonania ćwiczeń

laboratoryjnych nr 9 i 10

Wykonawca:

Prowadzący: laboratorium :

- Sebastian Polkowski gr. I6A1S1

- dr inż. Leszek Grad

Ćwiczenie 9

Projektowanie filtrów cyfrowych typu FIR, cz.I – filtracja dźwięku.

Ćwiczenie obejmuje zadania:

1.

Zaprojektowanie filtru FIR dolnoprzepustowego oraz pasmowo przepustowego.

2.

Dokonanie filtracji sygnału dźwiękowego metodami:

-

z definicji splotu;

-

korzystając z DFT;

-

wykorzystując funkcję conv;

3.

Wyświetlenie charakterystyki amplitudowej filtru.

4.

Zbadać wpływ rzędu filtru na jego jakość.

Obowiązuje znajomość materiału z zakresu projektowania filtrów cyfrowych typu FIR.

Rozwiązanie

Do rozwiązania zadania utworzyłem skrypt w programie Matlab. W celu przedstawienia
rozwiązania zamieszczam poszczególne linie skryptu oraz opis ich działania.

Na początku skryptu wczytuję dźwięk, nad którym będę pracował. Następnie do zmiennej M
zapisuję długość sygnału. Na potrzeby dalszych obliczeń do zmiennej tp zapisuję okres
próbkowania, który jest równy odwrotności częstotliości próbkowania sygnału.

[u,fs]=wavread(

'dzwiek.wav'

);

M = length(u);
tp = 1/fs;

Następnie tworzę 3 zmienne, które będą parametrami filtru, którymi będzie można sterować. Są

to: N – rząd filtru, fc – częstotliwość graniczna, fc2 – druga częstotliwość graniczna (dla filtru

pasmowo przepustowego).

fc = 3000;
fc2 = 400;
N = 303;

Nastęnie obliczam parametry ν

c

potrzebne do znalezienia odpowiedzi impulsowych filtrów.

Parametry te zapisuję do zmiennych niC oraz niC2. Następnie do zmiennej pomocniczej N2

zapisuję połowę wartości N, ponieważ będę z niej potem często kożystał.

omegaC = 2*pi*fc;
omegaC2 = 2*pi*fc2;
niC = omegaC * tp;
niC2 = omegaC2 * tp;
N2 = (N-1) / 2;

Następnie wyznaczam wektory h i h2, które oznaczają kolejno: odpowiedź impulsową filtru

dolnoprzepustowego oraz pasmowo przepustowego. Kożystam z następującego wzoru:

h=zeros(N,1);
h2=zeros(N,1);

for

k = -N2:-1

h(N2+k+1) = sin(k*niC)/(pi*k);
h2(N2+k+1) = sin(k*niC2)/(pi*k);

end

h(N2+1) = niC / pi;
h2(N2+1) = niC2 / pi;

for

k = 1:N2

h(N2+k+1) = sin(k*niC)/(pi*k);
h2(N2+k+1) = sin(k*niC2)/(pi*k);

end

h2 = h - h2;

Posiadając wektory odpowiedzi impulsowych oraz sygnał wejściowy mogę obliczyć sygnał

wyjściowy. Na początek realizuję tą czynność kożystając z definicji splotu:

Do wektora y zapisuję odpowiedź filtru dolnoprzepustowego, natomiast do zmiennej z –

pasmowo przepustowego.

y = zeros(M+N-1,1);
z = zeros(M+N-1,1);

for

k = 1:M+N-1

for

l = 1:N

k2 = k-l+1;

if

(k2<=0) s = 0;

elseif

(k2>M) s = 0;

else

s = u(k2);

end

y(k) = y(k) + h(l)*s;
z(k) = z(k) + h2(l)*s;

end

end

Wyznaczone odpowiedzi wyświetlam na wykresach. Niebieskim kolorem zaznaczyłem sygnał

wejściowy, a kolorem czerowym – wyjściowy.

figure(1);
subplot(2,1,1);
plot(u);
title(

'Wynik filtracji z definicji splotu: Filtr

dolnoprzepustowy'

);

hold

on

;

plot(y,

'r'

);

subplot(2,1,2);
plot(u);
title(

'Wynik filtracji z definicji splotu: Filtr pasmowo

przepustowy'

);

hold

on

;

plot(z,

'r'

);

Kolejnym krokiem będzie wyznaczenie odpowiedzi filtrów kożystając z DFT. Na początku

uzupełniam zerami wektory h, h2 oraz u.

h(end:N+M-1) = 0;
h2(end:N+M-1) = 0;
u(end:N+M-1) = 0;

Następnie wyznaczam transformaty DFT powyższych wektorów.

hdft = fft(h);
h2dft = fft(h2);
udft = fft(u);

Kolejnym krokiem będzie wymnożenie transformat: hdft z udft oraz h2dft z u. W wyniku tych

operacji otrzymuję wektory, których transformaty odwrotne dają odpowiedzi filtrów.

y2dft = hdft .* udft;
y2 = ifft(y2dft);
z2dft = h2dft .* udft;
z2 = ifft(z2dft);

Otrzymane odpowiedzi przedstawiam na wykresach.

figure(2);
subplot(2,1,1);
plot(u);
title(

'Wynik filtracji z wykozystaniem DFT: Filtr

dolnoprzepustowy'

);

hold

on

;

plot(y2,

'r'

);

subplot(2,1,2);
plot(u);
title(

'Wynik filtracji z wykozystaniem DFT: Filtr pasmowo

przepustowy'

);

hold

on

;

plot(z2,

'r'

);

0

1000

2000

3000

4000

5000

6000

7000

8000

9000

10000

-0.04

-0.02

0

0.02

0.04

Wynik filtracji z wykozystaniem DFT: Filtr dolnoprzepustowy

0

1000

2000

3000

4000

5000

6000

7000

8000

9000

10000

-0.04

-0.02

0

0.02

0.04

Wynik filtracji z wykozystaniem DFT: Filtr pasmowo przepustowy

Ostatnią metodą wyliczenia splotu wektorów odpowiedzi impulsowej oraz sygnału wejściowego

jest skożystanie z funkcji wbudowanej w środowisko MatLab – conv.

y3 = conv(u,h);
z3 = conv(u,h2);

Jak w poprzednicg przypadkach otrzymane odpowiedzi przedstawiam na wykresach.

figure(3);
subplot(2,1,1);
plot(u);
title(

'Wynik filtracji z wykozystaniem funkcji conv: Filtr

dolnoprzepustowy'

);

hold

on

;

plot(y3,

'r'

);

subplot(2,1,2)
plot(u);
title(

'Wynik filtracji z wykozystaniem funkcji conv: Filtr

pasmowo przepustowy'

);

hold

on

;

plot(z3,

'r'

);

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x 10

4

-0.04

-0.02

0

0.02

0.04

Wynik filtracji z wykozystaniem funkcji conv: Filtr dolnoprzepustowy

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x 10

4

-0.04

-0.02

0

0.02

0.04

Wynik filtracji z wykozystaniem funkcji conv: Filtr pasmowo przepustowy

Następnie przedstawiam charakterystyki częstotliwościowe zaprojektowanych filtrów. Są to

transformaty odpowiedzi impulsowych.

figure(4);
subplot(2,1,1);
plot(abs(fft(h)));
title(

'Charakterystyka amplitudowa filtru

dolnoprzepustowego'

);

subplot(2,1,2);
plot(abs(fft(h2)));
title(

'Charakterystyka amplitudowa filtru pasmowo

przepustowego'

);

Ostatnim elementem ćwiczenia będzie porównanie odpowiedzi filtru dolnoprzepustowego dla

różnych wartości N (rząd filtru).

Ćwiczenie 10

Projektowanie filtrów cyfrowych typu FIR, cz.II – filtracja obrazu.

Ćwiczenie obejmuje zadania:

1.

Dokonanie filtracji obrazu metodami:

Z definicji splotu dwuwymiarowego;

Wykożystując funkcję conv2

filtrami o odpowiedziach impulsowych:





















=

1

1

1

1

1

1

1

1

1

9

1

h

filtr dolnoprzepustowy (uśrednianie)





















−

−

−

−

=

0

2

0

2

8

2

0

2

0

h

filtr górnoprzepustowy (detekcja krawędzi)

2.

Dokonanie

operacji

wyostrzenia

obrazu

poprzez

złożenie

wyniku

filtracji

górnoprzepustowej z obrazem pierwotnym ( suma).

3.

Wy

świetlenie charakterystyk amplitudowych powyższych filtrów.

Rozwiązanie

Na początku zadania wczytuję obraz oraz jego mapę kolorów do zmiennych u oraz map.

Nastęnie obliczam wymiary obrazu i zapisuję je do zmiennych sz oraz w.

[u,map] = imread(

'gray.bmp'

);

[sz w] = size(u);

Na początku wykonam filtrację dolnoprzepustową. W tym celu tworzę wektor h, w którym

przechowam odpowiedź impulsową filtra. Nastęnie do zmiennej y zapisuję przefiltrowany obraz,

który otrzymałem wykonując splot macieży u oraz h. Powyższą czynność zrealizowałem

kozystając z funkcji MatLaba conv2(). Jako jeden z parametrów podałem macież u zrzutowaną

na typ double, aby zachować zgodność typów.

h = 1/9 * ones(3);
y = conv2(double(u),h);

Następnie wyświetlam otrzymany obraz oraz jego pierwotną postać dla porównania.

figure(1);
colormap(map);
subplot(2,1,1);
title(

'Obraz pierwotny'

);

image(u);
subplot(2,1,2);
title(

'Odpowied

ź

filtru dolnoprzepustowego'

);

image(y);

Kolejnym krokiem będzie wykonanie podobnych czynności dla filtru górnoprzepustowego. Na

początku tworzę macież h2 (odpowiedź impulsową) na podstawie danych z treści zadania.

Następnie do zmiennej z zapisuję obraz przefitrowany jako wynik funkcji conv2().

h2 = [0 -2 0; -2 8 -2; 0 -2 0];
z = conv2(double(u),h2);

Nastęnie wyświetlam pierwotny obraz wraz z jego przefiltrowaną postacią, na której widzimy

wyostrzone krawędzie.

figure(2);
colormap(map);
subplot(2,1,1);
image(u);
title(

'Obraz pierwotny'

);

subplot(2,1,2);
image(z);
title(

'Odpowiedz filtru górnoprzepustowego'

);

Nastęnie wykonam prostą operację, której wynikiem będzie uzyskanie wyostrzonego obrazu. W

tym celu zsumuję macież u zawierającą dane o obrazie pierwotnym z macieżą z, która

przechowuje obraz z wyostrzonymi krawędziami. Na początku zmniejszam macież z obcinając

odpowiednią liczbę elementów, aby jej wymiary odpowiadały wymiarom macieży z obrazem

pierwotnym. Nastęnie do zmiennej x przypisuję macież będącą sumą dwóch wcześniej

wspomnianych. Przed tym, macież u rzutuję do typu double.

z = z(1:375,1:500);
x = double(u) + z;

Uzyskany wyostrzony obraz przedstawiam na wykresie.

figure(3);
colormap(map);
subplot(2,1,1);
image(u);
title(

'Obraz pierwotny'

);

subplot(2,1,2);
image(x);
title(

'Obraz wyostrzony'

);

Ostatnim elementem ćwiczenia jest wyświetlenie charakterystyk amplitudowych filtrów. W tym

celu obliczam transformaty dwówymiarowe DFT i wyświetlam za pomocąfunkcji surf, która służy

do rysowania trójwymiarowych wykresów w formie płaszczyzny.

figure(4);
surf(abs(fft2(h)));
title(

'Charakterystyka amplitudowa filtru dolnoprzepustowego'

);

figure(5);
surf(abs(fft2(h2)));
title(

'Charakterystyka amplitudowa filtru górnoprzepustowego'

);

Wnioski

Ćwiczenie 9 polegało na wykonaniu operacji filtracji za pomocą zaprojektowanych przez siebie

filtrów dolno- oraz pasmowo przepustowego na dźwięku. W zadaniu tym mogliśmy

manipulować niektórymi parametrami: rzędem filtru, oraz częstotliościami granicznymi.

Zwiększając rząd filtru można było zaobserwować zwiększenie przesunięcia wymuszenia

względem odpowiedzi. Bardziej widoczny wpływ na wyjście filtra miały częstotliwości graniczne.

Zmniejszając pasmo przenoszenia można było zaobserować pogorszenie sięjakości dźwięku.

Ćwiczenie 10 polegało na dokonaniu filtracji dolno- i górnoprzepustowej. W odróżnieniu do

poprzedniego ćwiczenia odpowiedzi impulsowe filtrów zostały zadane w treści ćwiczenia.

Pierwszy filtr dawał w odpowiedzi zamazany obraz. Efektem działania drugiego filtru był obraz o

wyostrzonych krawędziach. W dalszej części sumując obraz pierwotny z obrazem o

wyostrzonych krawędziach otrzymałem wyostrzony obraz.