%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.000005; %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_rek,historie,pocetIter ] = rekBezVyl( y, F, gama, lambda, maxIter, tolerance, zacDir, velDir );

% Výpočet SNR

SNR_jenDir = 20*log10((norm(y(zacDir:konDir),2))/(norm(y(zacDir:konDir)-y_rek(zacDir:konDir),2))) %SNR jen v oblasti chybějící části
SNR_cele = 20*log10((norm(y,2))/(norm(y-y_rek(1:length(y)),2))) %SNR pro celý signál

% Zobrazení výsledků v grafech

figure(1)
plot([y,y_rek(1:length(y))])
title('Původní a zrekonstruovaný signál')
legend('původní signál', 'zrekonstruovaný signál')
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',':')

figure(2)
plot(historie(1:pocetIter))
title('Historie objektivní funkce v průběhu iterací')
xlabel('počet iterací')

figure(3)
plotframe(F,frana(F,y),'dynrange',dynRozsah)
title('Spektrogram původního signálu')
xlabel('vzorky')
ylabel('frekvence (normalizována)')
xlim([zacDir-3*velDir konDir+3*velDir])
line([zacDir-1 zacDir-1], get(gca,'YLim'), 'Color','black','LineStyle',':')
line([konDir+1 konDir+1], get(gca,'YLim'), 'Color','black','LineStyle',':')

figure(4)
plotframe(F,frana(F,y_rek),'dynrange',dynRozsah)
title('Spektrogram zrekonstruovaného signálu')
xlabel('vzorky')
ylabel('frekvence (normalizována)')
xlim([zacDir-3*velDir konDir+3*velDir])
line([zacDir-1 zacDir-1], get(gca,'YLim'), 'Color','black','LineStyle',':')
line([konDir+1 konDir+1], get(gca,'YLim'), 'Color','black','LineStyle',':')

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