%% B2B31CZS - Cvičení 3
%% Úkoly:
%% DFT - periodogram, prosakování
% Generujte sinusový signál:

f = 1000; % Hz
fs = 8000; % Hz, Vzorkovací frekvence
t_max = 0.5; % s, Délka trvání

% Vyberte ze signálu segment o délce N=512 vzorků a sledujte jeho periodogram (funkce fft). Popište správně frekvenční osu a srovnejte různé způsoby zobrazení (plot s různými parametry, stem, apod.).
% Sledujte periodogram sinusového signálu pro segmenty různých délkách: N=512, 510, 508 vzorků resp. N=512, 515, 520 vzorků. Vysvětlete efekt délky okna (segmentu) na výsledcích, které uvidíte.
% Sledujte periodogram sinusovky o kmitočtu f=1231 Hz opět pro různé krátkodobého segmentu N=512, 510, 508 vzorků. Vysvětlete pozorované výsledky.

N = 512; % 510, 508, 515, 520 Měňte délku okna
% f = 1231; % Změňte frekvenci sinusovky

% Zopakujte předcházející analýzy na stejném signálu, který bude ale převáhován Hammingovým oknem (funkce hamming).

w = hamming(N);

% (Pokud budete rychlí, opakujte předchozí analýzu s jinými okny: Hannovo, Blackmanovo a trojúhelníkové váhovací okno - funkce hanning, blackman, triang.)
%% 1. Vliv váhování v detekci harmonických složek v signálu softwarově definovaného rádia
% Načtěte signál x_spl.mat obsahující směs harmonických složek v záznamu délky N = 250 vzorků. Signál je vzorkovaný kmitočtem fs = 1 kHz a je uložen v MATLABovském formátu. Pro načtení do MATLABu použijte funkci load.

load(x_spl.mat) % Nejdříve stáhněte a přidejte do aktuální složky/Matlabovské cesty
fs = 1000; % Hz
N = 250;

% Pro daný signál sledujte:
% Amplitudové spektrum signálu dané délky N = 250 bez váhování v lineárním i logaritmickém měřítku
% Amplitudové spektrum signálu váhovaného Hammingovým oknem v lineárním i logaritmickém měřítku
% Amplitudové spektrum signálu váhovaného Blackmanovým oknem v lineárním i logaritmickém měřítku
%% 2. Odhad krátkodobého výkonového spektra řečového signálu
% Načtěte signál vm7.bin, který je vzorkovaný kmitočtem fs = 16 kHz a je uložen v binární formě bez hlavičky. Pro načtení do MATLABu použijte funkci loadbin.m.

sig = loadbin("vm7.bin")
fs = 16000; % Hz

% Vykreslete výkonové spektrum v dB pro první dostupný segment délky N=512 vzorků se správným škálováním frekvenční osy.
% Srovnejte výsledek pro neváhovaný segment a segment váhovaný Hammingovým oknem. Sledujte vliv prosakování.
% (V případě dostatku času opakujte ilustrativně i pro další signály vf0.bin, vm0.bin, vf1.bin, vm1.bin, vf2.bin, vm2.bin, vf3.bin, vm3.bin, vf4.bin, vm4.bin, vf5.bin, vm5.bin, vf6.bin, vm6.bin, vf7.bin, vm7.bin, vf8.bin, vm8.bin, vf9.bin, vm9.bin.)
%
%% Doplnění nulami při výpočtu DFT (interpolace)
% Vygenerujte velmi krátký úsek sinusového signálu s parametry:

f = 1231; % Hz
fs = 8000; % Hz, Vzorkovací frekvence
N = 512; % Počet vzorků

% Počítejte spektrum pro délku diskrétní Fourierovy transformace NFFT = 512, 1024, 2048, 4096, 8192 a dosažené výsledky srovnejte.

NFFT = 512; % 1024, 2048, 4096, 8192

%% 1. Zpřesnění odhadu harmonických komponent v záznamu
% Načtěte signály sig1.mat a sig2.mat obsahující krátké záznamy směsi sinusových signálů za přítomnosti Gaussovského bílého šumu. Záznamy jsou vzorkovaný kmitočtem fs = 200 Hz a jsou uloženy v MATLABovském formátu. Pro načtení do MATLABu použijte funkci load.

load sig1.mat
load sig2.mat
fs = 200; % Hz

% Pro oba zvolené signály určete:
% Počet sinusových komponent a jejich kmitočty
% Spektrum spočítané ze segmentu dostupného záznamu o délce N=40 vzorků resp. ze segmentu doplněného nulami na délku N=512 vzorků
%% 2. Vliv interpolace v detekci harmonických složek v signálu softwarově definovaného rádia
% Pro dříve realizovaný výsledek spektrální analýzy signálu x_spl.mat obsahující směs harmonických složek v záznamu délky N = 250 vzorků, viz výše, sledujte vliv interpolace při výpočtu amplitudového spektra s různým typem váhováním a pro různé délky DFT, tj. zobrazte:
% Interpolované amplitudové spektrum signálu bez váhování v logaritmickém měřítku, které je počítané pomocí DFT délky NFFT = 2048 (doplněno nulami)
% Interpolované amplitudové spektrum signálu váhovaného Hammingovým oknem a následně také Blackmanovým oknem v lineárním i logaritmickém měřítku, které je počítané pomocí DFT délky NFFT = 2048 (doplněno nulami), určete také pro délku NFFT = 512 a NFFT = 4096 (doplněno nulami)
%
%% Detekce periodicity v DFT
% Pozorujte periodogramy následujících signálů, které obsahují jednu, dvě, tři resp. čtyři periody znělé hlásky vm0-1-per.bin,  vm0-2-per.bin, vm0-3-per.bin, vm0-4-per.bin. Pro načtení do MATLABu použijte opět funkci loadbin.m. Vzorkovací frekvence tohoto signálu je 16 kHz. Periodogramy počítejte vždy pro délky dané délkou záznamu. Pozorované výsledky vysvětlete.

loadbin("vm0-1-per.bin")
fs = 16000; % Hz

% Sledujte, jak se mění spektrum signálu vm0-4-per.bin, počítáme-li DFT následujících délek (tj. se zkrácením posloupnosti nebo při doplnění nulami):
% a) 479, 481
% b) 470, 490
% c) 256, 512
% d) 1024
% e) 2048
%% 1. Sledujte periodicitu v reálném signálu (znělá hláska) vm0-real-per.bin
% Pokuste se určit základní periodu resp. základní frekvenci daného kvaziperiodického signálu. Počítejte periodogram:
% a) pro NFFT = 512
% b) pro NFFT = 4096
%% 2. Sledujte spektrum neperiodického signálu (neznělá hláska) real-non-per.bin a srovnejte s předchozí variantou kvaziperiodického signálu
% Počítejte opět periodogram:
% a) pro NFFT = 512
% b) pro NFFT = 4096
%% V případě dostatku času...
% Periodicita ve spektru řečového signálu různé délky: Pro signál vm7.bin z prvního úkolu sledujte vliv délky okna na detekci periodicity ve spektru signálu. Pozorujte výkonové spektrum v dB pro první segmenty délky N=2000, 1024, 512, 256, 128, 64 vzorků, které vždy váhujte Hammingovým oknem stejné délky. Sledujte vždy analyzovaný váhovaný signál a vypočítané výkonové spektrum. Všimněte si souvislosti počtu pozorovaných period v dostupném segmentu signálu a zřetelnosti harmonických komponent v napočítaném spektru.
% Určete výšku hlasu na základě DFT i pro jiné záznamy odkazované v prvním bodě pro dnešní resp. minulé cvičení a srovnejte s výsledky dosaženými pomocí odhadu na bázi autokorelace. Pracujte vždy s vybraným krátkodobým segmentem délky N=512 a pro odhad na bázi DFT doplňte nulami, jako v předcházejícím kroku (tj. volte NFFT=4096).
% Určete výšku tónu u hudebních signálů z minulého cvičení a srovnejte s výsledky dosaženými pomocí autokorelace - cembalo_d_dur_2.wav, fletna_d_dur_6.wav, housle_d_dur_5.wav, kytara_d_dur_1.wav, piano_d_dur_4.wav, varhany1_d_dur_3.wav, varhany2_d_dur_8.wav. POZN.: Dostupné signály vzorkované kmitočtem 48 kHz bude vhodné převzorkovat minimálně na 16 kHz (pro analyzované signály by vzhledem k frekvenčnímu rozsahu bylo možné i převzorkování na 8 kHz).
%
%% Domácí práce
%% Numerický výpočet FT, FS, DTFT, DFS pomocí DFT
%% Numerický výpočet Fourierovy transformace spojitého neperiodického signálu pomocí DFT.
% Zobrazte v subplot(221) spojitý neperiodický obdélníkový signál s amplitudou A=5, šířkou pulzu to=0.04 s v časovém intervalu -0,1 - 0,4 s. Volte malý krok zobrazení, tj. nejvýše 0.1 ms.
% Zobrazte v subplot(222) též teoreticky vypočtené spojité spektrum, viz doplňkový studijní matediál - DFT-priklady.pdf. Spektrum zobrazte pro úhlový kmitočet omega v pásmu +/-3000 rad/s.
% V subplot(223) zobrazte navzorkovaný konečný počet vzorků uvedeného signálu v časovém intervalu od tmin=0 do tmax=0,4 s. Použijte vzorkovací kmitočet fs=500 Hz.
% V subplot(224) zobrazte vypočtený periodogram pomocí DFT výše zmiňovaného signálu.
% Navzorkovaný signál zobrazte rovněž jako body odlišnou barvou do původního obrázku spojitého průběhu, tj. do subplot(221).
% Spektrum vypočtené pomocí DFT zobrazte pro srovnání rovněž do zobrazení teoreticky vypočteného spojitého neperiodického spektra, tj. subplot(222).
% Opakujte rychle pro fs=820 Hz, fs=200 Hz, fs=2000 Hz.
% Pro fs=820 Hz opakujte pro různé délky vzorkovaného signálu, tj. tmax=0,8 s, tmax=0,08, tmax=0,06, tmax=0,03.
% Srovnejte teoretický výpočet spektra spojitého periodického obdélníkového signálu v časovém intervalu 0-X s a jeho spektra v pásmu +/-YY rad/s.
% Srovnejte teoretický výpočet spektra diskrétního neperiodického obdélníkového signálu vzorkovaného kmitočtem fs=xx Hz v časovém intervalu 0-X s a jeho spektra v pásmu +/-fs Hz.
% Srovnejte teoretický výpočet spektra diskrétního periodického obdélníkového signálu vzorkovaného kmitočtem fs=xx Hz v časovém intervalu 0-X s a jeho spektra v pásmu +/-fs Hz.