function [cepst]=vadcepstralni(y,Fvz,delka_ramce,init,p,alfa,graf1,graf2)
%% Soucast detektoru recove aktivity VAD
%funkce pro detekci recove aktivity na principu kepstralniho detektoru
%Priklad spusteni funkce:   vadcepstralni('psenice',256,30,1,0.5,1,1)
%vstupni parametry:
%vst_signal - zvukovy soubor ve formatu *.wav, ktery chceme analyzovat
%delka_ramce - pocet vzorku, vstupniho signalu, ktery se ulozi do jednoho ramce (segmentu)
%volitelne hodnoty(128, 256)
%init - pocet inicializacnich ramcu potrebnych k pocatecnimu nastaveni prahove hodnoty
%volitelne hodnoty(10 - 30)
%
%2008   Michal Adamec VUT Brno (v ramci Diplomove prace - Moderni rozpoznavace reci)


%[y,Fvz]=wavread(vst_sig);                               %y vraci vzorkova data. Amplitudove hodnoty jsou v rozsahu [-1,+1]
signal      = y;                             
delka_okna  = delka_ramce;                              %delka prislusneho typu vahovaciho okenka bez prekryti
prekryti    = delka_okna/2;                             %prekryti oken pri vypoctu
tvar_okna   = 'hamming';
delka_sig   = length(signal);                           %delka vektoru signalu y
cas_delka_ramce = prekryti / Fvz * 1000;                %vypocet casove delky jednoho ramce v [ms] 
p;         %= 1;                                        %casova konstanta exponencialniho prumerovani, viz obnova prahu
t           = (0:delka_sig-1)/Fvz;                      %vektor poctu vzorku signalu y. Je potrebny k vykresleni grafu.       
pocet_ramcu = floor((delka_sig-delka_okna)/prekryti);	%zjisteni poctu ramcu signalu
alfa;      %= 0.5;                                      %konstanta, viz obnova prahu

%*********podminka pro volbu okenkove funkce***************
if strcmp(tvar_okna,'hamming')
    okno = hamming(delka_okna);
else    
    okno = ones(delka_okna);
end


%******inicializace - zjisteni kepstra pozadi z pocatecnich ramcu**
prah_ramcu   = init;                                    %zadavame pocet ramcu, kterych bude vyuzito pro vypocet prahove hodnoty
for i   = 1:prah_ramcu   
    z   = (i-1) * prekryti + 1;  
    e   = abs(signal(z:z+delka_okna-1).* okno);
    h   = length(e);                                    %delka ramce                       
    X   = (sum(e))/h;                                   %prumerna hodnota v ramci              
    c(i) = abs(rceps(X));                               %realne kepstrum prumerne hodnoty
    
end
c0    = abs(sum(c)/prah_ramcu);                         %kepstrum pozadi

%***cyklus pro vypocet vzdalenosti kepstra pozadi a realnych kepster
%vypoctenych v inicializacnich ramcich

for b       = 1:prah_ramcu  
    d(b)    = abs(c0 - c(b));
    var(b)  = (c(b) - c0)^2;                             
end
d           = sum(d)/init;
D           = sum(var)/30;                              %rozptyl

tl = d + alfa*sqrt(D);                                  %prahova hodnota

%***** Vytvoreni casove osy *******************************************
for b=1:pocet_ramcu
    cas(b)      = b * prekryti/Fvz;
end

vystup(init)    = 0;

%******** pocitani v jednotlivych ramcich (oknech s prekrytim)********
for k = init:pocet_ramcu
    s           = (k-1) * prekryti + 1;    	            % zacatek okna
    ramec       = abs(signal(s:s+delka_okna-1).* okno); %vyber hodnot ze signalu, ramec je vazen okenkem

    %**** vypocet kepstralni vzdalenosti realneho kepstra v jednotlivych ramcich od kepstra pozadi
    suma        = sum(ramec);
    g           = length(ramec);
    C           = suma /(g);  
    ce(k)       = abs(rceps(C));                        %vypocet realneho kepstra v jednotlivych ramcich

    de(k)       = abs(c0 - ce(k));                      %kepstralni vzdalenost realneho kepstra v jednotlivych ramcich od kepstra pozadi

%*******detekce****
if de(k) > tl
        vystup  =   1;  
           tl   = tl;
    else
        vystup  = 0;
            
        %***** Obnova prahove hodnoty kepstra *****
        d = p*d + (1-p)*de(k);
        tl = d + alfa*sqrt(D);
        
    end
    vvv(k)= vystup;                                     % vektor vystupni hodnoty 0/1
    prah(k) = tl;
    cepst(k,1) = median([vvv(k-2) vvv(k-1) vvv(k)]);    %vyhlazeni prekmitu pomoci funkce median
   
end


%********graficke vyjadreni****************************************
cara = init * cas_delka_ramce/1000;
konec_inicializace = [-1:0.01:1.1];

% figure
% plot(t,signal,cas,cepst,cara,konec_inicializace,cas,de,cas,tl);
% 
if graf1==1
figure(5)
plot(t,signal,cas,cepst,cara,konec_inicializace,'c'), title('signal a detekce')
xlabel('Cas (s)')
ylabel('Amplituda')
end
if graf2==1
figure(6)
subplot(211)
plot(t,signal,cas,cepst,cas,de,'r',cas,prah,'k',cara,konec_inicializace,'c')
xlabel('Cas (s)')
ylabel('')
legend('signal', 'detekce', 'kepstralni vzdalenost', 'prahova hodnota','konec inicializace') 

subplot(212)
plot(cas,de,'r',cas,prah,'k')
xlabel('Cas (s)')
ylabel('')
legend('kepstralni vzdalenost', 'prahova hodnota')
end