%% B2M31CZS cvičení - Návrh FIR a IIR filtrů z toleračního pole
%% Úkoly:
%% Návrh IIR pásmové propusti z tolerančního pole
% Navrhněte pásmovou propust IIR splňující požadavky následujícího tolerančního pole (využijte funkce buttord, cheb1ord, cheb2ord, ellipord):
% Hranice propustného pásma: 400 Hz, 3300 Hz
% Hranice nepropustného pásma: 200 Hz, 3500 Hz
% Vzorkovací kmitočet: 16000 Hz
% Maximální zvlnění v propustném pásmu: 1 dB
% Minimální potlačení v nepropustném pásmu: 40 dB
% Srovnejte vždy dosažené optimální řády navržených filtrů.

fpass1 = 400; % Hz
fpass2 = 3300; % Hz

fstop1 = 200; % Hz
fstop2 = 3500; % Hz

fs = 16000; % 16kHz
Rp = 1; % dB
Rs = 40; % dB

% Návrh IIR pásmové propusti
% help buttord, cheb1ord, cheb2ord, ellipord


% Opakujte návrh pro vzorkovací kmitočet fs = 44100 Hz a zvažte možnosti nalezení stabilního návrhu.
%
% Ověřte návrh filtrací. Signály ke zpracování jsou v souboru speech_8_16_44.mat (binární formát v MATLABu, obsahuje 3 signály uložené v proměnných sig8, sig16 a sig44 s fs = 8, 16 resp. 44.1 kHz.

load speech_8_16_44.mat % Obsahuje proměnné sig8, sig16 a sig44

fs = 8000; % = 16000, 44100 Hz, podle signálu

% Filtrace řečového signálu



% Zobrazte frekvenční charakteristiky navržených filtrů v dB pro všechny uvedené vzorkovací kmitočty fs = 8, 16 a 44.1 kHz.
% Srovnejte řády navržených filtrů a ověřte jejich stabilitu.
% Pro signál sig16 vzorkovaný kmitočtem 16 kHz realizujte filtraci stabilním elliptickým IIR filtrem a pozorujte spektrogramy vstupního a výstupního signálu.
%
%% Návrh FIR filtrů z tolerančního pole
% Sledujte tvar a spektrum Kaiserova okna pro různé N a β (funkce kaiser).
% Zobrazte Kaiserova okna délky stejné délky M = 101 vzorků pro různé hodnoty β = 0, 2, 8, 15 a zobrazte také jejich výkonová spektra v dB v rozlišení NFFT = 1024.
% Analogicky zobrazte Kaiserova okna délky se stejným β = 8 pro různé délky M = 31, 51, 101, 201 vzorků a rovněž i jejich výkonová spektra v dB v rozlišení NFFT = 1024.

% Zobrazení Kaiserova okna
M = [31, 51, 101, 201]; % vzorků, také: M = 31, 51, 101, 201
beta = [0, 2, 8, 15]; % také: beta = 2, 8, 15
NFFT = 1024;

figure("Name","Vizualizace Kaiserova okna");
TL = tiledlayout(2,4); % Vizualizace 2*4 subploty
title(TL, "Vizualizace Kaiserova okna");

for i = 1:length(beta)
    % Výpočet výkonového spektra Kaiserova okna s M = M(3) = 101 a beta = beta(i) v dB s NFFT = 1024 sem:
    
    

    % Vykreslení
    nexttile
    title(sprintf("Kaiserovo okno: M=%d, β=%d", M(3), beta(i)))
    


end

for j = 1:length(beta)
    % Výpočet výkonového spektra Kaiserova okna s M(j) a beta = beta(3) = 8 v dB s NFFT = 1024 sem:
    
    
    
    % Vykreslení
    nexttile
    title(sprintf("Kaiserovo okno: M=%d, β=%d", M(j), beta(3)))
    


end

%
% Navrhněte pásmovou propust Kaiserovou metodou (funkce kaiserord nebo fir1) splňující požadavky z první úlohy. (POZOR! Při návrhu FIR filtrů se parametry tolerančního pole zadávají mírně odlišně od funkcí pro návrh IIR filtrů).
% Sledujte dosažené řády navržených filtrů pro různé vzorkovací kmitočty a srovnejte s řády výše navržených IIR filtrů pro stejnou úlohu.

fpass1 = 400; % Hz
fpass2 = 3300; % Hz

fstop1 = 200; % Hz
fstop2 = 3500; % Hz

fs = 16000; % 16kHz
Rp = 1; % dB
Rs = 40; % dB

% Návrh pásmové propusti Kaiserovou metodou
% help kaiserord



%
% Opakujte návrh pro FIR filtr s rovnoměrným zvlněním podle Parks-McClellanovy aproximace (funkce firpmord nebo firpm) a pozorujte zvlnění ve frekvenční chrakteristice a dosažené řády navržených filtrů.

% Návrh pásmové propusti podle Parks-McClellanovy aproximace
% help firpmord
% help firpm



% Zobrazte frekvenční charakteristiky navržených filtrů v dB pro všechny uvedené vzorkovací kmitočty fs = 8, 16 a 44.1 kHz.
% Srovnejte řády navržených filtrů (porovnejte je rovněž s optimálními řády pro IIR filtry).

% Zobrazení frekvenčních charakteristik v dB pro fs = 8, 16 a 44.1 kHz
fs = 8000; % Také 16000, 44100



%
%% Srovnání zkreslení výstupního signálu při filtraci IIR a FIR filtry
% Pro výše uvedený řečový signál sig16 vzorkovaný frekvencí 16 kHz sledujte časové průběhy vstupního signálu a výstupních signálů po filtraci FIR filtry (navrženými Kaiserovou metodou a pomocí Parks-McClellanovy approximace) či IIR filtrem (se silně nelineární fázovou charakteristikou - stabilní elliptický filtr).
% Zobrazte ve třech obrázcích vždy vstupní signál sig16 a výstup každého z výše uvedených 3 filtrů a pozorujte zpoždění signálu resp. změnu tvaru signálu.
% Zobrazte také frekvenční charakteristiky použitých FIR a IIR filtrů a srovnejte jejich fázové charakteristiky.

% Srovnání zkreslení FIR-Kaiser, FIR-PM, a IIR-Ellip filtrů



%
%% Návrh úzkopásmové pásmové zádrže
% Navrhněte pásmovou zádrž ve frekvenčním pásmu 48\mathrm{Hz}<f<52\mathrm{Hz}  pro potlačení síťového rušení v EKG signálu pro vzorkovací frekvenci f_s =200\mathrm{Hz}.
% EKG: ekg1.asc, f_s =200\mathrm{Hz}, použijte funkci load pro načtení signálu do MATLABu.
% EKG se síťovým rušením 50Hz: ekg50.asc, f_s =200\mathrm{Hz}.

load ekg1.asc
% load ekg50.asc

fs = 200; % Hz