function [v]=g729(vst_sig,Bef,Bel,BZC,Blsf,graf)

%funkce pro detekci recove aktivity
%Priklad spusteni funkce:   g729('franek',0.1,0.1,0.1,0.1,1)
%vstupni parametry:
%vst_signal - zvukovy soubor ve formatu *.wav, ktery chceme analyzovat
%Bef - autoregresivni koeficient pro aktualizaci klouzaveho prumeru, hodnoty 0 az 1
%Bel - autoregresivni koeficient pro aktualizaci klouzaveho prumeru, hodnoty 0 az 1
%BZC - autoregresivni koeficient pro aktualizaci klouzaveho prumeru, hodnoty 0 az 1
%Blsf - autoregresivni koeficient pro aktualizaci klouzaveho prumeru, hodnoty 0 az 1
%
%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]
%%%%%%%%% pocatecni inicializace %%%%%%%%%%%%
signal  = y;
window  = 240;                                              %delka prislusneho typu vahovaciho okenka bez prekryti
frame_past = 120;
frame_present = 80; %delka
frame_future = 40;  %prekryti
r=0;
xr = 0;
IVD = 0;
R= 0;
pocitadlo = 0;
Efp = 0;

% pro vyhlazeni
S1vd = 1;          % pro inicializaci tretiho kroku, vyhlazene rozhodnuti o aktivite v minulem ramci
S2vd = 1;          % vyhlazene rozhodnuti o aktivite v predminulem ramci
FVD1 = 1;          % logicky parametr
Ce = 0;            % pocitadlo vyhlazovaciho procesu
Cs = 0;
Cn = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%TABULKA KONSTANT (itu.G729)%%%%%%%%%%%%%%%%%%%%%
% /2^15 a /2^31 prevod z 32-bit a 16-bit cisla
Ni = 32;
N0 = 128;
K0 = 0;
K1 = -53687091/2^31;
K2 = -671088640/2^31;
K3 = -93952410/2^31;
K4 = -134217728/2^31;
K5 = -161061274/2^31;
a1 = 23488/2^15;
a2 = -30504/2^15;
a3 = -32768/2^15;
a4 = 26214/2^15;
a5 = 0;
a6 = 28160/2^15;
a7 = 0;
a8 = 16384/2^15;
a9 = -19065/2^15;
a10 = 0;
a11 = 22400/2^15;
a12 = 30427/2^15;
a13 = -24576/2^15;
a14 = 23406/2^15;
N1 = 4;
N2 = 10;
T1 = 671088640/2^31;
T2 = 738197504/2^31;
T3 = 26843546/2^31;
T4 = 40265318/2^31;
T5 = 40265318/2^31;
T6 = 40265318/2^31;
b1 = 28521/2^15;
b2 = 19446/2^15;
b3 = -32768/2^15;
b4 = -19661/2^15;
b5 = -30802/2^15;
b6 = -19661/2^15;
b7 = 30199/2^15;
b8 = -22938/2^15;
b9 = -31576/2^15;
b10 = -173671/2^15;
b11 = -27034/2^15;
b12 = 29959/2^15;
b13 = -29491/2^15;
b14 = -28087/2^15;

posun       = window/3;        % 80                        %posun oken pri vypoctu 80 vzorku
delka_sig   = length(signal);                              %delka vektoru signalu y
cas_delka_ramce = posun / Fvz * 1000;                      %vypocet casove delky jednoho ramce v [ms] 10ms
t = (0:delka_sig-1)/Fvz;                                   %vektor poctu vzorku signalu y. Je potrebny k vykresleni grafu.

pocet_window = floor((delka_sig-window)/window);  % 78     % pocet oken, pro vazeni nesymetrickym oknem
pocet_posunu = floor((delka_sig)/posun);    %234           %zjisteni poctu ramcu, ve kterych se pocita energie

%filtrovani horni propusti s fm = 140 Hz
b=[0.46363718 -0.92724705 0.46363718];
a=[1 -1.9059465 0.9114024];
signal=filter(b,a,signal);                                  

oknoa = 0.54-0.46*cos((2*pi*([1:200]-1))/399);
oknob = cos((2*pi*(([201:240]-1)-200))/159);
okno = [oknoa oknob];
clear oknoa;
clear oknob;
clear b;
clear a;
sloupcu = floor(delka_sig/80);

%%% rozdeleni signalu na sloupce po 80 vzorcich
for j=1:sloupcu  %1:234
    jj  = (j-1) * 80 + 1;                                   % posunujeme se po 80 radcich
    r=0;
    for i=jj:jj+80-1
        r=r+1;
        signal_sl(r,j)=signal(i,1);
    end
end

%%% prochazeni signalu po 80 vzorcich (10ms), vazeni signalu 240 vzorku
for k = 1:pocet_posunu-3                                    % -3 aby nedoslo k preteceni pres pocet vzorku

    %%%%%%% rozdeleni na framy (past,present,future) - vytvoreni LP
    %%%%%%% windows
    if k>2 && k<pocet_posunu
        LP_window=0;
        %%%% past
        LP_window(1:40) = signal_sl(41:80,k-2);
        LP_window(41:120) = signal_sl(:,k-1);
        %%%% present
        LP_window(121:200) = signal_sl(:,k);
        %%%% future
        LP_window(201:240) = signal_sl(1:40,k+1);

        LP_window = LP_window.*okno; % radek x radek
        LP_window = LP_window';                            %vysledny ramec po nasobeni okenkovou funkci
    else                                                   %%% dokud neni k>2, nemuzu ziskat past frame, neni z ceho
        LP_window=0;
        LP_window(1:80) = signal_sl(:,k);
        LP_window(81:160) =   signal_sl(:,k+1);
        LP_window(161:240) = signal_sl(:,k+2);
        %LP_window =signal(:,1);        alternativa
        LP_window = LP_window.*okno;
        LP_window = LP_window';
    end
    length(LP_window);
    length(okno);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% VYPOCET SPEKTRALNICH PARU LSF %%%%%%%%%%%%% autokorelace %%%
    xr = xcorr(LP_window);                      % autokorelace, koeficienty udavaji podobnost sign. samotneho se sebou, kdyz ho posuneme o k vzorku (k=10)
    wlag=exp(-0.5*((2*pi*60*[1:10])/Fvz).^2);
    wlagf=exp(-0.5*((2*pi*60*[1:12])/Fvz).^2);
    % pro filtraci
    xrf = [1.0001*xr(240) wlagf.*xr(241:252)'];
    xr = [1.0001*xr(240) wlag.*xr(241:250)'];

    LP = levinson(xr);
    LSF(:,k) = poly2lsf(LP);
    energie(:,k) = LP_window' * LP_window;
    Ef(:,k) = 10*log10(energie(:,k));          %  sum(Ef)/length(Ef)*2

    %%%%%%%%%% Uzkopasmova energie El %%%%%%
    % pomoci fdatool 5.3. podle jineho zdroje DP 1kHz: fdatool: LP,
    % butterworth, Least-squares, order 13, fs = 8000, Fpass=400, Fstop = 1000
    h = [
        -0.035206032388628
        -0.022452379503980
        0.010602845885909
        0.061167885940769
        0.119217524826979
        0.200085319792819
        0.200085319792819
        0.170236577633659
        0.119217524826979
        0.061167885940769
        0.010602845885909
        -0.022452379503980
        -0.035206032388628
        ];
    R=toeplitz(xrf);   % toeplitzova matice pro matici 13 x 13
    %%%%%%%%%%
    El(:,k) = (10*log10(1/240*h'*R*h));

    %%%%%%% Střední počet průchodů signál nulou %%%%%%%%%%%%
    % pocitam vzdy v aktualnim ramci LP_window(121:200)
    M = 80;
    sgn = (abs(sign(LP_window(120+[1:M]))-sign(LP_window(120+[1:M]-1))));
    ZC(:,k) = 1/(2*M)*sum(sgn);                                          % normalizovana hodnota stredniho poctu pruchodu signalu nulou

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%INICIALIZACE KLOUZAVÉHO PRŮMĚRU PARAMETRŮ Z CHARAKTERISTIK%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%ŠUMOVÉHO POZADÍ%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%prvni krok inicializace%%%%
    % pro prvnich Ni ramcu:
    if k < Ni
        % pocatecni nastaveni En, ZCp, LSFp
        LSFp(:,k) = LSF(:,k);                       %klouzavy prumer spektralnich parametru
        ZCp(:,k)  = ZC(:,k);                        %klouzavy prumer stredniho poctu pruchodu nulou
        En(:,k)   = Ef(:,k);
        pocitadlo = pocitadlo+1;
        if Ef(:,k) > (-15)                          % energie ramce je vetsi nez 15dB viz
            IVD(1,k) = 1;
        else
            IVD(1,k) = 0;
        end
    end

    % prumer LSF,ZC a En v Ni ramcich
    if k == Ni
        LSFp = sum(LSFp')./pocitadlo;
        ZCp  = sum(ZCp)   /pocitadlo;
        En   = sum(En)    /pocitadlo;

        %%%%%druhy krok inicializace%%%%
        %%%%% inicializace prumeru energii Elp, Efp
        if En < T1
            Efp = En + K0;
            Elp = En + K1;    
        elseif  T1 <= En <= T2
            Efp = En + K2;
            Elp = En + K3;
        else
            Efp = En + K4;
            Elp = En + K5;
        end
    end

    %%%%%%%%%%%% vypocet rozdilovych parametru %%%%%%%%%%%%%%
    if k > Ni
        delta_S(1,k) = sum((LSF(:,k)-LSFp').^2);
        delta_Ef(1,k) = Efp - Ef(:,k);
        delta_El(1,k) = Elp - El(:,k);    % vychazi nejak divne
        delta_ZC(1,k) = ZCp - ZC(:,k);

        if Ef(:,k) > (-15)      % kdyz je Ef vetsi nez 15 dB, projde vicebodovym rozhodovacim pravidlem
            % jinak Ef je rovno 0
            %%%%%%%%%%% Vicebodove rozhodovaci pravidlo %%%%%%%%%%%
            IVD(1,k) = 0;           % pri inicializci je IVD = 0
            %1
            if delta_S(1,k) > a1 * delta_ZC(1,k)+b1
                IVD(1,k) = 1;
                %2
            elseif delta_S(1,k) > a2 * delta_ZC(1,k)+b2
                IVD(1,k) = 1;
                %3
            elseif delta_Ef(1,k) < a3 * delta_ZC(1,k)+b3
                IVD(1,k) = 1;
                %4
            elseif delta_Ef(1,k) < a4 * delta_ZC(1,k)+b4
                IVD(1,k) = 1;
                %  5
            elseif delta_Ef(1,k) < b5
                IVD(1,k) = 1;
                %6
            elseif delta_Ef(1,k) < a6 * delta_S(1,k)+b6 
                IVD(1,k) = 1;
                % 7
            elseif delta_Ef(1,k) > b7
                IVD(1,k) = 1;
                %8
            elseif delta_El(1,k) < a8 * delta_ZC(1,k)+b8    
                IVD(1,k) = 1;
                %9
            elseif delta_El(1,k) < a9 * delta_ZC(1,k)+b9
                IVD(1,k) = 1;
                %10
            elseif delta_El(1,k) < b10
                IVD(1,k) = 1;
                %11
            elseif delta_El(1,k) < a11 * delta_S(1,k)+b11
                IVD(1,k) = 1;
                %12
            elseif delta_El(1,k) < a12 * delta_Ef(1,k)+b12
                IVD(1,k) = 1;
                %13
            elseif delta_El(1,k) < a13 * delta_Ef(1,k)+b13
                IVD(1,k) = 1;
                %14
            elseif delta_El(1,k) < a14 * delta_Ef(1,k)+b14
                IVD(1,k) = 1;
            end
        else
            IVD(1,k)=0;
        end
        cass(k)= k * posun/Fvz;    % pro graf

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%Vyhlazeni rozhodnuti VAD %%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if k>2
            v_flag = 0;                                     % indikuje aktivitu vyhlazovaciho procesu
            S0vd(1,k) = IVD(1,k);                           % vyhlazene rozhodnuti o aktivite v aktualnim ramci
            S1vd(1,k) = IVD(1,k-1);
            S2vd(1,k) = IVD(1,k-2);

            %%%%%%%% prvni krok %%%%%%%%%%
            if IVD(1,k) == 0 && S1vd(1,k) == 1 && Ef(1,k) > (Efp + T3)
                S0vd(1,k) = 1;
                v_flag = 1;
            end

            %%%%%%%% druhy krok %%%%%%%%%
            if (FVD1 == 1) && (IVD(1,k) == 0) && (S1vd(1,k) == 1) && (S2vd(1,k) == 1) && (abs(Ef(1,k) - Ef(1,k-1)) <= T4)
                S0vd(1,k) = 1;
                v_flag = 1;
                Ce = Ce +1;

                if (Ce <= N1)
                    FVD1 = 1;
                else
                    FVD1 = 0;
                    Ce = 0;
                end

            else FVD1 = 1;
            end

            %%%%%%% treti krok %%%%%%%%
            if (S0vd(1,k) == 0)
                Cs = Cs + 1;
            end

            if (S0vd(1,k) == 1) && (Cs > N2) && ((Ef(1,k) - Ef(1,k-1)) <= T5)
                S0vd(1,k) = 0;
                Cs = 0;
            end

            if (S0vd(1,k) == 1)
                Cs = 0;
            end

            %%%%%%% ctvrty krok %%%%%%%
            if (Ef(1,k) < (Efp + T6)) && (k > N0) && (v_flag == 0)
                S0vd(1,k) = 0;
            end
            IVD(1,k)=S0vd(1,k);

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%moje vyhlazeni%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            vvv(k,1)=IVD(1,k);                                              % vektor vystupni hodnoty 0/1
            vv(k,1) = round(median([vvv(k-3) vvv(k-2) vvv(k-1) vvv(k)]));   % vyhlazeni prekmitu pomoci funkce median
            v(k-2,1) = vv(k,1);
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%% Aktualizace klouzaveho prumeru %%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if k > N0
            Cn = Cn + 1;   % pocitadlo aktualizaci prumeru
            Emin = min(Ef(1,k-N0:k));
            zEmin(1,k) = Emin;
        end
        if (Ef(1,k) < (Efp + T6))
            Efp = Bef(1) * Efp + ((1-Bef(1))*Ef(1,k));
            Elp = Bel(1) * Elp + ((1-Bel(1))*El(1,k));
            ZCp = BZC(1) * ZCp + ((1-BZC(1))*ZC(1,k));
            LSFp = Blsf(1) * LSFp + ((1-Blsf(1))*LSF(:,k)');
        end
        if (k > N0) && (Efp < Emin)
            Efp = Emin;
            Cn = 0;
        end
    end
end

v(length(cass)-1,1)=0;
v(length(cass),1)=0;


if graf==1
    figure(7)
    plot(t,signal,cass,v)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


