%ltfatstart %je třeba mít spuštěn toolbox LTFAT

clc
clear
close all

tic

% Výběr souboru se signálem, určení délky chybějící části

soubor = 'signaly/sinus.wav';
velDir = 430; %počet chybějících vzorků, pozice chybějící části v signálu je náhodná

% Určení parametrů slovníku

velOkna = 2048; %šířka okna
posunOkna = 128; %velikost posunu okna
pocetKanalu = 2048; %počet kanálů

% Určení parametrů Douglas-Rachfordova algoritmu

gama = 0.05; %musí být >0
lambda = 1; %je z intervalu [epsilon,(2 - epsilon)], kdy epsilon je z intervalu ]0,1[
maxIter = 3000; %maximální počet iterací
tolerance = 0.00005; %tolerance rozdílu objektivní funkce pro ukončení iterací

% Parametr spektrogramu

dynRozsah = 80; %dynamický rozsah spektrogramu

% Následuje vlastní výpočet

[y,fs] = audioread(soubor);
[~,n] = size(y);
if n > 1, y = (y(:,1)+y(:,2))/2; end %pokud je nahrávka stereo, vytvoří mono

zacDir = randi([15*velDir,length(y)-16*velDir]); %náhodný první chybějící vzorek
konDir = zacDir+velDir-1; %poslední chybějící vzorek

F = frame('dgtreal',{'hann',velOkna},posunOkna,pocetKanalu); %vytvoření objektu Gaborovy transformace
F = frametight(F); %z F vytvoří 1-těsný frame, pokud už byl těsný, pouze normalizuje
[A,B] = framebounds(F) %zobrazí meze framu, je-li A=B=1 frame je 1-těsný

[ y_rekBezVyl,historieBezVyl,pocetIterBezVyl ] = rekBezVyl( y, F, gama, lambda, maxIter, tolerance, zacDir, velDir );

[ y_rekSVyl,historieSVyl,pocetIterSVyl ] = rekSVyl( y, F, gama, lambda, maxIter, tolerance, zacDir, velDir, velOkna, posunOkna, pocetKanalu );

% Výpočet SNR

SNR_jenDirBezVyl = 20*log10((norm(y(zacDir:konDir),2))/(norm(y(zacDir:konDir)-y_rekBezVyl(zacDir:konDir),2))) %SNR jen v oblasti chybějící části, rekonstrukce bez vylepšení
SNR_celeBezVyl = 20*log10((norm(y,2))/(norm(y-y_rekBezVyl(1:length(y)),2))) %SNR pro celý signál, rekonstrukce bez vylepšení

SNR_jenDirSVyl = 20*log10((norm(y(zacDir:konDir),2))/(norm(y(zacDir:konDir)-y_rekSVyl(zacDir:konDir),2))) %SNR jen v oblasti chybějící části, rekonstrukce s vylepšením
SNR_celeSVyl = 20*log10((norm(y,2))/(norm(y-y_rekSVyl(1:length(y)),2))) %SNR pro celý signál, rekonstrukce s vylepšením

rodilSNR_jenDir = SNR_jenDirSVyl-SNR_jenDirBezVyl
rodilSNR_cele = SNR_celeSVyl-SNR_celeBezVyl

% Zobrazení časových průběhů originálního signálu i obou rekonstrukcí v grafu

figure
plot([y,y_rekBezVyl(1:length(y)),y_rekSVyl(1:length(y))])
title('Původní, zrekonstruovaný signál a zrekonstruovaný signál s vylepšením')
legend('původní signál', 'zrekonstruovaný signál', 'zrekonstruovaný signál s vylepšením')
xlabel('vzorky')
ylabel('úroveň signálu')
line([zacDir-1 zacDir-1],[-1 1],'Color','black','LineStyle',':')
line([konDir+1 konDir+1],[-1 1],'Color','black','LineStyle',':')

toc
cas_min = toc/60 %doba trvání výpočtu v minutách