%% B2B31CZS cvičení 11 - Číslicová filtrace ve frekvenční oblasti, Metoda sčítání přesahů (OLA) s obecným oknem
%% Úkoly:
%% OLA s obecným oknem - základní implementace algoritmu
% POZNÁMKA: Tato část je opakování, následující úkoly byly součástí minulého cvičení, úkoly 1 - 2.
% Vytvořte posloupnost reprezentující impulzní odezvu filtru klouzavého průměru (MA) řádu M=100, tj. délky N_1 =101 vzorků.
% Vytvořte posloupnost délky N_2 =4000 vzorků reprezentující signál x\left\lbrack n\right\rbrack, která bude obsahovat součet:
% Sinusového signálu s nulovým fázovým posuvem a parametry A=1,f=12\mathrm{Hz},f_s =8000\mathrm{Hz}
% Bílého šumu s normálním rozložením s nulovou střední hodnotou a jednotkovým rozptylem.
% Realizujte filtraci výše uvedeným filtrem pro potlačení přítomného náhodného šumu na bázi lineární konvoluce (funkce conv) těchto dvou posloupností a zobrazte:
% Impulsní odezvu filtru klouzavých průměrů h\left\lbrack n\right\rbrack (stem)
% Vstupní signál x\left\lbrack n\right\rbrack (plot)
% Výstupní signál y\left\lbrack n\right\rbrack, tj. výsledek lineární konvoluce y\left\lbrack n\right\rbrack =x\left\lbrack n\right\rbrack *h\left\lbrack n\right\rbrack (plot)

M = 100; % Řád MA filtru
N1 = 101; % Délka imp. odezvy MA filtru
N2 = 4000; % Délka signálové směsi
A = 1; % Amplituda sinusového signálu
f = 12; % Hz
fs = 8000; % Hz

% Vytvoření signálů


% Filtrace lineární konvolucí


%
%% Realizace metody sčítání přesahů (OLA) s obecným oknem pro filtraci dlouhého signálu ve frekvenční oblasti
% Parametry metody a jednotlivé dílčí kroky postupu jsou následující:
% Délka okna: \mathrm{wlen}=512
% Krok segmentace (příp. překryv okna - overlap): \mathrm{wstep}=\mathrm{wlen}/2
% Typ váhovacího okna \mathrm{wtype}: Hannovo
% Výpočet počtu zpracovávaných oken: \mathrm{wnum}
% Generace nulového výstupní signálu (preallocation)
% Hlavní cyklus pro zpracování i-tého segmentu (výběr segmentu, modifikace spektra, přičtení k výstupní posloupnosti)
% Korekce pro nastavení jednotkového zesílení: \mathrm{wnorm}
% Vyzkoušejte jak bude vypadat výsledek při použití různých oken (Hammingovo, Blackmanovo, obdélníkové) a různých úrovní překryvu (25%, 50%, 75%).
% Zobrazte a srovnejte výstupní signály vypočtené na bázi lineární konvoluce a s využitím metody OLA s obecným oknem při použití Hammingova okna s 50% (případně 75%) překryvem.

wlen = 512;
wstep = wlen/2;
wtype = "hann";

% wnum = ;
% y = ; % Prealokace výstupního vektoru
% Hlavní cyklus:


% Korekce jednotkového zesílení:


% Srovnání výstupních signálů lineární konvoluce a OLA s obecným oknem


%
%% Implementace pásmové filtrace ve frekvenční oblasti
% Realizujte ve frekvenční oblasti pásmovou filtraci do telefonního pásma pro signál vzorkovaný kmitočtem f_s =16\mathrm{kHz} s využitím výše uvedené metody OLA s Hammingovým oknem a 75% překryvem. Délku okna volte 32 ms!
% Filtraci realizujte vytvořením masky nulující frekvenční komponenty komplexního spektra mimo propustné pásmo 300\mathrm{Hz}\le f\le 3400\mathrm{Hz}.
% Zobrazte vytvořenou masku filtru (ideální frekvenční charakteristiku).
% Zobrazte časový průběh a spektrogram vstupního signálu SA107S06.CS0 (f_s =16\mathrm{kHz}, bez hlavičky, 16-bit PCM, signál lze načíst pomocí funkce loadbin.m, délku okna pro krátkodobou analýzu při výpočtu spektrogramu volte \mathrm{tlen}=32\mathrm{ms}).
% Zobrazte časový průběh a spektrogram výstupního signálu po filtraci do telefonního pásma.
% Vstupní a výstupní signály srovnejte rovněž poslechem.

% Vytvoření a zobrazení masky filtru


% Časový průběh a spektrogram vstupního signálu
fs = 16000; % Hz
tlen = 0.032; % s
sig = loadbin("SA107S06.CS0"); % Signál


% Filtrace do telefonního pásma


% Zobrazení časového průběhu a spektrogramu výstupního signálu


% Opakujte při použití překryvu pouze 50% a informativním poslechem ověřte kvalitu výstupního signálu.
% Opakujte při použití pravoúhlého okna a poslechem ověřte vliv nespojitostí na hranicích segmentů.
% Opakujte případně s dalšími okny a různými úrovněmi překryvu a informativním poslechem ověřte kvalitu výstupního signálu.
% Opakujte pro řečové signály s různými vzorkovacími kmitočty uložené 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, pro načtení použijte příkaz load speech_8_16_44.mat).
% Zobrazte a srovnejte časové průběhy a spektrogramy širokopásmového vstupního signálu sig44 vzorkovaného kmitočtem \mathrm{fs}=44\ldotp 1\mathrm{kHz} a výstupního signálu omezeného do telefonního pásma.
%
%% Implementace spektrálního prahování a spektrálního odečítání pro potlačení bílého šumu ve směsi sinusovek
% Realizujte potlačení aditivního šumu v signálu metodou spektrálního prahování a OLA s obecným oknem pro modelovaný signál x\left\lbrack n\right\rbrack =s\left\lbrack n\right\rbrack +b\left\lbrack n\right\rbrack:
% Čistý signál s\left\lbrack n\right\rbrack je součet 2 sinusových signálů s nulovým fázovým posuvem s následujícími parametry:
% Amplituda prvního sinu: A_1 =0\ldotp 8
% Frekvence prvního sinu: f_1 =100\mathrm{Hz}
% Amplituda druhého sinu: A_2 =0\ldotp 5
% Frekvence druhého sinu: f_2 =427\mathrm{Hz}
% Vzorkovací kmitočet (v obou případech): f_s =8000\mathrm{Hz}
% Délka signálu: t_0 =0\ldotp 5s
% Aditivní šum b\left\lbrack n\right\rbrack je bílý šum s:
% Normálním rozložením
% Nulovou střední hodnotou: \mu =0
% Rozptylem: \sigma^2 =1\ldotp 5
% Délka signálu: t_0 =0\ldotp 5s

% Vytvoření vstupních signálů a jejich směsi
% Siny:
A1 = 0.8;
f1 = 100; % Hz
A2 = 0.5;
f2 = 427; % Hz
fs = 8000; % Hz
t0 = 0.5; % s


% Šum:
m = 0;
s = 1.5;
% t0 = 0.5;


% Úroveň prahu pro spektrální prahování zvolte vhodně nad spočítanou úrovní bílého šumu v amplitudovém spektru.
% Amplitudové spektrum bílého šumu pro spektrální odečítání odhadněte z výkonu bílého šumu určeného průměrováním v horní polovině frekvenčního pásma.
% Po odečtení odhadnutého amplitudového spektra aditivního šumu vynulujte záporné hodnoty a pro převod do časové oblasti použijte fázi vstupního (zašuměného) signálu.

% Spektrální prahování



% Zobrazte krátkodobé amplitudové spektrum vstupního zašuměného signálu v prvním zpracovávaném segmentu.
% Vykreslete do odhadu spektra hladinu odpovídající amplitudovému spektru bílého šumu a zvolenou úroveň prahu.
% Zobrazte časové průběhy vstupního signálu a signálu s potlačeným šumem.
% Vyčíslete SNR pro vstupní zašuměný signál a pro výstupní signál po spektrálním prahování (SNR z prostředního úseku po odstranění počátečního a koncového neúplného výstupu, např. pro vzorky 1000:3000).

% Zobrazení výsledků



%
%% PRÁCE NA DOMA:
%% Implementace spektrálního odečítání
% Realizujte potlačení aditivního šumu v signálu SA107S06_auto2.CS0 (věta,  f_s =16\mathrm{kHz}, bez hlavičky, 16-bit PCM) metodou spektrálního odečítání a OLA s obecným oknem.
% Odhad výkonu šumu určete průměrováním amplitudového spektra v počátečním úseku bez řečové aktivity.
% Po odečtení amplitudového spektra aditivního šumu negativní složky vynulujte a pro zpětný převod do časové oblasti použijte fázi vstupního (zašuměného) signálu.
% Zobrazte krátkodobé amplitudové spektrum vstupního zašuměného signálu v prvním zpracovávaném segmentu.
% Vykreslete do odhadu spektra hladinu odpovídající amplitudovému spektru šumu.
% Zobrazte časové průběhy a spektrogramy vstupního signálu a signálu s částečně potlačeným šumem.
% Signály si ilustrativně poslechněte.
% Další možné signály k filtraci:
% SA110992_auto1.CS0 - (slovo) f_s =16\mathrm{kHz}, bez hlavičky, 16-bit PCM
% ma034014_auto3.ils - (věta ve stojícím automobilu) f_s =16\mathrm{kHz}, bez hlavičky, 16-bit PCM
%
%% Realizace potlačení aditivního šumu ve vlastním nahraném řečovém signálu
% Nahrajte si vlastní promluvu s rušivým pozadím vzorkovanou frekvencí 16 kHz. Při nahrávání zachovejte minimálně 1 vteřinovou pauzu na začátku promluvy a amplitudové spektrum šumu pozadí odhadněte průměrováním krátkodobých spekter v počátečním úseku bez řečové aktivity. Realizujte spektrální odečítání dle postupu výše.
% Zobrazte časové průběhy a spektrogramy vstupního (zašuměného) i výstupního signálu s potlačeným rušivým pozadím.
% Vstupní i výstupní signály srovnejte rovněž ilustrativním poslechem.