🚨 N'oubliez pas de mettre votre / vos nom(s) 🚨¶

TP de « Signaux aléatoires »¶

RoboVox - Algorithme d'analyse-synthèse de signaux de parole par modélisation autorégressive¶

In [ ]:
import matplotlib
import IPython.display as ipd
import numpy as np
import soundfile as sf 
import matplotlib.pyplot as plt
import scipy as sp
import librosa
import librosa.display
from utils import plot_waveform, plot_spectrum, plot_spectrogram
%matplotlib inline
%load_ext autoreload
%autoreload 2

On commence par charger un signal de parole.

In [ ]:
x_all, fs = sf.read('./data/aeiou_8k.wav')

print(fs)

plot_waveform(x_all, fs)

ipd.Audio(x_all, rate=fs)

Exercice 1 : Analyse d'une portion de signal stationnaire¶

Les voyelles dans le signal précédent durent environ 0,5 seconde et commencent à :

  • 0,25 seconde pour /a/
  • 1,2 seconde pour /e/
  • 2,2 secondes pour /i/
  • 3,1 secondes pour /o/
  • 4,1 secondes pour /u/

Dans la cellule suivante :

  • On définit une variable x qui contient une portion supposée stationnaire du signal de parole x_all, correspondant à une seule voyelle.
  • On définit une variable x_dft qui contient la transformée de Fourier discrète (TFD) de x.
  • On visualise la forme d'onde et le spectre du signal en utilisant les fonctions plot_waveform et plot_spectrum du fichier utils.py.
  • On écoute la voyelle isolée.
In [ ]:
t0vec = [int(0.25*fs), int(1.2*fs), int(2.2*fs), int(3.1*fs), int(4.1*fs)]

ind_vowel = 0

x = x_all[t0vec[ind_vowel]:t0vec[ind_vowel]+int(0.5*fs)]
T = x.shape[0]

x_dft = np.fft.fft(x)

plot_waveform(x, fs)
plot_spectrum(x_dft, fs)

ipd.Audio(x, rate=fs)

À partir du spectre, vous devriez pouvoir identifier la hauteur de la voix (la fréquence fondamentale liée à la vibration des cordes vocales, associée au premier pic du spectre) ainsi que les formants (les fréquences de résonances du conduit vocal). Vous pouvez changer la voyelle extraite de x_all, la fréquence fondamentale devrait rester la même tandis que les formants devraient changer.

Question 1 : Fonction d'autocovariance empirique¶

On note $x(t)$, $t=0,...,T-1$, les $T$ échantillons du signal observé. Compléter la fonction autocov ci-dessous pour calculer la fonction d'autocovariance empirique $\hat{R}(k)$ pour $k = 0,...,P$, définie par :

$$ \hat{R}(k) = \frac{1}{T} \sum_{t=0}^{T-1-k} x(t) x(t+k). $$

Utilisez le fait que vous traitez des numpy array, ne faites pas de boucle sur le temps !

In [ ]:
def autocov(x, P):
    """
    Calcule la fonction d'autocovariance empirique du signal.

    Entrées : 
        x : signal d'entrée, array numpy de dimension (T,)
        P : nombre de coefficients à calculer pour la fonction d'autocovariance
    Sorties : 
        r : coefficients de la fonction d'autocovariance empirique pour k=0,...,P, array numpy de dimension (P+1,)
    """
    
    r = np.zeros(P+1)

    ### A COMPLETER
    
    return r
In [ ]:
# Test
P = 16
r = autocov(x, P)
print(r)

Pour information, voici la sortie que j'obtiens pour la première voyelle ('a') :

[ 6.73034926e-03 5.46079846e-03 3.06426071e-03 8.91140319e-04 -5.41530431e-04 -7.76751448e-04 -3.51933180e-04 -2.00306322e-04 -2.27295902e-04 1.91554574e-04 8.54538583e-04 1.73829031e-03 2.57478411e-03 2.62144382e-03 1.96124498e-03 1.05283183e-03 4.20450903e-05]

Question 2 : Résolution des équations de Yule-Walker¶

Compléter la fonction yule_walkerci-dessous permettant d'estimer les paramètres du modèle $AR(P)$ à partir du signal $X(t)$ en résolvant les équations de Yule-Walker :

\begin{equation} \begin{pmatrix} \hat{R}(0) & \hat{R}(1) & \cdots & \hat{R}(P) \\ \hat{R}(1) & \hat{R}(0) & \cdots & \hat{R}(P-1) \\ \vdots & \vdots & \ddots & \vdots \\ \hat{R}(P) & \hat{R}(P-1) & \cdots & \hat{R}(0) \end{pmatrix} \begin{pmatrix} 1 \\ a_1 \\ \vdots \\ a_P \end{pmatrix}= \begin{pmatrix} \sigma^2 \\ 0 \\ \vdots \\ 0 \end{pmatrix}. \end{equation}

On pourra utiliser les fonctions suivantes :

  • sp.linalg.toeplitz pour former une matrice Toeplitz ;
  • np.linalg.inv pour calculer l'inverse d'une matrice ;
  • @ pour le produit matriciel sur des arrays Numpy ;
  • np.concatenate pour concatener deux arrays Numpy.
In [ ]:
def yule_walker(x, P):

    a = np.zeros(P)
    sigma2 = 1
    
    # A COMPLETER

    return a, sigma2
In [ ]:
# TEST
P = 16
a, sigma2 = yule_walker(x, P)
print(a)
print(sigma2)

Pour information, voici la sortie que j'obtiens pour la première voyelle ('a') :

[-1.44087875 0.85827083 -0.38642303 0.60496031 -0.60492978 0.13583144 -0.12726464 0.61183839 -0.72054679 0.38558782 -0.26901521 0.1381512 -0.06034668 0.01878405 -0.10915575 0.11817264] 0.0009767540907426792

Question 3 : Calcul du résiduel¶

  • a) Complétez la cellule suivante pour calculer le résiduel $e(t) = x(t) - \hat{x}(t)$ avec $\hat{x}(t) = - \sum_{i=1}^P \hat{a}_i x(t-i)$ (on considérera $x(t)=0$ pour $t < 0$).

  • b) En quoi le spectre du résiduel est-il différent de celui du signal original ?

  • c) Écoutez les résiduels obtenus à partir de différentes voyelles, qu'observez-vous ?

In [ ]:
def residuel(x, a):

    T = x.shape[0]
    e = np.zeros(T)
    x_pred = np.zeros(P)

    for t in np.arange(T):
        # A COMPLETER
        pass

    return e

res = residuel(x, a)

res_dft = np.fft.fft(e)

plot_waveform(e, fs)
plot_spectrum(e_dft, fs)
ipd.Audio(e, rate=fs)

Réponse à la question 3¶

Ecrivez ici votre réponse.

Question 4 : Calcul de l'enveloppe spectrale¶

  • a) Calculez et tracez le périodogramme du signal défini par :

    $$ \frac{1}{T} \Big| \sum_{t=0}^{T-1} x(t) e^{-\imath 2 \pi \nu t} \Big|^2 $$

    En pratique, vous calculerez la transformée de Fourier discrète (TFD) en utilisant la fonction np.fft.fft. La TFD d'ordre $F = T$ correspond à une discrétisation de la transformée de Fourier à temps discret (TFTD), où la fréquence réduite $\nu$ est discrétisée de la façon suivante : $\nu_f = f/F$ avec $f \in \{0,...,F-1\}$. Le signal étant réel, sa TF(T)D est à symétrie hermitienne, vous ne tracerez donc que la partie non redondante du spectre contenue dans les $T//2 + 1$ premiers coefficients.

  • b) Calculez et tracez sur la même figure l'enveloppe spectrale normalisée définie par :

    $$ \frac{\hat{\sigma} ^2}{ \Big| 1 + \sum_{j=1}^P \hat{a}_j e^{-\imath 2 \pi \nu j} \Big|^2}, $$

    où $\{\hat{a}_j\}_{j=1}^P$ et $\hat{\sigma} ^2$ sont les estimations des paramètres du modèle AR obtenues précédemment. Notez que le dénominateur de l'enveloppe spectrale correspond à la TFTD de la séquence $\{1, a_1, ..., a_P\}$. En pratique, comme précédemment, on calculera la TFD à l'aide de la fonction np.fft.fft.

  • c) Qu'observez-vous ?

In [ ]:
perio_x_db = np.zeros(T//2+1)
spec_env_db = np.zeros(T//2+1)

# A COMPLETER

freq  = np.arange(0, T//2+1)*fs/T

plt.figure(figsize=(12,5))
plt.plot(freq, perio_x_db)
plt.plot(freq, spec_env_db)
plt.xlabel('fréquence (Hz)')
plt.ylabel('puissance (dB)')
plt.show()

Réponse à la question 4¶

Ecrivez ici votre réponse.

Question 5 : Synthèse¶

  • a) Complétez la cellule suivante pour générer récursivement le signal de parole (on considérera $x(t)=0$ pour $t < 0$) :

    $$ x(t) = \hat{x}(t) + w(t), \qquad \hat{x}(t) = - \sum_{i=1}^P \hat{a}_i \, x(t-i), $$

    où la source $w(t)$ est soit :

    • le signal résiduel que vous avez calculé précédemment ;
    • un bruit blanc gaussien de variance $\sigma^2$ ;
    • un train d'impulsions périodiques d'amplitude $\sigma$.

  • b) Exécutez la cellule encore après afin de visualiser la forme d'onde, le spectre, et écouter le signal synthétisé. Décrire ce que vous observez / entendez pour les différents types de source.

In [ ]:
source = 'non-voisee' # au choix : 'residuel', 'non-voisee', or "voisee"
f0 = 120
T0 = int(1/f0*fs)

if source=='voisee':
    w = np.zeros(T)
    w[0:T:T0] = np.sqrt(sigma2*T0)
elif source=='non-voisee':
    w = np.sqrt(sigma2)*np.random.randn(T)
elif source=='residuel':
    w = res  
    
def synthese(w, a):
    
    T = w.shape[0]
    x_gen = np.zeros(T)

    # A COMPLETER

    return x_gen

x_gen = synthese(w, a)

x_gen_dft = np.fft.fft(x_gen)

plot_waveform(x_gen, fs, title='Forme d\'onde du signal synthétisé')

plot_spectrum(x_gen_dft, fs, title='Spectre du signal synthétisé')

ipd.Audio(x_gen, rate=fs)

Réponse à la question 5¶

Ecrivez ici votre réponse.

Exercice 2 : Analyse à court terme d'un signal non stationnaire¶

Dans les expériences précédentes, nous avons analysé et synthétisé un court signal de parole stationnaire. Cependant, en pratique, nous devons traiter des signaux de parole non stationnaires de longueur arbitraire, comme le signal chargé dans la cellule suivante. Vous trouverez dans le dossier data d'autres signaux de parole à explorer. N'hésitez pas à vous enregistrer pour travailler sur votre propre voix !

In [ ]:
fs = 16000
x_all, fs_x = sf.read('./data/speakerM1_uttNum1459.wav')

if fs_x != fs:
    # resample to 8 kHz, if necessary
    ratio = float(fs) / float(fs_x)
    n_samples = int(np.ceil(x_all.shape[-1] * ratio))
    x_all = sp.signal.resample(x_all, n_samples, axis=-1)
    
x_all = x_all - np.mean(x_all)

ipd.Audio(x_all, rate=fs)

Nous allons utiliser le modèle source-filtre pour analyser, transformer, et resynthétiser ce signal de parole. Ce processus d'analyse/synthèse doit être effectué sur de courtes trames recouvrantes du signal, où l'on peut le supposer stationnaire afin d'être cohérent avec le modèle $AR(P)$. Cette procédure de fenêtrage est la même que celle utilisée dans la transformée de Fourier à court terme.

Le processus global est le suivant, pour toutes les trames temporelles :

  1. Extraire une courte trame du signal de parole et la multiplier par une fenêtre d'analyse.

  2. Analyse / Encodage :

    a. Calculer les paramètres du modèle source-filtre, c'est-à-dire les coefficients du filtre dans le modèle AR.
    b. Calculer le signal résiduel et son énergie.
    c. Détecter si le signal résiduel est voisé ou non voisé. S'il est voisé, calculer la fréquence fondamentale.

  3. Synthèse / Décodage :

    a. Synthétiser le signal source et normaliser son énergie pour qu'elle corresponde à celle du résidu.

    • Si le signal dans la fenêtre d'analyse courante est non voisé, la source est modélisée par un bruit blanc gaussien.

    • Si le signal est voisé, la source est modélisée par un train d'impulsions périodiques dont la fréquence fondamentale correspond à celle du résidu.

    b. Synthétiser un signal de parole suivant le modèle source-filtre, en utilisant le signal source synthétique et les coefficients du filtre dans le modèle AR.

  4. Reconstruire le signal complet par addition-recouvrement (overlap-add).

Vous devez compléter la cellule suivante pour implémenter ce processus ; vous serez guidé par les commentaires dans le code.

Remarque : L'étape 2.c n'est pas implémentée dans le code ci-dessous. Pour simplifier la procédure, nous synthétisons une source qui est voisée ou non voisée pendant toute la durée du signal. Ce signal source est contrôlé par les variables voiced et f0. Si on fixe une source voisée de $f0$ constant, on obtient une voix robotique, car dépourvue de variations de hauteur. Cependant, lorsqu'un signal est périodique, sa fonction d'autocovariance présente des maxima séparés par la période de la fondamentale. La fréquence fondamentale pourrait donc être mesurée en trouvant le maximum d'une estimation de la fonction d'autocovariance normalisée du résidu. Cette fonction pourrait également être utilisée pour détecter l'activité de voisement (par exemple, le signal est considéré comme voisé si le maximum de la fonction est supérieur à 0,5). Cette partie est laissée comme exercice bonus.

In [ ]:
L = int(0.030*fs) # longueur de la fenêtre d'analyse
H = L//2 # pas d'analyse
P = 16 # ordre du modèle AR

voiced = True # booléen indiquant si la voix est voisée ou non
f0 = 120 # fréquence fondamentale en Hz (si voisée)
T0 = int(1/f0*fs) # prédiode associée

T = x_all.shape[0]

N = int(np.fix( (T-L)/H)) # nombre de trames

win = np.sin(np.arange(.5,L-.5+1)/L*np.pi); # fenêtre d'analyse sinusoidale

x_gen_all = np.zeros(T) # signal synthétisé
e_all = np.zeros(T) # résiduel

# Boucle sur les trames
for n in np.arange(N):

    # sélectionner une petite portion de signal et la multiplier par une fenêtre d'analyse
    n1 = ?
    n2 = ?
    x = x_all[n1:n2]*win
    
    # calculer les coefficients du filtre du modèle AR
    a = ?
    
    # calculer le résiduel et son énergie moyenne (moyenne du carré des coefficients)
    e = ?
    residual_energy = ?

    # synthétiser le signal source
    if voiced:
        w = ?
    else:
        w = ?

    # ajuster son énergie
    w = w / np.sqrt(np.mean(w**2)) * np.sqrt(residual_energy)

    # synthétiser le signal de parole
    x_gen = ?

    # effectuer l'addition-recouvrement
    x_gen_all[n1:n2] += x_gen*win

    # on stocke le résidu, au cas où
    e_all[n1:n2] += e*win
In [ ]:
plot_waveform(x_gen_all, fs, title="Forme d\'onde du signal transformé")

ipd.Audio(x_gen_all, rate=fs)

Exercice bonus : Synthèse croisée¶

La modélisation AR de la parole permet des applications au-delà du codage. Un bon exemple est celui de la synthèse croisée.

Après avoir résolu les équations de Yule-Walker, nous obtenons une paramétrisation du filtre (les coefficients $a_i$) et une estimation de la source (le résiduel). Si nous appliquons cette méthode sur deux signaux différents, nous pouvons utiliser la source d'un signal et le filtre de l'autre pour créer un signal hybride. C'est ce qu'on appelle la synthèse croisée, qui donne des résultats particulièrement intéressants lorsqu'on combine le filtre provenant d'un signal de parole avec une source sonore spectralement riche, comme une porte qui craque, des bulles qui éclatent, ou une voiture de course qui accélère.

Exercice bonus : utilisez les signaux sources chargés dans la cellule suivante pour effectuer une synthèse croisée avec le signal de parole précédent.

In [ ]:
fs = 16000

# w_all, fs_w = sf.read('./data/creak.wav')
w_all, fs_w = sf.read('./data/racing.wav')
# w_all, fs_w = sf.read('./data/bubbles.wav')

if w_all.ndim > 1:
    w_all = w_all[:,0]

if fs_w != fs:
    # resample to 8 kHz, if necessary
    ratio = float(fs) / float(fs_w)
    n_samples = int(np.ceil(w_all.shape[-1] * ratio))
    w_all = sp.signal.resample(w_all, n_samples, axis=-1)
    
w_all = w_all - np.mean(w_all)

w_all = w_all[:x_all.shape[0]]

ipd.Audio(w_all, rate=fs)
In [ ]:

Références¶

  • Linear Predictive Coding is All-Pole Resonance Modeling, article by Hyung-Suk Kim, Center for Computer Research in Music and Acoustics, Stanford University.
  • Time-varying linear predictive coding of speech signals, Doctoral dissertation at the MIT by Mark Gilbert Hall (1977).