%% B2B31CZS cvičení 12 - Zpracování signálu ve více frekvenčních pásmech
%% Úkoly:
%% I) Rozdělení signálů do jednotlivých pásem
% Rozdělte signál SA107S06.CS0 (fs = 16 kHz, bez hlavičky, 16-bit PCM, pro načtení do MATLABu použít funkci loadbin.m) do dvou stejně širokých frekvenčních pásem, tj. 0-\mathrm{fs}/4 a\mathrm{fs}/4-\mathrm{fs}/2.
% Filtry navrhněte jako základní FIR filtry metodou okna, tj. pomocí funkce fir1. Řád filtru volte M = 50.

sig = loadbin("SA107S06.CS0");
fs = 16000; % 16kHz

%% KROK 1: Vytvoření signálu v dolním kmitočtovém pásmu
% Navrhněte filtr pro filtraci do DKP 0-\mathrm{fs}/4 a filtrujte výše odkazovaný signál.
% Výsledek:
% Vykreslete frekvenční charakteristiku navženého filtru do DKP.
% Vykreslete časový průběh a spektrogram zpracovávaného signálu (pro výpočet spektrogramu volte délku krátkodobého segmentu \mathrm{tlen}=32\mathrm{ms}.
% Vykreslete do jednoho obrázku spektrogramy původního a filtrovaného signálu.

% Vytvoření signálu v DKP
M = 50;


% Časový průběh a spektrogram
tlen = 0.032;


%% KROK 2: Vytvoření signálu v horním kmitočtovém pásmu
% Navrhněte filtr pro filtraci do HKP\mathrm{fs}/4-\mathrm{fs}/2 a filtrujte výše odkazovaný signál.
% Filtr do HKP navrhněte transformací filtru do DKP.
% Výsledek:
% Vykreslete do jednoho obrázku frekvenční charakteristiky obou navžených filtrů (do DKP i HKP) a ověřte splnění podmínky perfektní rekonstrukce.
% Srovnejte spektrogramy signálů filtrovaných do DKP i HKP.

% Vytvoření signálu v HKP


% Vykreslení frekvenční charakteristiky filtru a spektrogramu filtrovaného
% signálu


%
%% II) Změna vzorkovacího kmitočtu
%% KROK 3: Proveďte decimaci 1:2 vypuštěním sudých vzorků signálů v obou frekvenčních pásmech.
% Pozorujte spektrogramy všech signálů před decimací a po decimaci v obou frekvenčních pásmech - vysvětlete!
% Poznámka: Nepoužívat funkci decimate, která nerealizuje jenom prosté vypuštění vzorku!
% Výsledek:
% Zobrazte časové průběhy a spektrogramy signálů v DKP a HKP před decimací a po decimaci.

% Decimace v DKP a HKP - nepoužívejte funkci decimate!


% Zobrazení čas. průběhů a spektrogramů před a po decimaci


%% KROK 4: Proveďte interpolaci 1:2 pro obnovení původního vzorkovacího kmitočtu v následujících dílčích krocích:
% Proložte decimované signály v jednotlivých pásmech nulami a sledujte jejich spektrogramy.
% Proveďte interpolaci filtrací pro obě pásma. Pro filtraci použijte stejné filtry jako pro rozdělení vstupního signálu do dvou frekvenčních pásem. Sledujte opět spektrogramy výsledných signálů.
% Sledujte výkony signálů před a po decimaci s následnou interpolací. Vysvětlete rozdíl a navrhněte správnou korekci úrovně signálu.
% Výsledek - Srovnejte spektrogramy v jednotlivých pásmech u signálů:
% Po decimaci
% Po 1. kroku interpolace (proložení nulami)
% Po 2. kroku interpolace (filtraci)

% Interpolace v DKP a HKP


% Srovnání spektrogramů


%
%% III) Složení signálu z obou zpracovávaných pásem
% Signály z jednotlivých pásem  sečtěte a vyčíslete míru zkreslení na bázi výkonu chybového signálu (tj. po rozdílu vstupního a výstupního signálu).
% POZOR! Nezapomeňte při výpočtu chybového signálu na zpoždění generované FIR filtrem.
% Výsledek:
% Časový průběh a spektrogram finálního složeného signálu.
% Určete {\mathrm{SNR}}_e =10\cdot \log \left(\frac{P_s }{P_e }\right) pro zpracovávané signály (tj. poměr výkonů původního signálu ku chybovému signálu po rekonstrukci v dB - nezapoměňte na správnou synchronizaci signálů pro výpočet chybového signálu).

% Složení signálu


% Časový průběh a spektrogram složeného signálu


% Určení SNR


% Zvažte chybu vzniklou decimací a interpolací, srovnejte s chybovým signálem signálu složeného z pásmových signálů bez decimace a interpolace (vyčíslete opět SNR_e pro tento případ).
%
%% IV) Banka filtrů s perfektní rekonstrukcí
% Pro filtraci do dolního kmitočtového pásma použijte filtr s koeficienty:
% b0 = [ -0.0059, 0.0129, 0.0013, -0.0274, 0.0086, 0.0510,  -0.0338, -0.1001, 0.1243, 0.4688, 0.4688, 0.1243, -0.1001, -0.0338,  0.0510, 0.0086, -0.0274, 0.0013, 0.0129, -0.0059 ];
% Pro filtraci do horního kmitočtového pásma použijte filtr vycházející z výše uvedeného prototypu dolní propusti.
% Pozorujte splnění podmínky perfektní rekonstrukce, tj. |H_0 |^2 +|H_1 |^2 =1, pro banku filtrů vycházející z tohoto prototypu a pro FIR filtry navržené v předchozím kroku metodou okna.

% Filtrace s pomocí filtru pro perfektní rekonstrukci
b0 = [ -0.0059, 0.0129, 0.0013, -0.0274, 0.0086, 0.0510,...
    -0.0338, -0.1001, 0.1243, 0.4688, 0.4688, 0.1243, -0.1001,...
    -0.0338,  0.0510, 0.0086, -0.0274, 0.0013, 0.0129, -0.0059 ];


%
%% PRÁCE NA DOMA:
% Proveďte rozklad včetně decimace a interpolace pro 4 lineární frekvenční pásma. Zaměřte se zejména na následující problémy:
% Jaký bude stupeň decimace?
% Jak se mění spektra signálů po decimaci v jednotlivých pásmech?
% Jak se změní postup interpolace? Sledujte spektra signálů po proložní nulami a určete vhodnou následnou filtraci.
% Vlastní příklady a motivace pro zpracování signálu ve více frekvenčních pásmech.