sift(X,Fs,oldV,oldF,B,aff)

% [T0,V] = sift(X) cette fonction calcule par tranche de 30 ms la valeur de To du vecteur
% acoustique X et indique si la tranche analyse faisait partie d'un son vois ou pas.

function [F0,V,aff] = sift(X,Fs,oldV,oldF,B,aff)

Y = filter (B,1,X);    % filtrage d'entre
E = decimate (Y,4);    % sous-chantillonnage

%*****************************************************************************************%
%                                                                                         %
%                 Calcul des coefficients du filtre inverse                               %
%                                                                                         %
%*****************************************************************************************%

nu = 0.95;            % Praccentuation
B2 = [1 -nu];
Y2 = filter (B2,1,E);

Y3 = Y2.*hamming (length(Y2));

coef = LPC (Y3,4);    % Calcul des coefficients PARCOR
rf=lpcar2rf(coef);

%*****************************************************************************************%
%                                                                                         %
%                         Filtrage par le filtre inverse                                  %
%                                                                                         %
%*****************************************************************************************%

E2 = filter (rf,1,E);

%*****************************************************************************************%
%                                                                                         %
%                             Calcul de la corrlation                                    %
%                                                                                         %
%*****************************************************************************************%

absE2 = abs(E2);    % Calcul du seuil et slection des valeurs plus grandes que ce seuil
maxE2 = max (absE2);
seuil = 0.6*maxE2;
len = length (E2);
R = zeros (len);

for i=1:len,
   if absE2(i)>seuil
        R(i)=E2(i)- sign(E2(i))*seuil;
   end;
end;

C = xcorr(R);    % Calcul de l'autocorrlation

%*****************************************************************************************%
%                                                                                         %
%                 Chercher le max de l'autocorrlation entre 60 Hz et 320Hz.              %
%                                                                                         %
%*****************************************************************************************%

per1 = 4/Fs;                    % Calcul de la priode correspondant  1'intervalle
Binf = floor (1/(400*per1));    % Calcul du nombre d'intervalles correspondant  une priode de 1/400Hz
Bsup = floor (1/(60*per1));     % Calcul du nombre d'intervalles correspondant  une priode de 1/60Hz

limit = len+Bsup;              % Calcul du nombre max d'intervalles  prendre en compte
Tlimit=length (C);              % en prvoyant le cas o Bsup serait plus grand que la taille du vecteur

if limit > Tlimit
    bornesup = Tlimit;
else
   bornesup = limit;
end;

Cxx = abs(C(len+Binf:len+Bsup));

[Amax,Imax]= max (Cxx);         % Localisation du max.
Cxx(Imax-1:Imax+1)=0;           % Mise  zro du 1er max(max +environs)afin de localis le 2nd max
[Amax2,Imax2]= max (Cxx);       % Localisation du second max

% Vrification de l'historique des prcdents vecteurs afin de prendre le bon maximum

if (oldV(2)*oldV(1) == 0) & ((Imax/Imax2)> 1.6) & (Amax2 > (0.6*Amax))
    Imax = Imax2;
    Amax = Amax2;
elseif (abs(oldF(2)-oldF(1))< 0.2*oldF(2)) & (abs(oldF(2)-(1/(per1*Imax2)))< abs (oldF(2)-(1/(per1*-Imax)))) & ((Imax/Imax2)> 1.2) & (Amax2 > (0.4*Amax))
    Imax = Imax2;
    Amax = Amax2;
end;

Imax = Imax + Binf;
%*****************************************************************************************%
%                                                                                         %
%     Interpolation quadratique entre les 2 max. Sinon la prcision serait de 100Hz       %
%                                                                                         %
%*****************************************************************************************%

a = 1/2*((C(Imax+1)+C(Imax-1))-2*C(Imax));
b = 1/2*(C(Imax+1)-C(Imax-1));
Xmax = b/(2*a);

if abs(Xmax)<1                 % Xmax plausible donc modifier
   Imax = Imax+Xmax;
    Amax = Amax + (b*b)/(4*a);
end;

% Xmax pas possible donc on ne modifie pas

Arel = Amax/C(len);

%*****************************************************************************************%
%                                                                                         %
%                                 Seuil de dcision                                       %
%                                                                                         %
%*****************************************************************************************%

if len < 52            % Dtermination des valeurs de seuil selon les cas
    So = 0.248;        % Ce sont des valeurs empiriques
    S1 = 0.8;
    FS2 = 0.6;
else
   So = 0.348;
    S1 = 0.9;
    FS2 = 0.6;
end;

if Imax > 16
    seuil = So-0.003*Imax;
else
   seuil = S1-0.0375*Imax;
end;

seuil2 = seuil*FS2;
F0 = 0;

if Arel>seuil        % Test du maximum
    F0 = 1/(per1*Imax);
elseif oldV(2)== 1
   if Arel>seuil2
        F0 = 1/(per1*Imax);
   end;
end;

if F0 == 0
    V = 0;
else
   V = 1;
end;

%***********************************************************************************%
%                                                                                   %
%         Affichage des diffrentes tapes du traitement de la fentre              %
%                                                                                   %
%***********************************************************************************%

if aff == 1
    disp ('Imax interpol vaut :');
    disp(Imax);
    disp ('Pitch');
    disp(F0);
    
    figure;
    subplot (4,2,1);
    plot (X);
    title ('Fentre non traite');
    subplot (4,2,2);
    plot (Y);
    title ('Filtre passe-bas');
    subplot (4,2,3);
    plot (E);
    title ('Dcimation');
    subplot (4,2,4);
    plot (Y2);
    title ('Pr-amplification');
    subplot (4,2,5);
    plot (C(len+Binf:len+Bsup));
    title ('Autocorrlation');
    subplot (4,2,6);
    plot (rf,'p');
    title ('Coefficients de PARCOR');
    subplot (4,2,7);
    plot (E2);
    title ('Rsidu du LPC4');
    subplot (4,2,8);
    plot (abs(C(len+Binf:len+Bsup)),'b');
    title ('Dcision');
    aff = 0;
    pause;
end;