%% B2B31CZS - Cvičení 4 - Odhady spektrálních charakteristik náhodných a nestacionárních signálů
%% Úkoly:
%% Modelujte signály pro účely následné spektrální analýzy
% s1 - 2 sinusovky různých kmitočtů

fs = 200; % Hz
tmax = 20 % s

f1 = 14 % Hz
A1 = 0.5;

f2 = 27.5; % Hz
A2 = 0.4;

% tt = ;
% s1 = ;

% b1 - bílý šum s Gaussovským rozložením
% b2 - bílý šum s rovnoměrným rozložením

% Pro fs = 200Hz a tmax = 20s
% Nb1 = ; 
% Nb2 = ;

% b1 = randn();
% b2 = rand();

% x1 - směs sinusovky s1 a bílého šumu b1 s výsledným SNR= -3 dB

% Koeficient 'k' pro směs s daným SNR zkopírujte z předešlých cvičení
SNR = -3; % dB - Požadované SNR

% x1 = ;

% Vytvořte požadovanou směs x1 a zobrazte první 2 segmenty o délce 512 vzorků.

% Zde vykreslete ze směsi x1 první dva segmenty o délce 512 vzorků
wlen = 512;

figure("Name","První dva segmenty signálu x1")


%
%% Vyhlazené odhady spektrální výkonové hustoty (PSD)
% Určete pomocí Welchovy metody (funkce pwelch) odhad výkonové spektrální hustoty (PSD) pro signál x1.

help pwelch

% Konkrétně sledujte:
% Výkonové spektrum v dB dvou různých segmentů o délce 512 vzorků (vypočtené přímo pomocí funkce fft, váhujte Hammingovým váhovacím okénkem)

% Vyberte ze signálu x1 dva různé segmenty - nemusí být nutně sousedící
% seg1 = x1();
% seg2 = x1();

% w = hamming();

% Výkonové spektrum segmentů
% Gseg1 = ;
% Gseg2 = ;

% Vyhlazený odhad PSD počítaný z celého signálu (funkce pwelch), váhujte opět Hammingovým váhovacím okénkem

% Vyhlazený odhad PSD
% Px1 = pwelch();

% Nevyhlazený odhad PSD počítaný z jednoho segmentu signálu (funkce pwelch), váhujte opět Hammingovým váhovacím okénkem

% Nevyhlazený odhad PSD
% Psegx1 = pwelch();

% Opakujte předchozí kroky a sledujte vliv váhování dalšími okénky (zejména pravoúhlým, Blackmanovým)

% w = rectwin();
% w = blackman();

%
% Vykreslete do jednoho obrázku výše uvedené 3 odhady pro modelovaný signál x1 (při váhování Hammingovým oknem)

figure("Name","Srovnání spektrálních odhadů signálu x1")



%
% Vykreslete do jednoho obrázku výše uvedené 3 odhady pro reálný EEG signál e1, viz e1.mat (fs = 200 Hz, binární MATLABovský mat-file, k načtení použijte funkci load)

fs = 200; % Hz
load("e1.mat") % Reálný EEG signál


figure("Name","Srovnání spektrálních odhadů EEG signálu e1")



%
%% Spektrogram stacionárního a nestacionárního signálu
%  Vykreslete spektrogram modelovaného signálu x1 (funkce spectrogram) a sledujte vliv jednotlivých volitelných parametrů při výpočtu spektrogramu:

help spectrogram 

% Sledujte zejména měnící se časově-frekvenční rozlišení pro různé délky krátkodobého okna (volte délku okna 512, 256, 128, 64 resp. 512, 1024, 2048, 4096).
% Nastavte korektní zobrazení frekvence na vertikální osu y (parametr Fs)
% Nastavte počet bodů DFT pro kratší segmenty a sledujte rozdíly pro výpočet s doplněním resp. bez doplnění nulami (parametr WINDOW pokud zadáte celé číslo)
% Sledujte vliv použitého váhovacího okna: rectwin (obdélníkové), bartlett (trojúhelníkové), hann (Hannovo), blackman (Blackmanovo - ostřejší)

% Délka segmentu/okna
% 512, 256, 128, 64 - srovnejte
% 512, 1024, 2048, 4096 - srovnejte
wlen = 512;

% Váhovací okno
wvtool(hamming(wlen)); % Nástroj pro zobrazení tvaru okna a jeho spektra
% w = hamming(); % blackman, hann, rectwin, bartlett

NFFT = 512;

fs = 200; % Hz

% Spektrogram pro x1



%
% Opakujte další signály s1, b1, b2 resp. e1 s optimálním nastavením pro každý signál

% Spektrogramy pro s1, b1, b2, e1



%
% Analyzujte spektrální charakteristiky nestacionárního řečového signálu r1 uloženého v souboru sm2.bin, který obsahuje promluvu "1 0 6 4 7" (muž) vzorkovanou kmitočtem fs = 16 kHz. (Binární signál lze načíst pomocí fce loadbin.m)
% Zobrazte spektrogramy signálu r1 pro segmenty o následujících délkách: 32 ms, 4 ms, 256 ms (použijte implicitní 50% překryv segmentů a Hammingovo okénko pro váhování) a sledujte měnící se časově-frekvenční rozlišení ve spektrogramu.

r1 = loadbin("sm2.bin"); % Mužský hlas, promluva "1 0 6 4 7"
fs = 16000; % Hz

tseg1 = 0.032; % 32 ms
tseg2 = 0.004; % 4 ms
tseg3 = 0.256; % 256 ms

% Nseg1 = ;
% Nseg2 = ;
% Nseg3 = ;


figure("Name","Spektrogramy signálu r1 pro různě dlouhá okna")

% Srovnejte spektrogram signálu r1 a dlouhodobý vyhlazený odhad PSD (funkce pwelch), obojí pro optimální délku krátkodobého segmentu 32 ms. Vysvětlete proč je nevhodné počítat PSD pro uvedený signál.

tseg1 = 0.032; % 32 ms
% Pr1 = pwelch();

figure("Name","Porovnání spektrogramu a vyhlazeného odhadu PSD pro segment 32ms")

%
% V případě dostatku času sledujte spektrogramy dalších promluv:
% sf2.bin - (promluva "1 0 6 4 7" - žena)
% sm1.bin - (promluva "Mobilní hlasová schránka" - muž - déšť)
% sf1.bin - (promluva "Mobilní hlasová shránka" - žena)
%
%% Intervaly spolehlivosti odhadu PSD
% Zobrazte vyhlazený odhad PSD pro signál x1, přičemž počítejte:
% Pro dříve vytvořený signál se SNR = -3 dB
% Pro modifikovanou směs se SNR = -10 dB
% Pro modifikovanou směs se SNR = 40 dB

% Směs se SNR -3 dB použijte z předchozích částí

% Vytvořte další dvě směsi s danými SNR níže
SNR2 = -10; % dB
SNR3 = 40; % dB



% Zobrazte PSD a intervaly spolehlivosti při hladině významnosti 98% pro výše uvedené 3 varianty signálu x1, tj. směsi sinusovky a šumu při SNR=-3, -10  resp. 40 dB.

figure("Name","PSD a intervaly spolehlivosti pro směs x1 se SNR -3, -10 a 40dB")



% Sledujte míru spolehlivost odhadu PSD pro uvedené 3 signály.
%
%% Ve volném čase
% Zobrazte případně 3-rozměrný prostorový spektrogram, viz help funkcí surf a spectrogram pro všechny výše analyzované signály.

help surf
help spectrogram