Détermination du pitch par la méthode du SIFT

1. Problème de la détermination du pitch

    Le problème de la détermination du pitch est délicat dans le cas d'un signal vocal. En effet, ce dernier est un signal non-stationnaire, c'est-à-dire qu'il ne se reproduit pas de manière régulière, aucune période du signal n'est totalement identique à la précédente. Par exemple, il n'est pas question de calculer plusieurs fois la période du signal puis de faire une moyenne des résultats obtenus pour en obtenir une estimation plus précise. De plus, le signal vocal est constitué de zones voisées et non-voisées, les premières sont obtenues à partir d'une vibration "régulière" des cordes vocales. Ce sont donc des zones où il existe une fréquence fondamentale. Les deuxièmes sont, elles, obtenues sans faire intervenir les cordes vocales, à partir d'un écoulement turbulent de l'air dans la bouche, telle une source de bruit.  Il est donc illusoire d'espérer trouver une fréquence fondamentale dans une telle zone. Il est impératif de prendre cela en considération sous peine de trouver des valeurs de pitch totalement incohérentes, en forçant l'algorithme a trouver des valeurs de pitch où il n'y en a pas. Un dernier problème qui est à la fois un désavantage et un avantage c'est que la fréquence fondamentale est en constante évolution, nous retombons donc sur le problème du début. Mais c'est aussi un avantage car il n'évolue pas de manière brusque, les valeurs précédentes ont donc de l'influence sur les valeurs précédentes. La prise en compte de cette caractéristique va donc nous permettre d'amélioré notre estimation du pitch et d'écarter les valeurs aberrantes.

    Je vais donc ici vous présenter une technique assez complète de détermination du pitch : le SIFT. Bien que ne donnant pas de résultats parfait, cette dernière prend en compte assez de paramètres que pour donner de bons résultats....  

2. Méthode du SIFT

2.1. Principes de cette méthode

 La méthode du SIFT (Simplified Inverse Filtering) se base sur une estimation du maximum de l'autocorrélation d'un signal filtré au préalable par le filtre LPC inverse.

Pourquoi c'est deux techniques combinées?

    Premièrement il est logique d'utilisé l'autocorrélation pour trouver la période de répétition, étant donné que cette dernière multiplie entre eux les échantillons de la fenêtre décalés entre eux. Il est donc logique de retrouvé un maximum de l'autocorrélation pour un décalage correspondant à la période fondamentale du signal. On peut donc ainsi facilement déterminer le nombre d'échantillons correspondant à la période du signal.

    Maintenant que le problème de la détermination du pitch à partir d'un signal acoustique est réglé, reste à réglé un dernier problème lié au signal vocal  lui-même. En effet, la composition spectrale du signal vocal est composé de formants qui peuvent entraîné l'accentuation de fréquences autres que la fréquence fondamentale, ce qui "salit" en quelque sorte le signal duquel on veut retrouvé le pitch. On peut d'ailleurs remarqué sur ce graphique que dans ce cas-ci le premier à effectivement tendance à cacher la fréquence fondamentale.

  

     C'est pourquoi  la méthode du SIFT procède comme son nom l'indique à un filtrage inverse du signal, dont le but est d'atténuer l'influence des formants. En effet, grâce à une analyse LPC on peut retrouvé les coefficients du filtre correspondant aux formants du signal, si l'on inverse ce filtre on annule donc l'effet des formants sur le signal.

    Nous sommes donc avec cette méthode dans les conditions idéales pour réaliser une recherche du pitch du signal.

2.2. Réalisation pratique  

    Voici donc le schéma de la réalisation pratiquer du SIFT :

    

      On remarque d'emblée la présence d'un filtre d'un filtre passe-bas, dont le but est encore une fois d'atténuer l'influence des fréquences autres que la fondamentale. Le filtre coupe à 1000Hz étant donné qu'une valeur de pitch est comprise entre environ 80 et 450 Hz. Mais son rôle le plus important reste tout de même d'éviter les recouvrements de spectres lors du sous-échantillonnage qui suit le filtre. En effet le théorème de Shannon nous dit qu'il faut pour préserver la qualité d'un signal échantillonné à au moins deux fois la fréquence maximale du spectre du signal. Or si nous arrivons avec un signal échantillonné à 8000 Hz, dont la fréquence maximale de spectre sera de 4000Hz, on le sous échantillonne à 2000 Hz, la fréquence maximale que l'on peut garder dans le signal sera de 1000 Hz sinon on aura un recouvrement de spectre qui va rendre le signal totalement inutilisable... Le but de ce sous-échantillonneur est simplement de diminuer la charge de calcul, en diminuant le nombre de points sur lesquels se porte l'analyse. Ensuite comme prévu, on retrouve une analyse LPC (limitées à 4 pôles, encore une fois pour la charge de calcul), avec, bien entendu, une pré-accentuation et une pondération de la fenêtre par une fenêtre de Hamming. Cette pondération ayant pour but d'amoindrir les fortes variations du signal sur les bords de la fenêtre, variations qui entraînent une mauvaise estimation des coefficients du filtre si elles n'étaient pas atténuées. On comprends facilement le but de l'opération en observant l'exemple ci-dessous.

    On voit donc bien que les variations sur les bords de la fenêtres ont été nettement diminuées.

     Une fois ces pré-traitements effectués, on effectue l'analyse LPC et on passe les coefficients PARCOR  ainsi obtenus au filtre inverse, le signal à la sortie de celui-ci est donc débarrassé, en partie, de l'influence des formants, comme déjà expliqué plus haut. Ensuite, comme prévu on calcule le vecteur d'autocorrélation de la fenêtre afin d'obtenir un maximum pour un décalage équivalent à la période du signal. Afin d'encore accentué ce maximum, nous allons encore traité le signal avant de calculer l'autocorrélation. Le traitement appliqué est simple, pour retrouver la période d'un signal on compte le temps écoulé entre 2 maxima ou minima du signal, nous n'allons donc garder dans ce dernier que les valeurs dépassant un certains seuil, les autres sont mises à zéros. Le vecteur de corrélation ainsi obtenu présentera des pics pour des valeurs de décalage multiple de la période du signal et aura une valeur nulle pour les autres décalages.

    Nous pouvons donc maintenant calculer le vecteur d'autocorrélation, vecteur duquel on extrait le maximum. Malheureusement, le positionnement de ce dernier est trop imprécis pour être utilisable. En effet, afin de gagner en temps de calcul, nous avons sous échantillonné le signal, en l'occurrence nous n'avons garder qu'un échantillon sur 4, donc, si la fréquence d'échantillonnage de base était de 8000 Hz, nous obtenons donc une fréquence d'échantillonnage de 2000 Hz. Or, le positionnement du maximum que nous obtenons se fait par valeur entière de la période d'échantillonnage (le signal se répète tous les X échantillons), ce qui entraînerait une précision de 2 Hz pour des petites valeurs du pitch (80 Hz) et 77 Hz pour les grandes valeurs du pitch (400 Hz). Ce qui reviendrait donc à perdre tous les efforts consentis jusqu'à présent pour obtenir une valeur la plus précise possible du pitch. C'est pourquoi on essaie d'interpoler la position du maximum, en supposant que l'évolution de la courbe d'autocorrélation est quadratique entre deux points. Cette interpolation nous permet donc de rattraper le manque de précision du positionnement du maximum tout en gardant l'avantage du sous-échantillonnage...

    La dernière étape consiste à prendre compte du fait que certaines valeurs du pitch sont plus probable que d'autres et que l'autocorrélation ne cible pas aussi bien la période qu'on ne pourrait le croire, en effet on peut très bien se trouver avec un maximum important correspondant à une fréquence fondamentale de 350 Hz (valeur possible mais peu probable) et un second maximum plus faible correspondant lui à un fondamental de 200 Hz. On serait donc tenté de choisir comme pitch la valeur de 400 Hz alors qu'il est fort probable que ce soit plutôt 200 Hz qui soit la bonne valeur. Afin d'éviter ces erreurs, on ajoute un seuil de décision variant avec la fréquence, c'est à dire un seuil bas pour les valeurs les plus fréquentes et augmentant avec la fréquences pour atteindre son maximum pour les valeurs les moins probables! ce seuil de décision nous permets donc de déterminer, avec plus de sécurité, la valeur du pitch. Une fois cette valeur valeur obtenue, on effectue un dernier test qui consiste à vérifier qu'il est compris dans des valeurs possibles (entre 80 et 400Hz), s'il passe ce test, on sort la valeur du pitch et on signale que c'est une fenêtre voisée (V); sinon on met la valeur de F0 à zéro et on signale que la fenêtre était non-voisée (NV).

    On se rend donc déjà compte, avec le SIFT, de la complexité de la détermination de la fréquence fondamentale d'un signal vocal. D'ailleurs, comme dans beaucoup de domaines difficilement maîtrisés par l'ingénieur, il existe des règles empiriques. C'est le cas ici pour les valeurs de seuils de décision finale.   

  

 

3. Explication sommaire du programme SIFTECH et SIFT

    La dernière étape de ce projet  était la réalisation de d'un fichier m sous Matlab réalisant l'algorithme du SIFT. Comme vous le verrez, il y a encore quelques différences avec les principes présentés ci-dessus, mais l'idée générale de la réalisation reste la même.

    Pour faciliter les choses, le programme a été séparé en deux parties : la première "SIFTECH" qui s'occupe de calculer les coefficients du filtre d'entrées, de prendre l'ensemble du signal et de le découper en fenêtre de 30 ms qui sont données à traité au deuxième programme "SIFT" qui ,lui ,s'occupe du SIFT à proprement dit. L'affichage des résultats étant géré en partie par chaque programme. Le premier programme permet la visualisation des résultats généraux : valeurs du pitch, voisé/non voisé, affichage du signal de base. Il permet aussi de définir sur quel fenêtre on veut avoir des informations plus précises quant au traitement : évolution du signal à travers les différents traitements.

    Il faut encore savoir que pour réaliser ce programme sous Matlab, j'ai fait appel à une toolbox non inclue dans les toolbox de base : la Voicebox. Cette dernière étant disponible sur ce site

 3.1. SIFTECH

     Voici, le code utilisé pour réaliser cette fonction somme toute assez simple :

function [F0s,Vs]=siftech (Y,Fs,fen)

close all;

Y = Y/max(abs(Y));       % normalise le vecteur d'entrée
N = 10e-3*Fs;             % Calcule le nombre d'échantillon à prendre avoir une fenêtre de 10ms

% Calcule les coefficients du filtre d'entrée

[n,fo,mo,w] = remezord ([900 1300],[1 0],[0.01 0.01],Fs);
B = remez (n,fo,mo,w);

% Initialisation des variables utiles au programme

long = floor(length (Y)/N);    % Détermine le nombre de fenêtres de 10 ms disponibles

F0last = [0 0];                 % Initialise des vecteurs utiles à la fonction SIFT
Vlast = [0 0];

F0s = [0];                      % Initialise le vecteur des fréquences fondamentales
Vs = [0];                       % Même chose pour le voisé/non voisé

af = 0;                         % Mise à zéro du flag demandant l'affichage des détails
                                % du traitement à la fonction SIFT

fen = fen +1;                   % Adaption du choix de la fenêtre dont il faut afficher les détails
                                % pour qu'il corresponde aux valeurs utilisées dans la boucle

for i=2:(long-2),               % Boucle permettant la création de fenêtres de 30 ms
  
  
  if i == fen                 % Test permettant d'afficher le détail de la bonne fenêtre
       
af = 1;
   end;

ech = Y ((i-1)*N:(i+2)*N);    % Création des fenêtres de 30 ms
[F0,V,af]= sift (ech,Fs,Vlast,F0last,B,af);    % Appel de la fonction SIFT

F0s = [F0s F0];                % Création des vecteurs de sorties
Vs = [Vs V];

Vlast = Vs (i-1:i);            % Création des vecteurs "historique"
F0last = F0s (i-1:i);          % nécessaire à la fonction SIFT

end
;                           % Fin de la boucle

F0s =F0s (2:length (F0s));    % On retire aux vecteurs de sortie le premier élément qui
Vs = Vs (2: length (Vs));      % aucune signification physique

figure;                        % Affichage des résultats
subplot(3,1,1);
plot (F0s);
title ('Evolution du pitch');
subplot (3,1,2);
plot (Vs);
title ('Détermination des zones voisées(1 = voisée / 0 = non voisée)');
subplot (3,1,3);
plot (Y);
title ('Signal analisé');

    Pour ce qui est de cette partie il me semble que le programme et ses commentaires se suffisent à eux-mêmes. Seul, une description des variables d'entrées est peut-être nécessaire : les deux premières valeurs représentent le signal à analyser (Y) et le fréquence d'échantillonnage de ce dernier (Fs), le dernier élément représente le numéro de la fenêtre de 30 ms dont on souhaite voir le traitement en détails.


  3.2. SIFT

   
    Nous allons donc maintenant passer à la partie principale du programme celle qui calcule réellement le pitch de la fenêtre :     

% [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 analysée 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'entrée
E = decimate (Y,4);    % sous-échantillonnage

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

nu = 0.95;            % Préaccentuation
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 corrélation                                    %
%                                                                                         %
%*****************************************************************************************%

absE2 = abs(E2);    % Calcul du seuil et sélection des valeu
rs 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'autocorrélation

%*****************************************************************************************%
%                                                                                         %
%                 Chercher le max de l'autocorrélation entre 60 Hz et 320Hz.              %
%                                                                                         %
%*****************************************************************************************%

per1 = 4/Fs;                    % Calcul de la période correspondant à 1
'intervalle
Binf = floor (1/(400*per1));    % Calcul du nombre d'interval
les correspondant à une période de 1/400Hz
Bsup = floor (1/(60*per1));     % Calcul du nombre d'interval
les correspondant à une période de 1/60Hz

limit = len+Bsup;              % Calcul du nombre max d'interval
les à prendre en compte
Tlimit=length (C);              % en prévoyant 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 à zéro du 1er max(max +environs)afin de localisé le 2nd max
[Amax2,Imax2]= max (Cxx);       % Localisation du second max

% Vérification de l'historique des précédents 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 précision 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 décision                                       %
%                                                                                         %
%*****************************************************************************************%

if
len < 52            % Détermination 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 différentes étapes du traitement de la fenêtre              %
%                                                                                   %
%***********************************************************************************%

if aff == 1
    disp (
'Imax interpolé vaut :');
    disp(Imax);
    disp (
'Pitch');
    disp(F0);
   
    figure;
    subplot (4,2,1);
    plot (X);
    title (
'Fenêtre non traitée');
    subplot (4,2,2);
    plot (Y);
    title (
'Filtre passe-bas');
    subplot (4,2,3);
    plot (E);
    title (
'Décimation');
    subplot (4,2,4);
    plot (Y2);
    title (
'Pré-amplification');
    subplot (4,2,5);
    plot (C(len+Binf:len+Bsup));
    title (
'Autocorrélation');
    subplot (4,2,6);
    plot (rf,
'p');
    title (
'Coefficients de PARCOR');
    subplot (4,2,7);
    plot (E2);
    title (
'Résidu du LPC4');
    subplot (4,2,8);
    plot (abs(C(len+Binf:len+Bsup)),
'b');
    title (
'Décision');
    aff = 0;
    pause;
end;

   On remarque bien que le programme suit bien les étapes décrites ci- dessus, jusqu'à la partie :"Vérification de l'historique des vecteurs..." Cette partie qui peut paraître un peu obscure sert tout simplement à essayer de trouver, parmi les deux maxima calculer plus haut, lequel est le meilleur candidat à représenter le pitch réel. Pour cela, on prend en considération les deux dernières valeurs obtenues. En effet, étant donné que l'on analyse un signal de parole, il existe quelque contraintes logiques pour la valeur du pitch. Par exemple, de fortes variations de pitch d'une fenêtre à l'autre sont peu probables, nos cordes vocales ont une certaine inertie qui leurs empêche ces variations importantes. Autre exemple, étant donné que l'on travaille avec des fenêtres de 30 ms décalées à chaque fois de 10 ms, une fenêtre de 30 ms suivant deux fenêtres qui ont été déterminées comme voisées a plus de chance d'être elle aussi voisée... Nous allons reprendre cette partie du code et expliqué les significations des différents tests.

    if  (oldV(2)*oldV(1) == 0) & ((Imax/Imax2)> 1.6) & (Amax2 > (0.6*Amax))
  
         ....
   
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))
  
         ...
   
    Nous allons commencer par expliquer les termes du "elseif", nous sommes donc dans le cas où les fenêtres précédentes ont été définies comme voisées. Le premier terme vérifie que le pitch était stable à 20 % pour les fenêtres précédentes. Ensuite on vérifie que le second maximum est plus proche de la dernière valeur de pitch détectée, que ne l'est le premier. Le troisième terme sert quant à lui à vérifier que les maxima sont assez différents que pour prendre le risque de prendre en compte le deuxième maximum à la place du premier. Le dernier test vérifie que le deuxième maximum est assez important pour être représentatif. Si toutes ces conditions sont remplies, on estime alors que le deuxième maximum est plus probable que le premier et on garde donc cette valeur pour effectuer les tests de seuils finaux.
  
    Revenons au "if" de
départ, nous sommes ici dans les cas où au maximum une des deux fenêtres précédentes étaient voisées. On revient donc d'une zone non voisée, ou on arrive dans une zone voisée. Analysons le deuxième terme : celui-ci vérifie si les deux maxima sont suffisamment différents l'un de l'autres. Le troisième terme sert comme la dernière fois à vérifier si le deuxième maximum est suffisamment intense. On remarque que dans ce cas ci, les conditions d'écartement et d'intensité sont toutes les deux plus sévère que dans le premier cas.
   
    Ces tests nous permettent de garder une certaine continuité dans la détermination du pitch. Continuité qui se rallie à la continuité présente dans le signal vocal.

    La suite du programme est, elle, plus classique. La seule différence existant avec ce qui a été décrit précédemment se trouve dans les seuil de décision. On retrouve une deuxième valeur de seuil qui ne peut être utilisée qu'au cas où l'échantillon précédent était voisé. Le but de cette seconde valeur de seuil est encore une fois d'assurer une continuité des valeurs obtenues...

3.3. Résultas obtenus avec le sample "par8.wav"

    Dans cette section, je vais vous présenter les résultats obtenus avec le sample de base utilisé pour les analyses d'un signal vocal dans le cadre de notre cour. Celui-ci étant également disponible sur ce site. Cet exemple permettra en plus de se familiariser avec une des fonctions disponibles avec la toolbox Voicebox servant à extraire d'un fichier wav le vecteur des différentes valeurs de celui-ci et sa fréquence d'échantillonnage, plus encore d'autres informations; mais celles-ci ne nous servirons pas dans cet exemple.

    Première étape, chargé le fichier wav :

    Le vecteur test contient donc les différentes valeurs des échantillons du signal original, fs représente la fréquence d'échantillonnage, mode donne le nombre de bits ayant servis à la représentation du signal, format contient quant à lui des informations sur le type de codage (PCM,...), à nouveau la fréquence d'échantillonnage, nombre de bits reçu par seconde,...

    si, par exemple, on ne veut s'intéresser qu'à la partie "thé" du mot parenthèse, il nous suffit de restreindre le la zone de test entre les échantillons 2800 et 4000 pour la fin de la prononciation du t et le début du son è. Une fois le vecteur créé, il ne nous reste plus qu'à rentré le vecteur test et la fréquence d'échantillonnage de ce dernier dans la fonction SIFTECH qui va alors nous ressortir les valeurs du pitch du signal.

    Cette fonction va donc nous fournir des résultats du style :

     On retrouve donc un vecteur de sortie F0 reprenant tourtes les valeurs de pitch de sortie et le vecteur définissant les zones voisées et non-voisées. On sort de plus deux graphiques l'un reprenant les résultats généraux de l'analyse : évolution du pitch, zones voisées, et affichage du signal de base.

    Le deuxième graphique représente les différentes étapes du traitement de la fenêtre dont on voulait le détail. ces traitements débouchant sur la dernière fenêtre où l'on décide si le maximum trouvé est bon ou non. (La courbe rouge représente le seuil d'acceptation).

 

4. Conclusion

    On voit bien que , comme prévu, la méthode du SIFT donne de bons résultats. Ce qui est logique étant donné tous les avantages cités avant, d'ailleurs quand on regarde le détails de l'analyse de la fenêtre, on voit que tous ces traitements aboutissent sur un pic de l'autocorrélation bien marqué de l'autocorrélation. Ceci débouchant sur une détermination "aisée" du pitch. Le seul défaut de cette algorithme, c'est qu'il est gourmant en temps d'exécution. Ceci est principalement du au calcul de la fonction d'autocorrélation. Comme dans beaucoup de problèmes de l'ingénieur il y a un compromis à faire entre la précision du résultat et la charge de calcul que cela nécessite. Je pense que cet algorithme du SIFT en est un bon exemple.