TP2 : systèmes et SLITs

Dans ce deuxième TP, nous revisitons quelques notions associées aux systèmes, et en particulier aux SLITs.

Note habituelle : Le cours porte uniquement sur les signaux analogiques, c'est-à-dire définis sur un espace de temps continu (et à valeurs dans un espace continu). Nos ordinateurs "parlent" de manière discrète, et il est donc difficile de manier des signaux analogiques. Ici, nous nous contenterons d'échantillonner très finement nos signaux, et considèrerons qu'il s'agit bien de signaux analogiques. Cependant, nous serons dès lors en train de manipuler des signaux échantillonnés. Il existe en réalité une théorie et un cadre mathématique pour traiter du passage des signaux analogiques aux signaux échantillonnés, que beaucoup d'entre vous rencontreront probablement plus tard dans votre cursus universitaire (ou ont peut-être déjà rencontré).

Comme précédemment, commençons par quelques imports qui nous seront utiles. Cette fois-ci, en plus des librairies extérieures usuelles, nous importons aussi utils, qui correspond à un module (i.e. un fichier .py) que vous avez reçu en même temps que ce notebook. Il contient certains des outils et fonctions du TP1, pour vous éviter d'avoir à les copier/coller ici, ainsi que quelques fonctions dont nous nous servirons pour effectuer une vraie application de vos connaissances dans la dernière section de ce TP.

%matplotlib inline

import numpy as np
from matplotlib import pyplot as plt
import scipy.io.wavfile as wavfile

import utils

Lors du TP1, nous avons étudié nos signaux sur la plage $t\in[-5, 5]$. Dans celui-ci, nous utiliserons $t\in[-7.5, 7.5]$ :

t_echant = np.linspace(-7.5, 7.5, 10001)

1. Représentation d'un système en python

Commençons par représenter le système suivant :

A vous : Ecrire une fonction system_1 qui prend en entrée un signal x et lui applique le système ci-dessus. Codez explicitement chacune des opérations représentées dans les blocs de la figure ci-dessus.

Note : Vous pouvez, selon vos préférences, vous inspirer de la fonction affine du TP1 et faire que system_1 ne prenne pas t comme argument et retourne une fonction python (qui, elle-même, prend t comme argument) ; ou, plus simplement, faire que system_1 recoive un signal x et un temps t comme arguments et retourne uniquement la valeur prise par le signal transformé.

def system_1(x):
    #TODO
    return

Considérons maintenant le système suivant : $$x(t) \rightarrow y_2(t) = \left|x(t) - x(t-t_0)\right|$$

A vous : Ecrire une fonction system_2 qui correspond à ce système.

def system_2(x):
    #TODO
    return

A vous : Appliquez les deux systèmes au signal d'entrèe de votre choix (si vous êtes en manque d'inspiration, vous trouverez les signaux $x_1$ et $x_2$ du TP1 dans utils). Représentez les deux signaux transformés sur la même figure, sur la plage $t\in[-7.5, 7.5]$.

#TODO

plt.plot()
plt.plot()
plt.show()
plt.close()

Question : Que constatez-vous ?

[votre réponse]

A vous : Testez empiriquement votre affirmation sur la plage $t\in[-7.5,7.5]$.

#TODO

2. Produit de convolution

Le produit de convolution, noté $x(t) * h(t)$, est un opérateur très important en traitement du signal. Comme vous l'avez vu, il permet notamment de caractériser uniquement un SLIT.

En python, l'une des possibilités pour effectuer des convolutions (discrètes) entre deux signaux uni-dimensionnels est d'utiliser la fonction convolve de numpy. Jupyter nous permet d'accéder directemment à la documentation d'objets python avec la syntaxe suivante :

?np.convolve

Aujourd'hui, nous utiliserons toujours le mode "same" afin de conserver la dimension de nos signaux.

A vous : Effectuez la convolution entre deux portes rectangulaires $\Pi_T(t)$ et la représenter.

#TODO

plt.plot()
plt.plot()
plt.show()
plt.close()

Dans le cas continu, appliquer une convolution par une impulsion de Dirac, $\delta(t)$, revient à appliquer l'identité, c'est-à-dire que pour tout $x$, on a

$$x(t)*\delta(t) = x(t)$$.

L'impulsion de Dirac continue, telle que vue en cours, est difficile à manier en python. A défaut, étudions l'effet de la convolution par un signal qui vaut $0$ partout et $1$ à $t=0$ :

def dir_dis(t):
    if t==0:
        return 1
    return 0

A vous : Représentez ce signal.

#TODO
plt.plot()
plt.show()
plt.close()

Bonus : Vous avez peut-être constaté que dans la deuxième cellule de ce notebook, nous avons défini t_echant en utilisant 10001 échantillons. Que se serait-il passé ici si nous en avions tiré 10000, comme lors du TP1 ?

Indices pour le bonus (cliquer pour afficher)

En cas de doute ou de confusion, vous pouvez essayer de représenter le signal ci-dessus sur la même plage, mais en n'utilisant que 10000 échantillons (avec un nouveau t_echant). Si le résultat vous surprend, essayez d'étudier de plus près ce qu'il s'est passé (par exemple en inspectant les valeurs de votre nouveau t_echant).

A vous : Etudiez l'effet de la convolution du signal de votre choix par dir_dis en représentant le signal choisi et le résultat de la convolution sur la même figure.

#TODO 

#TODO
plt.plot()
plt.show()
plt.close()

Question : que constatez-vous ?

[votre réponse ici]

3. Filtre moyenneur

Dans cette section, nous parlons d'un type de système bien particulier et très courant : les filtres. Aujourd'hui, revisitons le filtre moyenneur que vous avez rencontré en TD. Rappelons qu'il peut être défini par sa réponse impulsionnelle : $$h(t) = \Pi_T(t-T/2)$$.

A vous : définir une fonction moyenneur qui applique le filtre moyenneur à un signal d'entrée x en utilisant sa réponse impulsionnelle.

def moyenneur(x, T):
    #TODO
    return

A vous : L'appliquer aux signaux d'entrée $\Gamma(t)$ et $x_1(t) = \mathrm{cos}(2\pi f_0t)$, puis en représenter les sorties. Vérifier qu'elles correspondent à celle de votre TD.

#TODO

4. Une vraie application !

Dans cette dernière section, nous nous intéressons à la transformée de Fourier, et en particulier à son application pour appliquer un filtrage à un signal d'entrée. Cette section est inspirée des TPs de traitement du signal donnés ici même, ainsi que ceux donnés à l'Ecole Polytechnique, par Rémi Flamary.

Le fichier batterie.wav contient un court extrait d'enregistrement de batterie contenant des symballes et de la grosse caisse. L'un correspond à un signal de haute fréquence, l'autre de basse fréquence. Nous allons nous servir des connaissances que vous avez acquéries dans ce cours pour tenter de séparer les deux pistes du fichier.

Le signal d'entrée est bien sûr à valeurs réelles, mais comme nous allons travailler dans l'espace de Fourier, nous allons devoir manipuler des nombres complexes. Commençons donc par définir une petite fonction qui nous permet simplement de représenter une fonction à valeurs complexes :

def plot_complex(x, y, alpha=.6, **kwargs):
    """Represente une fonction à valeurs complexes.
    
    Parametres
    ----------
    x : iterable
        Les valeurs d'échantillonnage de la fonction (e.g. temps ou frequence).
        Correspond à l'axe des abscisses.
    y : iterable
        Valeurs complexes prises par la fonction à chaque point donné en ``x``.
    alpha : float, optional
        Pourcentage d'opacité de la courbe des valeurs imaginaires.
    Tous les autres arguments seront passés à `plt.plot`.
    """
    plt.plot(x, np.real(y), label="Partie Réelle", **kwargs)
    plt.plot(x, np.imag(y), label="Partie Imaginaire", alpha=alpha, **kwargs)
    plt.legend()
    plt.show()
    plt.close()

Après avoir écouté le fichier batterie.wav, importons-le en python. Notez que nous avons importé scipy.io.wavfile (sous l'alias wavfile) en début de notebook. Il s'agit d'un sous-module de scipy qui fournit des fonctions pour manipuler les fichiers .wav. Il contient en particulier deux fonctions, read et write qui permettent de lire des fichiers .wav et de les convertir en arrays numpy, et vice-versa. Commençons par regarder la syntaxe de la première :

?wavfile.read

On voit qu'elle retourne deux quantités ; le taux d'échantillonnage du signal, et le signal lui-même.

A vous : Utilisez la fonction wavfile.read pour lire notre fichier batterie (stockez le taux d'échantillonnage dans une variable samp_rate et le signal dans une variable batterie):

#TODO: samp_rate, batterie = ...

Inspectons un peu l'allule de batterie :

print(batterie)
print(batterie.dtype)
print(batterie.shape)
n_samp = len(batterie)

A vous : Représentez le signal.

plt.plot() #TODO
plt.show()
plt.close()

Appliquons maintenant la transformée de Fourier à notre signal. Comme toujours, comme nous travaillons ici avec des signaux discrets, nous ne pouvons pas simplement appliquer la définition, dans le cas continu, que vous avez vu en cours. Nous allons en fait appliquer un algorithme, la Fast Fourier Transform (FFT), adaptée aux signaux échantillonnés. Vous en apprendrez certainement plus sur la FFT plus tard ; pour ce TP, contentons-nous de considérer qu'il s'agit d'un moyen d'obtenir la transformée de Fourier d'un signal.

Notre module fait-maison utils contient deux fonctions, fourier et inverse_fourier, qui effectuent la transformée de Fourier, et la transformée de Fourier inverse, respectivement, pour des signaux semblables à celui que nous allons manipuler aujourd'hui. Notez que cette dernière ne conserve que la partie réelle du signal !

Appliquons la fonction utils.fourier et inspectons le résultat :

Fbatterie = utils.fourier(batterie)
print(Fbatterie.shape, Fbatterie.dtype)
print(Fbatterie)

Déterminons les fréquences correspondant à chaque échantillon de notre transformée de Fourier :

freqs = np.fft.fftshift(np.fft.fftfreq(n_samp, d=1/samp_rate))

A vous : Utilisez la fonction plot_complex d'utils pour représenter la transformée de Fourier de notre signal;

#TODO

Passons maintenant au filtrage ! Vous avez vu et étudié en cours les filtres passe-bas du premier ordre, caractérisés par la fonction de transfert $$H(\omega) = \frac{1}{1 + j\frac{\omega}{\omega_c}}$$.

Aujourd'hui, nous allons nous contenter du filtre passe-bas idéal, défini plus simplement par sa fonction de transfert :

$$H(\omega) = \begin{cases}1&\mathrm{si}&|\omega| < \omega_c,\\ 0&\mathrm{sinon.}\end{cases}$$

A vous : Ecrivez une fonction passe_bas appliquant la fonction de transfert ci-dessus.

Bonus : Ecrivez également une fonction appliquant la fonction de transfert d'un filtre passe-bas du premier ordre. Représentez-la.

def passe_bas(w, w_c=500):
    #TODO
    return 

A vous : Appliquez la fonction de transfert à toutes les fréquences de freqs. Représentez $H(\omega)$ sur cette plage de fréquence.

Hw = ... #TODO
#TODO

Appliquons maintenant le filtre à notre signal. Rappelez-vous que dans le domaine de Fourier, l'application d'un SLIT de fonction de transfert $H$ revient à la simple multiplication :

$$y(t) = \mathcal{F}^{-1}[Y(\omega)] = \mathcal{F}^{-1}[X(\omega)H(\omega)] = x(t)*h(t).$$

A vous : Appliquez le filtre passe-bas dans le domaine de Fourier pour obtenir $Y(\omega)$.

filt_Fbatterie = #TODO

A vous : Représentez, toujours dans le domaine de Fourier, le signal filtré.

#TODO

Appliquons maintenant la transformée de Fourier inverse :

filt_batterie = utils.inverse_fourier(filt_Fbatterie)

A vous : Représentez, sur le même graphe, notre signal avant et après filtrage.

Indice : Pour bien voir les deux signaux même là où ils se superposent, vous pouvez vous inspirer du code de plot_complex.

plt.plot()
plt.plot()
plt.legend()
plt.show()
plt.close()

Convertissons à présent notre signal filtré en fichier .wav afin de pouvoir l'écouter :

A vous : En utilisant wavfile.write, écrivez un fichier .wav contenant le signal filtré. Notez que le taux d'échantillonnage est le même que celui que nous avons récupéré pour batterie.

#TODO

A vous : Ecoutez le fichier ainsi produit. Que constatez-vous ?!

[votre réponse]

Nous venons en fait d'effectuer ce qu'on appelle la séparation de sources. Comme vous vous en doutez, c'est une application très utile du traitement du signal, par exemple (mais pas seulement !) dans le domaine de la musique.

Pour finir, puisque nous avons isolé l'une de nos deux sources, essayons à présent de récupérer une piste ne contenant que la seconde.

A vous : Retirez le signal filtré au signal original.

#TODO

A vous : Représentez le signal ainsi obtenu (éventuellement sur la même figure que le signal complet et celui filtré).

plt.plot()
plt.plot()
plt.plot()
plt.legend()
plt.show()
plt.close()

A vous : Ecrivez ce nouveau signal dans un fichier .wav et écoutez-le. Que constatez-vous ?

#TODO

[votre réponse]

Bonus : Repétez l'opération en changeant la fréquence de coupure $\omega_c$ pour essayer d'isoler les cymballes. Que constatez-vous ?

[votre réponse]