function [intenz]=vadintenz(vst_sig,delka_ramce,init,p,alfa,graf1,graf2)
 
%funkce pro detekci recove aktivity na principu kratkodobe intenzity
%Priklad spusteni funkce:   vadintenz('psenice',256,30,0,0.1,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;         %= 0.3;                                     %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, ve kterych se pocita energie
alfa;       %= 1;                                       %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 prahove hodnoty energie 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);                                       
    Mdl(i)  = (sum(e))/h;                               %vypocet prahove hodnoty kratkodobe intenzity pro jednotlive ramce
end

Md = sum(Mdl) / prah_ramcu;                             %prumerna hodnota prahove intenzity

%*******vypocet pocatecni prahove intenzity ************************
for x = 1:prah_ramcu 
    var(x) = (Mdl(x) - Md)^2;
end

Di   = sum(var);                                          %rozptyl
tl  = Md + alfa*sqrt(Di);                                 %pocatecni prahova hodnota intenzity


%**** Nastaveni pocatecni prahove intenzity v prahovych ramcich *****
for i = 1:prah_ramcu
    Mp(i) = tl;
end

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

%******** 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 prumerne intenzity v ramci
    suma        = sum(ramec);
    g           = length(ramec);
    prumer      = suma /(g);  
    Mn(k)       = prumer;                               % prumerna intenzita v ramci
    %**************************************
    Mp(k)       = tl;                                   % pocatecni prahova hodnota intenzity

%******** podminka pro detekci reci. Porovnavame prahovou hodnotu intenzity, s prumernou intenzitou **
    if Mn(k) > Mp(k)
        vystup  =   1;  
            Md  =   Md;
    else
        vystup  = 0;
            %***** Obnova prahove hodnoty intenzity *****
            Md  = (1-p)*Md + p*Mn(k);
            tl  = Md + alfa*sqrt(Di); 
        
    end
    vv(k,1)=vystup;                                      % vektor vystupni hodnoty 0/1
    intenz(k,1) = median([vv(k-4) vv(k-3) vv(k-2) vv(k-1) vv(k)]);             %vyhlazeni prekmitu pomoci funkce median
     
end

%********graficke vyjadreni****************************************

%normalizce vektoru Mn, aby nejvetsi hodnota intenzity byla rovna 1
Mnmax = max(Mn);                                      %zjisteni max hodnoty y vektoru energ
for j = 1:pocet_ramcu   
     normaliz(j) = Mn(j)/Mnmax;                       %provedeni normalizace
end

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

if graf1==1
figure(3)
plot(t,signal,cas,intenz,cara,konec_inicializace,'c'), title('signal a detekce')
xlabel('Cas (s)')
ylabel('Amplituda')
end

if graf2 == 1
figure(4)

subplot(211)
plot(t,signal,cas,intenz,cas,normaliz,'r',cas,Mp,'k',cara,konec_inicializace,'c')
xlabel('Cas (s)')
ylabel('')
legend('signal', 'detekce', 'intenzita', 'prahova int','konec inicializace') 

subplot(212)
plot(cas,Mn,'r',cas,Mp,'k')
xlabel('Cas (s)')
ylabel('')
legend('intenzita', 'prahova int')
end
