Laboratoire 4
=============

Ce laboratoire s'intéresse à la transformée de Fourier discrète, qui
est une opération mathématique permettant de calculer les composantes
fréquentielles d'un signal à partir d'un ensemble d'échantillons de
celui-ci.  En particulier, nous allons utiliser cette opération
pour visualiser le spectre d'un signal sonore.

La transformée de Fourier discrète
----------------------------------

Considérons un signal :math:`x` dont on connaît une séquence
d'échantillons :math:`x[0]`, :math:`x[1]`, ..., :math:`x[N - 1]`,
obtenus à des instants périodiques (ce qui signifie que le délai
séparant la prise de deux échantillons successifs est constant).

Par exemple, de nombreux dispositifs de traitement du son utilisent
une fréquence d'échantillonnage de 44.1 kHz, introduite à l'origine pour
le *Compact Disc (CD)*. Le signal :math:`x` correspond dans ce cas à
l'amplitude du signal sonore échantillonné 44100 fois par seconde.

La *transformée de Fourier discrète* (*Discrete Fourier Transform, DFT*)
de la séquence :math:`x[0]`, :math:`x[1]`, ..., :math:`x[N - 1]`
est la suite de valeurs :math:`X[0]`, :math:`X[1]`, :math:`X[2]`, ...,
définie par:

      :math:`\displaystyle X[k] = \sum_{n = 0}^{N-1} x[n]\, e^{-\frac{2\pi ikn}{N}},`

pour :math:`k` entier.

Principales propriétés
~~~~~~~~~~~~~~~~~~~~~~

1) Les valeurs :math:`X[0]`, :math:`X[1]`, :math:`X[2]`, ... sont des nombres
   complexes tels que :math:`X[k + N] = X[k]` pour tout :math:`k`. (En effet,
   on a :math:`e^{-2\pi i} = 1`, donc
   :math:`e^{-\frac{2\pi ikn}{N}} = e^{-\frac{2\pi i(k + N)n}{N}}`.)

   Par conséquent, pour connaître entièrement la transformée de Fourier
   discrète des
   échantillons :math:`x[0]`, :math:`x[1]`, ..., :math:`x[N - 1]`, il
   est suffisant de se limiter à la séquence
   :math:`X[0]`, :math:`X[1]`, ..., :math:`X[N - 1]`.

2) Si les échantillons :math:`x[0]`, :math:`x[1]`, ..., :math:`x[N -
   1]` sont des nombres réels, alors on a :math:`X[k] = (X[N - k])^*`
   pour tout :math:`k`, où :math:`z^*` dénote le conjugué de
   :math:`z`. En effet, on a:

       :math:`\displaystyle X[N - k] = \sum_{n = 0}^{N-1} x[n]\, e^{-\frac{2\pi i(N - k)n}{N}} = \sum_{n = 0}^{N-1} x[n]\, e^{\frac{2\pi ikn}{N}} = (X[k])^*`.

   Cela signifie que si :math:`x[0]`, :math:`x[1]`, ..., :math:`x[N - 1]`
   sont réels, comme c'est souvent le cas lorsque l'on échantillonne un
   signal physique, alors la transformée de Fourier discrète de ces
   échantillons est entièrement caractérisée par la séquence
   :math:`X[0]`, :math:`X[1]`, ..., :math:`X[\frac{N}{2}]` (en supposant
   pour simplifier que :math:`N` est pair).
   
3) On peut retrouver :math:`x[0]`, :math:`x[1]`, ...,
   :math:`x[N - 1]` à partir de :math:`X[0]`, :math:`X[1]`, ...,
   :math:`X[N - 1]` grâce à la *transformée de Fourier discrète inverse*,
   définie de la façon suivante:

      :math:`\displaystyle x[n] = \frac{1}{N} \sum_{k = 0}^{N-1} X[k]\, e^\frac{2\pi ikn}{N},`

   pour :math:`n` entier.
  
Cette dernière propriété permet d'expliquer intuitivement la
transformée de Fourier discrète: A une constante de proportionnalité près,
les valeurs :math:`X[0]`, :math:`X[1]`,
:math:`X[2]`, ...,
:math:`X[N - 1]` correspondent aux coefficients qu'il faut appliquer aux
fonctions :math:`e^\frac{2\pi i0n}{N}`, :math:`e^\frac{2\pi i1n}{N}`,
:math:`e^\frac{2\pi i2n}{N}`, ..., :math:`e^\frac{2\pi i(N - 1)n}{N}` afin
de reconstruire les échantillons :math:`x[n]`, où :math:`n` varie de
:math:`0` à :math:`N - 1`.

Dans le :doc:`laboratoire 1</labo1>`, nous avons vu que la fonction
:math:`f(n) = e^\frac{2\pi ikn}{N}` est périodique de période
:math:`\frac{N}{k}`, pour :math:`k > 0`. (Souvenez-vous des points que
vous faisiez tourner sur le cercle trigonométrique.)
Lorsque les échantillons :math:`x[0]`, :math:`x[1]`, ...,
:math:`x[N - 1]` sont réels, on en déduit que:

- :math:`X[1]` est le coefficient d'une composante périodique de fréquence
  :math:`\frac{1}{N}`.
- :math:`X[2]` est le coefficient d'une composante périodique de fréquence
  :math:`\frac{2}{N}`.
- ...
- :math:`X[\frac{N}{2}]` est le coefficient d'une composante périodique de fréquence
  :math:`\frac{1}{2}` (en supposant cette fois encore que :math:`N` est pair).
- :math:`X[0]` est le coefficient d'une composante constante (étant donné
  que l'on a :math:`e^{2\pi i0} = 1`).

Dans ce développement, les fréquences sont exprimées par rapport à celle
utilisée pour échantillonner le signal. Par exemple, dans le cas d'un
signal sonore échantillonné à 44.1 kHz et pour :math:`N = 1024`:

- :math:`X[0]` est le coefficient d'une composante constante.
- :math:`X[1]` est le coefficient d'une composante à 43.06 Hz.
- :math:`X[2]` est le coefficient d'une composante à 86.13 Hz. 
- ...
- :math:`X[512]` est le coefficient d'une composante à 22050 Hz.

La transformée de Fourier discrète permet donc de déterminer les
composantes fréquentielles d'un signal. Dans la suite de ce laboratoire,
nous allons l'utiliser pour analyser la composition d'un signal sonore.

*Notes:*

- Chaque coefficient :math:`X[k]` est une nombre complexe. Le module
  :math:`|X[k]|` de celui-ci fournit l'*amplitude* de la composante
  périodique correspondante, et son argument :math:`\arg(X[k])` la
  *phase* de cette composante (c'est-à-dire son décalage par rapport
  au temps).

- Un échantillonnage effectué à la fréquence :math:`F` ne permet d'extraire
  les composantes du signal que jusqu'à la fréquence :math:`\frac{F}{2}`
  (appelée `fréquence de Nyquist
  <http://en.wikipedia.org/wiki/Nyquist_frequency>`_).
  
- Une bonne `référence sur la transformé de Fourier discrète
  <http://mathworld.wolfram.com/DiscreteFourierTransform.html>`_ est
  disponible sur le site de `MathWorld
  <http://mathworld.wolfram.com/>`_.

Calcul
~~~~~~

Un intérêt majeur de la transformée de Fourier discrète est que cette
opération possède un algorithme de calcul efficace, le *Fast Fourier
Transform (FFT)*.

En effet, si l'on tente de calculer :math:`X[0]`, :math:`X[1]`, ...,
:math:`X[N - 1]` à partir de :math:`x[0]`, :math:`x[1]`, ..., :math:`x[N -
1]` en appliquant naïvement la formule

      :math:`\displaystyle X[k] = \sum_{n = 0}^{N-1} x[n]\, e^{-\frac{2\pi ikn}{N}},`

alors le nombre d'opérations à effectuer est :math:`O(N^2)`, ce qui est
prohibitif pour de grandes valeurs de :math:`N`. Nous allons montrer
que ce nombre d'opérations peut être considérablement réduit en
effectuant le calcul d'une façon plus judicieuse.

Supposons, pour simplifier, que :math:`N` est une puissance entière de 2.
On peut séparer, dans l'expression de :math:`\displaystyle X[k]`, les termes
d'ordre pair et impair:

    .. math::
       X[k] &= \big(x[0] + x[2]e^{-\frac{2\pi i\,2k}{N}}
               + x[4]e^{-\frac{2\pi i\,4k}{N}} + \cdots
	       + x[N-2]e^{-\frac{2\pi i\,(N-2)k}{N}}\big)\\
	    &~~~~
	       + \big(x[1]e^{-\frac{2\pi i\,k}{N}}
               + x[3]e^{-\frac{2\pi i\,3k}{N}} + \cdots
	       + x[N-1]e^{-\frac{2\pi i\,(N-1)k}{N}}\big)\\ 
	    &= \sum_{n = 0}^{\frac{N}{2}} x[2n] e^{-\frac{4\pi ikn}{N}}
	       + \sum_{n = 0}^{\frac{N}{2}} x[2n + 1]
		 e^{-\frac{2\pi ik(2n + 1)}{N}}

Dans la deuxième somme, on peut mettre :math:`e^{-\frac{2\pi ik}{N}}`
en évidence, ce qui fournit:

    .. math::
       X[k] = \sum_{n = 0}^{\frac{N}{2}} x[2n] e^{-\frac{4\pi ikn}{N}}
               + e^{-\frac{2\pi ik}{N}}
	       \sum_{n = 0}^{\frac{N}{2}} x[2n + 1] e^{-\frac{4\pi ikn}{N}}

Cette expression montre comment décomposer le calcul de la transformée
de Fourier discrète de :math:`x[0]`, :math:`x[1]`, ..., :math:`x[N -
1]` en deux calculs plus simples. Appelons

* :math:`X_1[0]`, :math:`X_1[1]`, ..., :math:`X_1[\frac{N}{2} - 1]`
  la transformée de Fourier discrète de :math:`x[0]`, :math:`x[2]`,
  ..., :math:`x[N - 2]`.

* :math:`X_2[0]`, :math:`X_2[1]`, ..., :math:`X_2[\frac{N}{2} - 1]`
  la transformée de Fourier discrète de :math:`x[1]`, :math:`x[3]`,
  ..., :math:`x[N - 1]`.

Pour :math:`k = 0, 1, ..., \frac{N}{2} - 1`, on a donc

  .. math::
     X[k] = X_1[k] + e^{-\frac{2\pi ik}{N}} X_2[k]

Rappelez-vous que :math:`X_1` et :math:`X_2` sont périodiques
de période :math:`\frac{N}{2}` (cf. Propriété 1). On a donc également,
toujours pour :math:`k = 0, 1, ..., \frac{N}{2} - 1`:

  .. math::
      X[k + \frac{N}{2}] &= X_1[k] + e^{-\frac{2\pi i(k + \frac{N}{2})}{N}}
      X_2[k]\\
      &= X_1[k] - e^{-\frac{2\pi ik}{N}} X_2[k]

Ces expressions nous indiquent comment calculer efficacement la transformée
de Fourier discrète: Etant donnés les :math:`N` échantillons  :math:`x[0]`, :math:`x[1]`, ..., :math:`x[N - 1]` (pour rappel, nous avons supposé que :math:`N`
est une puissance de 2), il suffit de:

1) Séparer les échantillons en deux tableaux :math:`x_1` et :math:`x_2`
   de taille :math:`\frac{N}{2}`, contenant respectivement les échantillons
   d'indice pair et impair.

2) Calculer récursivement les transformées de Fourier discrètes :math:`X_1`
   et :math:`X_2` de ces deux tableaux.

   Le cas de base de cette récursion correspond à :math:`N = 1`: La
   transformée de Fourier discrète d'une séquence comprenant un seul
   échantillon :math:`x[0]` est trivialement égale à :math:`X[0] = x[0]`.

3) Combiner :math:`X_1` et :math:`X_2` en un seul tableau :math:`X` de la
   façon suivante. Pour :math:`k = 0, 1, ..., \frac{N}{2} - 1`, on a:

   .. math::
     X[k] &= X_1[k] + e^{-\frac{2\pi ik}{N}} X_2[k]\\
     X[k + \frac{N}{2}] &=  X_1[k] - e^{-\frac{2\pi ik}{N}} X_2[k]

En suivant cette procédure, chaque étape de décomposition présente
un coût total proportionnel à :math:`N`, et le nombre de décompositions
nécessaires pour arriver à des tableaux de taille 1 est égal à
:math:`\log_2 N`. Le coût de cet algorithme est donc :math:`O(N \log N)`,
ce qui est bien meilleur que l'approche naïve discutée au début
de cette section.

Implémentation
~~~~~~~~~~~~~~

Pour implémenter la FFT, vous aurez besoin de manipuler des nombres
complexes. Pour ce faire, vous pouvez utiliser la bibliothèque que
vous avez développée au :doc:`laboratoire 1</labo1>`, mais il est sans
doute plus pratique de directement exploiter les facilités offertes
par Python:

- L'instruction ``import cmath`` met à votre disposition un module
  contenant de nombreuses fonctions de manipulation de nombres complexes
  (telles que ``cmath.exp()``, ``cmath.sqrt()``, ...).

- La partie imaginaire d'une constante complexe peut être exprimée
  grâce au suffixe ``j``. Par exemple, l'instruction ``z = -1 + 2j``
  affecte à la variable ``z`` le nombre complexe :math:`-1 + 2i`.

Vous disposez à présent de tous les éléments nécessaires à l'implémentation
de la FFT. Afin de faciliter cette implémentation, il est pratique
de définir une interface qui permet de ne jamais devoir modifier le tableau
contenant les échantillons.

Votre objectif consiste à définir une fonction ``fft()`` qui:

* accepte quatre arguments:

  - une liste ``valeurs`` de nombres complexes, représentant les
    échantillons dont on souhaite calculer la transformé de Fourier
    discrète.

  - deux entiers ``debut`` et ``pas`` donnant respectivement l'indice
    du premier échantillon à considérer dans le tableau ``valeurs``,
    et la différence d'indice entre deux échantillons successifs.

    Par exemple, si ``debut = 4`` et ``pas = 2``, cela signifie qu'il
    faut traiter les échantillons ``valeurs[4]``, ``valeurs[6]``,
    ``valeurs[8]``, ...

  - un entier ``nombre`` représentant le nombre d'échantillons à
    traiter. Ce nombre doit être une puissance de 2.

* retourne une liste (de taille ``nombre``) contenant la transformée
  de Fourier discrète des échantillons. Dans le cas d'une erreur
  (par exemple, si la valeur de ``nombre`` n'est pas une
  puissance de 2), la fonction doit retourner la valeur spéciale
  ``None``.

Premiers tests
~~~~~~~~~~~~~~

Nous allons maintenant effectuer différents tests visant à vérifier
que votre implémentation se comporte correctement, et également
à vous familiariser
avec les propriétés de la transformée de Fourier discrète.

Un premier test consiste à calculer la transformée de Fourier discrète
des échantillons suivants:

  .. math::
     x[0] &= -1\\
     x[1] &= 2 + 2i\\
     x[2] &= -2i\\
     x[3] &= -3
     
Vous devriez obtenir ce résultat (aux imprécisions de calcul près):

  .. math::
     X[0] &= -2\\
     X[1] &= 1 - 3i\\
     X[2] &= -4i\\
     X[3] &= -3 + 7i

Cela fonctionne-t-il? (Dans la négative, il faut afficher les résultats
intermédiaires de vos calculs afin de déterminer les causes du problème.
Vous pouvez également écrire une fonction qui calcule de façon naïve
la transformée de Fourier discrète, afin de pouvoir en comparer les
résultats avec votre implémentation de la FFT.)

Si vous avez obtenu le bon résultat, nous allons maintenant montrer comment
calculer la transformée inverse, c'est-à-dire recalculer :math:`x[0]`,
:math:`x[1]`, :math:`x[2]` et :math:`x[3]` à partir de :math:`X[0]`,
:math:`X[1]`, :math:`X[2]` et :math:`X[3]`. Si vous remontez au début de
cette page, vous constaterez que la formule de la transformée discrète inverse
ne diffère de celle de la transformée discrète qu'en deux points:

- Les coefficients :math:`e^{-\frac{2\pi ikn}{N}}` sont remplacés par
  :math:`e^{\frac{2\pi ikn}{N}}`.
- La formule comprend un facteur supplémentaire :math:`\frac{1}{N}`.

Pour calculer la transformée inverse à l'aide du même algorithme que celui
que vous venez d'implémenter, il suffit de remarquer que pour tout
:math:`k`, on a

  :math:`e^{\frac{2\pi ikn}{N}} = e^{-\frac{2\pi i(N - k)n}{N}}`.

Calculer la transformée de Fourier discrète inverse du tableau
:math:`X[0]`, :math:`X[1]`, :math:`X[2]` et :math:`X[3]` revient
donc à:

1) Echanger les valeurs aux indices :math:`1` et :math:`N - 1`,
   celles aux indices :math:`2` et :math:`N - 2`, et ainsi de suite
   (ce qui revient à retourner le tableau entre les indices
   :math:`1` et :math:`N - 1`), la valeur d'indice :math:`0` restant
   inchangée.
    
2) Calculer la FFT du tableau obtenu.

3) Diviser chaque composante du résultat par :math:`N`.

Vous pouvez appliquer cette procédure au tableau :math:`X` précédemment
calculé. Si vous calculez la FFT de

  .. math::
     X'[0] &= -2\\
     X'[1] &= -3 + 7i\\
     X'[2] &= -4i\\
     X'[3] &= 1 - 3i\\

et divisez le résultat par 4, vous devriez retrouver les échantillons
initiaux. Est-ce bien le cas?

Fonctions périodiques
~~~~~~~~~~~~~~~~~~~~~

Le test suivant vise à mettre en évidence la capacité de la transformée
de Fourier discrète à extraire les composantes fréquentielles d'un signal.

Vous allez rédiger un programme qui effectue les opérations suivantes:

1) Générer :math:`N = 1024` échantillons d'une fonction périodique
   dont la période est un diviseur entier de :math:`N`, par exemple

     :math:`x[n] = \sin \frac{\pi n}{8}`

   pour :math:`n \in [0, 1023]`.

2) Utiliser votre fonction ``fft()`` pour calculer la transformée de Fourier
   discrète :math:`X` de cette séquence d'échantillons.

3) Ouvrir une fenêtre de 1024 pixels de large sur 768 pixels de haut.

4) A chaque position horizontale :math:`n \in [0, 1023]` dans cette fenêtre,
   afficher le module de la composante :math:`X[n]` de la transformée de
   Fourier discrète calculée au point (2), sous la forme d'une barre
   verticale dessinée à partir du bas de la fenêtre.

   (Pour dessiner cette barre verticale, le plus simple consiste à utiliser
   la fonction ``pygame.draw.rect()``, en spécifiant un rectangle de 1 pixel
   de largeur.)

Vous devriez obtenir un résultat similaire à celui-ci:

.. figure:: figures/prog-11-1-screenshot.png
  :scale: 50%
  :align: center
  :alt: DFT d'un signal périodique.

  DFT d'un signal périodique.

Comme on le voit, la figure obtenue est symétrique, puisque que les
échantillons de départ ne contiennent que des valeurs réelles. Chaque
moitié contient un pic isolé correspondant à la fréquence
:math:`\frac{64}{N}` du signal. Vous pouvez bien sûr essayer de modifier
cette fréquence afin de voir de quelle façon elle influence le résultat.

L'étape suivante consiste à superposer plusieurs signaux périodiques,
par exemple en définissant à présent

  :math:`x[n] = \frac{1}{3} \sin \frac{\pi n}{8} + \cos \frac{\pi n}{4} + \frac{1}{2} \sin \frac{\pi n}{2}`.

Comme on le voit bien sur la figure suivante, la transformée de
Fourier discrète parvient sans problème à séparer ces trois
composantes périodiques, et même à mesurer précisément leur contributions
relatives:

.. figure:: figures/prog-11-2-screenshot.png
  :scale: 50%
  :align: center
  :alt: DFT d'une superposition de signaux périodiques.

  DFT d'une superposition de signaux périodiques.

Enfin, on peut chercher à analyser un signal dont la période n'est pas
un diviseur exact de la période d'échantillonnage, par exemple

  :math:`x[n] = \sin \frac{\pi n}{5}`

Dans ce cas, on voit apparaître un phénomène de "fuite" (*spectral
leakage*) dans lequel une partie de l'énergie du pic fréquentiel
déborde sur les fréquences voisines:

.. figure:: figures/prog-11-3-screenshot.png
  :scale: 50%
  :align: center
  :alt: Fuite spectrale.

  Fuite spectrale.

(Ce phénomène de fuite se manifestera chaque fois que l'on observera
un signal dans une fenêtre dont la largeur n'est pas un multiple
exact de sa période.)

Si vous avez atteint cette étape avec un programme qui fonctionne, déposez-le
dans le répertoire centralisé, en lui donnant le suffixe ``prog-11.py``.

Spectrogramme d'un son
----------------------

Nous allons à présent utiliser la transformée de Fourier discrète afin
de visualiser le spectre d'un signal sonore.

Fichiers sonores
~~~~~~~~~~~~~~~~

La première étape consiste
à obtenir des échantillons. Nous allons utiliser le module ``wave``
de Python, permettant de travailler avec des fichiers au format
``.wav``.

Principes:

- L'instruction ``import wave`` charge ce module.

- L'instruction ``fichier = wave.open("nomfichier.wav", "r")``
  ouvre le fichier nommé ``"nomfichier.wav"`` en lecture, et
  place dans la variable ``fichier`` un objet permettant
  de le manipuler.

- Après cette instruction, la commande ``fichier.getnchannels()``
  retourne le nombre de canaux sonores du fichier: 1 pour un signal
  mono, 2 pour un signal stéréo, ...

- La fonction  ``fichier.getframerate()`` retourne la fréquence
  d'échantillonnage du fichier, exprimée en Hertz.

- La fonction  ``fichier.getsampwidth()`` retourne le nombre
  d'octets utilisés pour encoder chaque échantillon. Ce nombre est
  identique pour chaque canal.

- La fonction ``fichier.getnframes()`` retourne le nombre d'échantillons
  contenus dans le fichier. Ce nombre est identique pour chaque canal.

- La fonction ``fichier.setpos(i)`` repositionne le fichier de façon
  à ce que la prochaine lecture retourne l'échantillon ``i`` (dans
  chaque canal). Le premier échantillon correspond à ``i = 0``.

- L'instruction ``donnees = fichier.readframes(N)`` lit ``N``
  échantillons dans chaque canal, à partir de la position courante
  dans le fichier, et les place dans un tableau d'octets ``donnees``.

  Si cette opération réussit, alors le nombre d'octets contenus dans
  ``donnees`` (correspondant à ``len(donnees)``) est égal à ``N`` fois
  le nombre de canaux fois le nombre d'octets d'un échantillon. S'il
  reste moins que ``N`` échantillons dans le fichier, alors le tableau
  retourné sera plus petit.

- Dans le cas d'un signal à ``M`` canaux et à 2 octets par échantillon,
  qui est le seul que nous considérerons dans le cadre de ce laboratoire, le fragment
  de code suivant permet de transformer le résultat de
  ``donnees = fichier.readframes(N)`` en un tableau d'échantillons::

     import struct
    
     echantillons = []
     for i in range(N):
         echantillons.append(struct.unpack_from("<h", donnees, i * 2 * M)[0])
  
  *Notes:*

  - Par simplicité, on ne lit que le premier canal.
  - Chaque échantillon est représenté sous la forme d'un entier dans
    l'intervalle [-32768, 32767].
  
Dans un premier temps, afin de valider ces mécanismes, vous pouvez écrire
un programme ``prog-12.py`` qui:

1) Ouvre un fichier sonore donné.

2) En extrait les paramètres (fréquence d'échantillonnage, nombre de
   canaux, ...), et s'assure que les échantillons sont bien représentés
   sur 2 octets chacun.

3) Lit les échantillons relatifs au premier canal, et les affiche simplement
   sur la console (via ``print()``).

Pour tester votre programme, une `archive
<https://www.montefiore.ulg.ac.be/~boigelot/cours/labmp/prog/sounds.zip>`_
contenant quelques fichiers
sonores est à votre disposition.

Visualisation du spectre
~~~~~~~~~~~~~~~~~~~~~~~~

Pour visualiser le spectre du signal, nous n'allons pas calculer la
transformée de Fourier discrète de l'ensemble des échantillons. En effet,
nous souhaitons pouvoir étudier comment les composantes fréquentielles
du signal évoluent avec le temps. La procédure que nous allons suivre est
la suivante:

1) Pour une certaine valeur de :math:`N`, par exemple :math:`N = 2048`,
   extraire :math:`N` échantillons successifs du fichier, à partir d'une
   certaine position dans celui-ci.

2) Calculer la transformée de Fourier discrète de ces échantillons.

3) Afficher cette transformée dans une fenêtre, sous la forme d'une ligne
   verticale dont la position horizontale correspondra à l'endroit
   dans le fichier où les échantillons ont été extraits. En d'autres
   termes, la dimension horizontale de la figure obtenue correspondra
   au temps, et la dimension verticale à la fréquence.

Marche à suivre:

1) Dans votre programme, créer une fenêtre de 1024 pixels de large
   et de 768 pixels de haut.

2) Comme précédemment, ouvrir un fichier sonore donné, et vérifier
   que ses paramètres sont corrects. Vérifier également qu'il contient
   au moins :math:`N = 2048` échantillons.

3) Pour chaque position horizontale ``i`` dans l'intervalle
   [0, 1023]:

   a) Calculer l'endroit dans le fichier où extraire
      :math:`N` échantillons, de façon à ce que ``i = 0`` et ``i = 1023``
      correspondent respectivement au début et à la fin du fichier.

   b) Extraire ces échantillons.

   c) A l'aide de votre fonction ``fft()``, calculer la transformée
      de Fourier discrète :math:`X[0]`, :math:`X[1]`, ...,
      de ces échantillons.

   d) Afficher cette transformée dans la fenêtre, de la façon suivante:
      Pour ``j = 0, 1, ..., 767``, la couleur du pixel situé à la coordonnée
      horizontale ``i`` et la coordonnée verticale ``j`` (depuis le bas
      de la fenêtre) doit reflèter le module :math:`|X[j]|` de
      :math:`X[j]`.

Pour réaliser l'affichage, le plus simple est d'exploiter les mêmes mécanismes
que dans les programmes ``prog-3.py``, ``prog-4.py`` et ``prog-5.py``
du :doc:`laboratoire 1</labo1>`.

Un point qui mérite quelques commentaires concerne le calcul de la couleur
à attribuer au pixel représenant une composante :math:`X[j]` de la
transformée de Fourier discrète. L'approche suivante donne de bons résultats:

1) Calculer :math:`c = 1000 \log (|X[j]| + 0.1)`.
2) Pour les valeurs de :math:`c` négatives, utiliser la couleur
   (0, 0, 0) (noir).
3) Pour les valeurs de :math:`c` comprises entre 0 et 3599,
   suivre une rampe linéaire depuis la couleur (0, 0, 0) (noir) vers
   la couleur (0, 0, 255) (bleu).
4) Pour les valeurs de :math:`c` comprises entre 3600 et 7199,
   suivre une rampe linéaire depuis la couleur (0, 0, 255) (bleu) vers
   la couleur (255, 0, 255) (mauve).
5) Pour les valeurs de :math:`c` comprises entre 7200 et 10799,
   suivre une rampe linéaire depuis la couleur (255, 0, 255) (mauve)
   vers la couleur (255, 0, 0) (rouge).
6) Pour les valeurs de :math:`c` comprises entre 10800 et 14399,
   suivre une rampe linéaire depuis la couleur (255, 0, 0) (rouge)
   vers la couleur (255, 255, 0) (jaune).
7) Pour les valeurs de :math:`c` comprises entre 14400 et 17999, 
   suivre une rampe linéaire depuis la couleur (255, 255, 0) (jaune)
   vers la couleur (255, 255, 255) (blanc).
8) Pour les valeurs de :math:`c` supérieures ou égales à 18000,
   utiliser la couleur (255, 255, 255) (blanc).

Voici un exemple de ce que vous devriez obtenir, pour un des fichiers
de l'archive:

.. figure:: figures/prog-12-screenshot.png
  :scale: 50%
  :align: center
  :alt: Spectrogramme.

  Spectrogramme.
   
Si vous y êtes arrivé, félicitations! Vous pouvez placer votre programme
dans le répertoire centralisé, avec le suffixe ``prog-12.py``.

Prenez-le temps de faire quelques expériences. (L'archive contient
des exemples de sons produits par différents instruments de musique.)
Etes-vous capable de mettre en correspondance ce que vous observez sur les
spectrogrammes avec les éléments de
la musique que vous écoutez?