Juliusz

% Autor opracowania: Juliusz Tarnowski
% Proszę o nierozpowszechnianie bez zgody autora.
% W przypadku pytań lub uwag proszę o wiadomość na messengerze.
% Opracowanie zawiera materiał niezbędny do zaliczenia pierwszego kolokwium z laboratoriów PPS.
% Jest to jednocześnie skrypt przystosowany do uruchomienia w Matlabie.
% W Octave przed uruchomieniem należy w linii poleceń wpisać: pkg load signal, aby załadować odpowiednią bibliotekę.
% Jeśli pojawi się błąd, należy ją doinstalować.

%% 1) Komendy ogólne
clear all;            % Służy do wyczyszczenia workspace'a, czyli wszystkich zdefiniowanych zmiennych
                      % Należy używać go na początku nowego skryptu, lub na początku "rozdziału" skryptu
close all;            % Służy do zamknięcia wszystkich okienek – inaczej trzeba zamykać je ręcznie.


%% 2) Definiowanie osi czasu.
% Parametry:
t_max = 5;            % Czas trwania sygnału. Wyrażony w sekundach.
fp = 2000;            % Częstotliwość próbkowania, ilość próbek na sekundę. Wyrażona w hercach.
N = 10001;            % Ilość próbek w całym sygnale. Należy pamiętać o próbce "numer 0".

% Wystarczy zdefiniować dwie spośród powyższych zmiennych. Trzecią można wyznaczyć:
t_max = (N-1)/fp;     % Czas trwania sygnału.
fp = (N-1)/t_max;     % Częstotliwość próbkowania.
N = (t_max * fp) + 1; % Liczba próbek.


%% 3) Metody generowania osi.
% Na PPSach definiowane są dwie podstawowe osie: czasu oraz częstotliwości.
% Pojawiają się też inne, np. oś przesunięcia (wyrażona w sekundach lub ilości próbek).
% Przy tworzeniu każdego sygnału, widma czy charakterystyki należy zastanowić się, czy została zdefiniowana odpowiednia oś.
% Należy wybrać taką metodę, jaka jest dla Ciebie najwygodniejsza, ale warto pamiętać o wszystkich.

t = 0 : 1/fp : (N-1)/fp;   % Podstawowa metoda generowania osi czasu.
%   ^   ^^^^   ^^^^^^^^^   % Komentarze czytaj od dołu do góry. ;)
%   |    |          \ ... do ostatniej chwili w sygnale.
%   |    |     ^^^^  ^^^
%   |    |      |     \ (Częstotliwość próbkowania)
%   |    |       \ (Ostatni indeks próbki - jedna została zużyta na "chwilę zero".)
%   |     \ ... co jeden okres próbkowania...
%    \ Tablica wartości od zera...

t = 0 : 1/fp : t_max;      % To samo, co powyżej, ale z wyliczonym maksymalnym czasem.

t = (0:N-1) ./ fp;         % Najbardziej zwięzła forma generowania osi czasu. Niewydajna obliczeniowo.
%   ^^^^^^^ ^^ ^^
%      |    |   \ ... przez częstotliwość próbkowania.
%      |     \ ... podzielony... (Kropka oznacza, że wykonujemy operację na próbkach wektora, a nie na całym wektorze.)
%       \ Wektor z liczbami naturalnymi od 0 do N-1 wartości...

t = linspace(0,t_max,N);   % Funkcja linespace() najczęściej jest wykorzystywana do generowania osi częstotliwości.
%   ^^^^^^^^ ^ ^^^^^ ^
%      |     |   |   \ Ilość próbek.
%      |     |    \ Wartość maksymalna.
%      |      \ Wartość początkowa. W tym przypadku to chwila "0", ale może to być np. chwila rozpoczęcia się fonemu, lub przesunięcie.
%       \ Nazwa funkcji. Generuje dokładnie N równomiernie rozłożonych wartości, między 0 a t_max.


%% 4) Najważniejsze rodzaje sygnałów
% Definiując parametry sygnału warto dodawać indeks (np. "_s"), aby później nie było kolizji oznaczeń.

%% 4.1) Szum biały, gaussowki, normalny.
% Podstawowy rodzaj szumu. W razie wątpliwości użyj właśnie jego.
W_1 = 3; % Wariancja/sigma. Pełni rolę pseudoamplitudy. Zgodnie z prawem trzech sigm 99.8% wartości będzie leżało nie dalej, niż w odległości 3*W_1 od wartości średniej.
sr_1 = 9; % Wartość średnia sygnału.

s_1 = W_1*randn(1,N) + sr_1;
%     ^^^ ^^^^^ ^^^    ^^^^
%      |    |    |       \ Wartość średnia.
%      |    |    \ Argumenty funkcji: 1 wiersz, N kolumn. Pomijając pierwszy argument wygenerujemy macierz rozmiaru N x N.
%      |     \ Funkcja generująca losowe wartości o rozkładzie normalnym. Posiada literkę "n".
%       \ Wariancja.

%% 4.2) Szum jednorodny.
% Dodatkowy rodzaj szumu, służący najczęściej do symulowania błędu kwantyzacji.
A_2 = 3; % Pseudoamplituda. Odległość między skrajnymi wartościami szumu.
sr_2 = 6; % Wartość średnia.

s_2 = A_2*(rand(1,N)-1/2) + sr_2; %
%     ^^^  ^^^^ ^^^ ^^^^    ^^^^
%      |     |   |   |        \ Wartość średnia.
%      |     |   |    \ Unormowanie wartości średniej szumu.
%      |     |   |      (Standardowo rand() ma wartość oczekiwaną równą 1/2, a wygodniej jest, przesuniemy ją do 0.)
%      |     |   \ Argumenty funkcji: 1 wiersz, N kolumn.
%      |      \ Funkcja generująca losowe wartości o rozkładzie jednorodnym [0,1]. Nie posiada literki "n".
%       \ Pseudoamplituda.

%% 4.3) Sinusoida.
% Parametry sinusoid:
A_3 = 5; % Amplituda, czyli wychylenie od zera.
f_3 = 24; % Częstotliwość sinusoidy. Musi być mniejsza od połowy częstotliwości próbkowania, inaczej pojawi się aliasing.
fi_3 = pi/4; % Przesunięcie fazowe. Wartość wyrażona w radianach, czyli np. [0,2*pi] lub [-pi,pi]

s_3 = A_3*cos(t*2*pi*f_3+fi_3);
%     ^^^ ^^^ ^ ^^^^ ^^^ ^^^^
%      |   |  |  |    |    \ Przesunięcie fazowe.
%      |   |  |  |     \ Częstotliwość.
%      |   |  |   \ Przemnożenie osi czasu przez 2*pi, żeby unormować częstotliwość (móc wyrażać ją w wygodnych dla człowieka hercach, a nie radianach).
%      |   |   \ Oś czasu.
%      |    \ Jeśli prowadzący nie zwraca na to uwagi, lepiej używać funkcji cos(), a nie sin().
%      |      Wynika to z własności transformat oraz ze wzoru Eulera: e^(i*x) = cos(x) + i*sin(x)
%       \ Amplituda.

% Sinusoida o zmiennej częstotliwości.
% Ten sygnał na PPS jest najczęściej wykorzystywany do testowania filtrów.
% Zachęcam posłuchać, jak brzmi ten sygnał na stronie: https://en.wikipedia.org/wiki/Chirp
% Parametry:
f0 = 1;                      % Częstotliwość dla czasu 0s.
t1 = t_max;                  % Koniec osi czasu.
f1 = 50;                     % Częstotliwość dla czasu t1.

s_4 = chirp(t, f0, t1, f1);  % Sinusoida o liniowo zmieniającej się częstotliwości.
%     ^^^^^ ^  ^^  ^^  ^^
%       |   |  |   |    \ Częstotliwość w zadanej chwili.
%       |   |  |    \ Wybrana w dowolny sposób chwila. Najwygodniej podać czas trwania sygnału (ostatni moment).
%       |   |   \ Częstotliwość w chwili 0.
%       |    \ Oś czasu. Jako jedyna jest wektorem, wszystkie pozostałe zmienne są liczbami.
%        \ Nazwa funkcji. Przyjmuje 4 argumenty.


%% 5) Wyświetlanie sygnałów
% Funkcje "porządkujące":
figure(1); % Wybór okienka.
subplot(2,2,1); % Wybór komórki.
%       ^ ^ ^
%       | | \ ... pierwsza komórka.
%       |  \ ... dwie kolumny...
%        \ Dwa wiersze...

% Funkcje "wyświetlające":
plot(t,s_1,'g--');               % Standardowy wykres, generowany linią ciągłą.
%    ^ ^^^ ^^^^^
%    |  |   \ Ozdobniki, czyli kolor i kształt wykresu.
%    |   \ Wartości sygnału. x oraz t muszą być takiego samego wymiaru.
%     \ Oś czasu.

% Dostępne kolory: 'r', 'g', 'b', 'c', 'm', 'y', 'k'.
% Dostępne ozgobniki: '-', '--', '.', 'o', '*' i inne.

subplot(222);
plot(t,s_3,'b', t,s_4,'r');    % Funkcja plot() przewiduje możliwość dodania większej ilości wykresów.

subplot(223);
stem(t,s_2,'ko');              % Wykres słupkowy, najczęściej wykorzystywany przy wizualizacji widma.

subplot(224);
nbits = 50;                    % Ilość słupków histogramu.
hist(s_1,nbits);               % Histogram, czyli wizualizacja zbioru wartości sygnału.
%    ^^^ ^^^^^
%     |    \ Drugi argument jest opcjonalny. Jeśli go nie podamy, Matlab automatycznie dobierze odpowiednią wartość.
%      \ Analizowany sygnał.

% Komendy specjalne:
% Poniższe komendy należy podać poniżej funkcji plot(). Wywołanie jej po subplot(), bez wywołania plot() nic nie zmieni.
hold on;  % Komenda pozwalająca na wyświetlanie wielu wykresów na jednej komórce.
hold off; % Antagonista powyższej komendy, używany znacznie rzadziej.
grid on;  % Aktywowanie siatki na wykresie. Ją również można dezaktywować.

% Funkcje generujące opis:
title('Tytuł wykresu, napisany pogrubioną czcionką, na górze.')
xlabel('opis osi x (zwykle to czas / częstotliwość)')
ylabel('opis osi y (zwykle to sygnał / widmo)')


%% 6) Wczytywanie nagrania
% Wczytanie nagrania wykorzystuje funkcję wavread(). Alternatywą dla niej jest audioread. Działają niemal identycznie.
[x, fpx] = wavread('mbi04wyczy.wav'); % Wczytanie sygnału i jego częstotliwości próbkowania.
%^  ^^^    ^^^^^^^ ^^^^^^^^^^^^^^^^
%|   |        |            \ Nazwa pliku (musi być w tym samym folderze co skrypt).
%|   |         \ Funkcja pozwalająca na wczytanie pliku. Można zastąpić ją funkcją audioread().
%|    \ Zmienna, do której przypiszemy częstotliwość próbkowania nagrania z pliku.
% \ Zmienna, do której przypiszemy sygnał z pliku.
Nx = length(x);                       % Ustalanie ilości próbek nagrania.
tx_max = (Nx-1)/fpx;                  % Czas trwania sygnału.
tx = 0 : 1/fpx : (Nx-1)/fpx;          % Generowanie osi czasu.


%% 7) Decymacja
% Decymacja to ograniczenie rozmiaru sygnału, poprzez usunięcie nadmiarowych próbek.
% Można to zrealizować na dwa sposoby: ręcznie, oraz z wykorzystaniem funkcji decimate().
% Decymacja ręczna nie jest kompletna, ponieważ pozwala na występowanie aliasingu.
% Funkcja decimate() posiada specjalny filtr antyaliasingowy.

dr = 6;                               % Rząd decymacji. Mogą pojawić się takie implementacje, gdzie musi być to potęga dwójki.
y_1 = x(1:dr:Nx);                     % Decymacja ręczna. W praktyce nie wykorzystywana.
%     ^ ^ ^^ ^^
%     | | |   \ (Kończymy na ostatniej.)
%     | |  \ ... co dr-tą próbkę.
%     |  \ (Zaczynamy od pierwszej.)
%      \ Oryginalny sygnał z którego pozostawiamy...

y_2 = decimate(x, dr);                % Decymacja sygnału z użyciem odpowiedniej funkcji.
%              ^  ^^
%              |   \ Rząd decymacji.
%               \ Sygnał do decymacji.

% Po decymacji należy wygenerować nową oś czasu:
Ny = length(y_2);                     % Ustalanie ilości próbek nagrania.
fpy = fpx / dr;                       % Nowa częstotliwość.
ty = 0 : 1/fpy : (Ny-1)/fpy;          % Generowanie osi czasu.


%% 8) Analiza częstotliwości
% Do analizy częstotliwości używamy funkcji fft().
% Wymaga ona podania jako argumentu liczby próbek widma, będących potęgą dwójki.
% Wynik będzie symetryczny, dlatego wyświetlać będziemy tylko połowę widma.
% Transformata Fouriera generuje wartości zespolone, więc przed narysowaniem wykresu należy przeliczyć je na liczby rzeczywiste.

Nfx = 2 ^ nextpow2(Nx);           % Liczba próbek widma.
%     ^   ^^^^^^^^ ^^
%     |      |      \ Ilość próbek sygnału.
%     |       \ Funkcja zwracająca wykładnik (ilość cyfr w systemie binarnym) pierwszej potęgi dwójki, która jest większa bądź równa zadanej liczbie.
%     |         (Na przykład dla wartości od 17 do 32 funkcja zwróci wartość 5, ponieważ 32=2^5.)
%      \ Dwójka, którą podniesiemy do potęgi wygenerowanej przez funkcję.

Nfx21 = Nfx/2 + 1;                % Zmienna pomocnicza, służąca wyświetlaniu tylko połowy widma (druga połowa zawiera takie same informacje).

fx = linspace(0, fpx/2, Nfx21);   % Oś częstotliwości.
%    ^^^^^^^^ ^  ^^^^^  ^^^^^
%       |     |    |      \ Ilość próbek, jakie chcemy wyświetlać.
%       |     |     \ Częstotliwość próbkowania (cz. Nyquista), czyli maksymalna częstotliwość, jaka interesuje nas w analizowanym widmie.
%       |      \ Zero, czyli minimalna interesująca nas częstotliwość.
%        \ Funkcja zwracająca tablicę wartości: Nf21 równie oddalonych próbek, których wartości zaczynają od 0, a kończą się na fp/2.

vx = fft(x, Nfx);          % Wyznaczenie widma zespolonego.
%    ^^^ ^  ^^^
%     |  |   \ ... oraz liczbę próbek widma, będącą potęgą dwójki.
%     |   \ ... sygnał...
%      \ Szybka transformata Fouriera, która jako argument przyjmuje...

avx = abs(vx);              % Widmo amplitudowe, czyli moduł transformaty Fouriera.
lavx = 20*log10(avx);       % Widmo amplitudowe w skali decybelowej / logarytmicznej.
%      ^^ ^^^^^ ^^^
%      |    |    \ Widmo, które przeliczamy na skalę decybelową.
%      |     \ ... na dekadę.
%       \ Dwadzieścia decybeli...
fvx = angle(vx);            % Widmo fazowe, czyli argument wartości transformaty Fouriera. Przyjmuje wartości ze skali [-pi, pi].
ufvx = unwrap(fvx);         % "Rozwinięte" widmo fazowe. Funkcja unwrap() nie pozwala na gwałtowne skoki wartości
                            % (Zamiast przeskoku na granicy wartości, będziemy kontynuować wykres.
                            % ... to trochę tak, jakbyśmy na zegarku zamiast wyświetlać po 23:59 godzinę 00:00, wyświetlilibyśmy 24:00.)

% Wizualizacja analizy widma:
figure(2);
subplot(311);
plot(tx,x);
title('Przebieg sygnału')
xlabel('Czas [s]');
ylabel('Wartość sygnału');

subplot(323);
plot(fx,avx(1:Nfx21));
title('Widmo amplitudowe')
xlabel('Częstotliwość [Hz]');

subplot(324);
plot(fx,fvx(1:Nfx21));
title('Widmo fazowe')
xlabel('Częstotliwość [Hz]');

subplot(325);
plot(fx,lavx(1:Nfx21));
title('Widmo amplitudowe w skali decybelowej')
xlabel('Częstotliwość [Hz]');

subplot(326);
plot(fx,ufvx(1:Nfx21));
title('Rozprostowane widmo fazowe')
xlabel('Częstotliwość [Hz]');


%% 9) Wycinanie fragmentu nagrania.
% Wycięcie sylaby "wy".
t0_wy = 0.11;                     % Chwila w której zaczyna się sylaba. Odczytujemy ją z wykresu.
t1_wy = 0.26;                     % Chwila w której kończy się sylaba.
pr0_wy = round(fpx*t0_wy);        % Numer próbki w której zaczyna się sylaba.
%   ^^   ^^^^^ ^^^ ^^^^^
%   |      |    |    \ ... przemnożona przez chwilę w której sylaba się zaczyna.
%   |      |     \ Częstotliwość próbkowania...
%   |       \ Zaokrąglenie. W przypadku nagrania nie ma znaczenia, czy w górę, czy w dół.
%   |         Dla sygnałów o niższej częstotliwości próbkowania warto użyć funkcji ceil() lub floor().
%    \ Polecam pamiętać o indeksach. Lub zmienić język programowania, np. na Pythona z NumPy. :P
pr1_wy = round(fpx*t1_wy);        % Numer próbki w której kończy się sylaba.
Nwy = pr0_wy - pr1_wy;            % Długość sylaby w próbkach.
wy = x(pr0_wy:pr1_wy);            % Wycięcie sylaby z nagrania.
twy = t0_wy : 1/fpx : t1_wy;      % Oś czasu dla sylaby.
%     ^^^^^   ^^^^^   ^^^^^
%       |       |       \ ... kończymy na chwili końcowej.
%       |        \ ... (nie zmieniając czêstotliwości próbkowania, bo sylaba jest fragmentem nagrania)...
%        \ Tym razem nie zaczynamy od zera, a od chwili początkowej sylaby...

t0_czy = 0.36;                    % Chwila w której zaczyna się sylaba. Odczytujemy ją z wykresu.
t1_czy = 0.5;                     % Chwila w której kończy się sylaba.
pr0_czy = round(fpx*t0_czy);      % Numer próbki w której zaczyna się sylaba.
pr1_czy = round(fpx*t1_czy);      % Numer próbki w której kończy się sylaba.
Nczy = pr0_czy - pr1_czy;         % Długość sylaby w próbkach.
czy = x(pr0_czy:pr1_czy);         % Wycięcie sylaby z nagrania.
tczy = t0_czy : 1/fpx : t1_czy;   % Oś czasu dla sylaby.

figure(3);
subplot(411);
plot(tx,x);
title('Słowo "wyczyść"')
xlabel('czas [s]');

subplot(423);
plot(twy,wy);
title('Sylaba "wy"')
xlabel('czas [s]');

subplot(424);
plot(tczy,czy);
title('Sylaba "czy"')
xlabel('czas [s]');


%% 10) Korelacja, autokorelacja
% Aby obliczyć korelację, należy podać ilość próbek, o jaką ma zostać przesunięty sygnał.
% Może być zadana na dwa sposoby:
kmax = 230;                      % Jako liczba naturalna.

tk_max = 0.15;                   % Lub jako czas, w tym przypadku 0.15 sekundy.
kmax = tk_max * fpx;             % Wyznaczenie ilości próbek.
%      ^^^^^^   ^^^
%        |       \ ... przemnożone przez ilość próbek na sekundę.
%         \ Maksymalne opóźnienie...

r = xcorr(wy,czy,kmax);          % Obliczanie autokorelacji.
%   ^^^^^ ^^ ^^^ ^^^^
%     |   |   |     \ ... ilość próbek, definiujące maksymalne opóźnienie.
%     |    \ ... dwukrotnie ten sam sygnał, oraz...
%      \ Funkcja do wyznaczania autokorelacji, która przyjmuje jako argumenty...

pr = -kmax : kmax;               % Oś przesunięcia (w próbkach).
tr = -kmax/fp : 1/fp : kmax/fp;  % Oś opóźnienia/przesunięcia (w sekundach).
%    ^^^^^^^^   ^^^^   ^^^^^^^
%       |        |        \ Maksymalne opóźnienie sygnału.
%       |         \ Czas próbkowania, wielkość co którą będą pojawiać się próbki czasu.
%        \ Ujemne opóźnienie sygnału.

figure(3);
subplot(413);                          % Wybór komórki dla trzeciego wykresu.
plot(tr, r);                           % Wykres autokorelacji.
xlabel('Przesunicie [s]');             % Opis poziomej osi.
title('Korelacja sylab "wy" i "czy"'); % Tytuł.

rwy = xcorr(wy,wy,kmax);               % Autokorelacja sylaby 'wy'.
rczy = xcorr(czy,czy,kmax);            % Autokorelacja sylaby 'czy'.

subplot(427);                          % Wybór komórki.
plot(tr, rwy);                         % Wykres autokorelacji.
xlabel('Przesunięcie [s]');            % Opis poziomej osi.
title('Autokorelacja sylaby "wy"');    % Tytuł.

subplot(428);                          % Wybór komórki.
plot(tr, rczy);                        % Wykres autokorelacji.
xlabel('Przesunięcie [s]');            % Opis poziomej osi.
title('Autokorelacja sylaby "czy"');   % Tytuł.


%% 11) Definiowanie filtra i filtracja
% Jest wiele rodzajów filtrów cyfrowych, tutaj omówione są filtry o skończonej odpowiedzi impulsowej.
% Filtr FIR (finite impulse response) wymaga:
% - Zdefiniowania długości odpowiedzi impulsowej.
% - Określenia zakresu filtracji w częstotliwości unormowanej.
%   (Częstotliwość unormowana jest zdefiniowana jako procent częstotliwości Nyquista.)
% - Określenia typu filtracji.

M = 401; % Długość odpowiedzi filtra. Powinna być to liczba nieparzysta.

fg = 3000;                   % Częstotliwość graniczna, wyrażona w hercach.
fgn = fg / (fpx/2);          % Częstotliwość graniczna unormowana.
%     ^^    ^^^^^
%     |       \ ... podzielona przez połowę częstotliwości próbkowania (częstotliwość Nyquista).
%      \ Częstotliwość graniczna...
hl = fir1(M-1, fgn, 'low');                     % Odpowiedź impulsowa filtra dolnoprzepustowego
%    ^^^^ ^^^  ^^^  ^^^^^^
%     |    |    |      \ Typ filtra. Są cztery typy, pozostałe zdefiniowano poniżej.
%     |    |     \ Częstotliwość graniczna
%     |     \ Długość odpowiedzi impulsowej. Dla zachowania spójności musi być pomniejszona o 1.
%      \ Funkcja generująca filtr.
hh = fir1(M-1, fgn, 'high');                    % Filtr górnoprzepustowy.

% Możliwe jest definiowanie filtrów, które modyfikują zakres widma.
fg_min = 5000;                                  % Dolny zakres częstotliwości granicznej.
fg_max = 7500;                                  % Górny zakres częstotliwości granicznej.
fgn_min = fg_min / (fpx/2);                     % Unormowana częstotliwość graniczna
fgn_max = fg_max / (fpx/2);                     % Jak wyżej.
hb = fir1(M-1, [fgn_min, fgn_max], 'bandpass'); % Filtr pasmowoprzepustowy.
hs = fir1(M-1, [fgn_min, fgn_max], 'stop');     % Filtr pasmowozaporowy.

ph = 0:(M-1);                                   % Oś próbek dla wykresu odpowiedzi impulsowej.
th = 0 : 1/fp : (M-1)/fp;                       % Os czasu.

% Aby zobaczyć działanie filtra odkomentuj wybraną linijkę.
h = hl;                         % 'low', filtr dolnoprzepustowy
% h = hh;                         % 'high', filtr górnoprzepustowy
% h = hb;                         % 'bandpass', filtr pasmowoprzepustowy
% h = hs;                         % 'stop', filtr pasmowozaporowy

y = filter(h,1,x);                % Filtracja. Sygnał przefiltrowany jest zapisany w zmiennej y.
%   ^^^^^^ ^ ^ ^
%     |    | | \ Filtrowany sygnał.
%     |    |  \ Jedynka. Wstawimy tu inną wartość, gdy będziemy korzystać z filtrów IIR (nieskończona odpowiedź impulsowa).
%     |     \ Odpowiedź impulsowa.
%      \ Funkcja realizująca filtrację.

vy = fft (y,Nfx);                 % Widmo zespolone sygnału po filtracji.
avy = abs(vy);                    % Widmo amplitudowe sygnału po filtracji.

% Filtry nie są sygnałami, tylko systemami, jednak możemy je analizować w dziedzinie częstotliwości.
% "Widmo" odpowiedzi impulsowej nazywamy transmitancją (transfer function, funkcja transmisji).
vh = fft (h,Nfx);                 % Transmitancja.
avh = abs(vh);                    % Moduł transmitancji.

figure(4);
subplot(321);
plot(tx,x);
xlabel('Czas [s]');
title('Sygnał przed filtracją');

subplot(322);
plot(fx,avx(1:Nfx21));
xlabel('Częstotliwość [Hz]');
title('Moduł widma sygnału');

subplot(323);
plot(th, h);
xlabel('Czas [s]');
title('Odpowiedź impulsowa');

subplot(324);
plot(fx,avh(1:Nfx21));
xlabel('Częstotliwość [Hz]');
title('Modul transmitancji');

subplot(325);
plot(tx, y);
xlabel('Czas [s]');
title('Sygnał po filtracji');

subplot(326);
plot(fx,avy(1:Nfx21));
xlabel('Częstotliwość [Hz]');
title('Moduł widma sygnału po filtracji');