%% B2M31CZS cvičení - Vlastnosti a návrh IIR filtrů
%% Úkoly:
%% IIR dolní propust
% Srovnejte vlastnosti IIR dolní propusti stejného řádu pro různé aproximace modulové charakteristiky:
% Butterworth - fc = 800 Hz, fs = 8 kHz, N = 6 (funkce butter)
% Chebyshev I (I. typu) - fc = 800 Hz, fs = 8 kHz, N = 6, Rp = 3 dB (funkce cheby1)
% Chebyshev II (II. typu - inverzní) - fc = 800 Hz, fs = 8 kHz, N = 6, Rs = 30 dB (funkce cheby2)
% Eliptický filtr (Cauerův) - fc = 800 Hz, fs = 8 kHz, N = 6, Rp = 3 dB, Rs = 30 dB (funkce ellip)
% Pozorujte dosaženou frekvenční charakteristiku navrženého filtru, tj. modulovou a fázovou charakteristiku (funkce freqz). Nastavte vhodně měřítko na ose 'y' (funkce ylim).
% Pozorujte rozložení nulových bodů a pólů přenosové funkce (funkce zplane) resp. impulzovou odezvu (funkce impz)- Zvažte, zda je navržený filtr stabilní.
% Vykreslete modulovou charakteristiku v dB pro všechny 4 navržené filtry do jednoho obrázku a srovnejte výsledek pro různé aproximace a stejný řád.

% Srovnání modulové charakteristiky v dB
figure("Name","Modulové charakteristiky různých typů filtrů")
TL = tiledlayout(2,2);
title(TL,"Modulové charakteristiky různých typů filtrů")

%

% Butterworth filtr
fc = 800; % Hz
fs = 8000; % 8kHz
N = 6; % Řád filtru
% help butter


%

% Chebyshev I. Typu
fc = 800; % Hz
fs = 8000; % 8kHz
N = 6; % Řád filtru
Rp = 3; % dB - Zvlnění v propustném pásmu (Ripple pass)
% help cheby1


%

% Chebyshev II. Typu (Inverzní)
fc = 800; % Hz
fs = 8000; % 8kHz
N = 6; % Řád filtru
Rs = 30; % dB - Zvlnění v nepropustném pásmu (Ripple stop)
% help cheby2


%

% Eliptický (Cauerův)
fc = 800; % Hz
fs = 8000; % 8kHz
N = 6; % Řád filtru
Rp = 3; % dB - Zvlnění v propustném pásmu (Ripple pass)
Rs = 30; % dB - Zvlnění v nepropustném pásmu (Ripple stop)
% help ellip


%
%% Filtrace bílého šumu
% Nagenerujte si Gaussovský bílý šum o délce N = 10000 vzorků.
% Vyfiltrujte tento šum každým ze 4 filtrů, které jste vytvořili v předchozích krocích.
% Vykreslete si obrázky, kde bude vedle sebe modulová frekvenční charakteristika použitého filtru v dB a spektrogram vyfiltrovaného bílého šumu (délku okna pro spektrogram volte 256 vzorků).

% Filtrace bílého šumu
N = 10000;
wlen = 256;


figure("Name","Filtrace bílého šumu různými filtry")
TL = tiledlayout(2,4)
title(TL,"Filtrace bílého šumu různými filtry")


%
% V případě volného času opakujte návrh pro změněné parametry a srovnejte výsledky.
% Zvyšte požadavek na potlačení v nepropustném pásmu na Rs=60
% Při požadovaném potlačení v nepropustném pásmu Rs=60 sledujte výsledek pro řády N=6, 10, 20
% Pozorujte také výsledek filtrace (časový průběh i spektrogram) pro filtr, který je nestabilní.
% Srovnejte analogicky pro návrh horní propusti se stejným mezním kmitočtem.
%
%% Návrh IIR pásmové propusti
% Navrhněte experimentálně pásmovou propust ve frekvenčním pásmu 300 < f < 3400, pro filtrování akustického signálu do telefonního pásma pro vzorkovací frekvence 8kHz, 16kHz a 44.1kHz.
% Potlačení v nepropustném pásmu by mělo být 40-60 dB, možné zvlnění v propustném pásmu 0.5-1 dB.
% Pro všechny filtry pozorujte frekvenční charakteristiku, impulsní odezvu resp. rozložení nulých bodů a pólů. V prvním kroku volte řád navrhovaných filtrů N = 6 a případně zvyšujte pro splnění výše zadaných podmínek.
% Pro ověření stability určete vždy maximální hodnotu modulu pólů dané přenosové funkce.

% IIR Pásmová propust
Rpmax = 0.5; % 0.5-1dB Zvlnění v propustném pásmu
N = 6; % Řád filtru
fmin = 300; % Hz
fmax = 3400; % Hz


%
%% Filtrace řečového signálu
% Pro stabilní filtry realizujte filtraci pro řečové signály 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. Ověřte dosažení požadovaného potlačení ve spektrogramu signálu.
% Pro každý filtrovaný signál s různým vzorkovacím kmitočtem vykreslete :
% Časové průběhy původního i filtrovaného signálu
% Spektrogramy původního i filtrovaného signálu, délku okna při výpočtu spektrogramu volte 32 ms
% Zobrazte spektrogramy všech 3 filtrovaných signálů (s délkou okna při výpočtu spektrogramu 32 ms) ve frekvenčním pásmu 0-4000 Hz (tj. upravte měřitko osy 'y' pomocí funkce ylim)

% Načtení signálu
load speech_8_16_44.mat % Načte proměnné sig8, sig16, sig44

% Filtrace
% Délka okna pro spektrogram = 32ms
% fs = ;
wlen = fs*0.032;
