Algorithme FFT - Fast Fourier Transform

Au début du XIXe siècle, dans ses travaux sur la chaleur, Joseph Fourier propose une méthode pour transformer une fonction périodique en une somme de fonctions trigonométriques. Ce principe sera ensuite abondamment développé et se révélera d'une importance capitale en particulier en physique et en traitement du signal.

C'est mathématiquement trop compliqué pour des lycéens, donc nous allons essayer de comprendre les idées générales sans trop entrer dans les détails.

Prenons la corde d'un instrument : violon, guitare… La corde est bloquée de deux côtés et on la fait vibrer. La corde a plusieurs façons possibles de vibrer. On parle de modes de vibrations.

Vous voyez que certains modes ont plus d'ondulations. Une onde plus courte correspond à une fréquence plus élevée. Donc on peut dire que la corde peut onduler selon plusieurs fréquences.

La matière et la tension de la corde sont réglées de façon que F corresponde à une certaine note. Supposons que F corresponde au La du diapason, c'est à dire $F = 440\,Hz$.

Hz = Hertz = unité de fréquence. 440 Hz signifie que l'on 440 vibrations par seconde.

Le second mode aura donc une fréquence de $2F = 880\,Hz$, le 3e de $3F = 1320\,Hz$…

La vibration de la corde entraîne des variations de pression dans l'air. Ces variations de pression se propagent, c'est le son. Elles arrivent dans notre oreille. La vibration d'air fait vibrer le tympan… et au bout la cochlée réalise un traitement qui s'apparente à la transformation de Fourier qui nous intéresse.

On peut aussi raisonner d'un point de vue électrique : le son arrive au microphone, fait vibrer la partie sensible du micro qui transforme ces variations de pression en signal électrique : une tension qui varie en fonction du temps. On peut ensuite, électroniquement, faire subir des traitements à ce signal : filtrage, amplification, etc.

Si la corde ne vibrait que selon le mode du haut à 440 Hz, on aurait alors une fréquence pure. Le signal serait un pur sinus et le son ressemblerait à celui ci.

Mais la corde mélange tous les modes possibles et la vibrations est la somme de plusieurs modes, plus ou moins forts. Par exemple le mode F aura une amplitude de 1, 2F une amplitude de 0.3, 3F un amplitude de 0.5 et 4F une amplitude de 0.2. Voici le signal résultant de la somme de ces 3 composantes :

Dans l'exemple précédent, le signal est compliqué, bien plus qu'un simple sinus. Nous ne voyons pas à l’œil qu'elle n'est que la somme de 4 sinus :

• un sinus d'amplitude 1 a la fréquence F
• un sinus d'amplitude 0.3 à la fréquence 2F
• un sinus d'amplitude 0.5 à la fréquence 3F
• un sinus d'amplitude 0.2 à la fréquence 4F

Et pourtant notre oreille le perçoit. En effet, dans la cochlée, le son se propage et agite des cils. Les cils mis en mouvements correspondent justement aux fréquences du signal. Ainsi notre cerveau reçoit que le son correspondant est composé de 4 certaines fréquences avec une certaines intensité.

Nous y sommes presque : La transformée de Fourier va permettre, partant du signal sonore compliqué, de retrouver les composantes du signal. La transformée de Fourier va nous donner le spectre, c'est à dire la décomposition :

$$F = 440\,Hz\quad F\to 1 \quad 2F\to 0,3 \quad 3F \to 0,5 \quad 4F \to 0,2$$

Ce signal est plus compliqué qu'un simple sinus. On peut appeler fondamental le mode de la première fréquence et harmoniques les autres. En musique, le fondamental a une importance particulière.

Dans le cadre musical, c'est la géométrie, la matière de l'instrument qui décide de son spectre. Le spectre est un élément important de ce qui fait qu'un son de violon n'est pas un son de piano et réciproquement. Cependant, d'autres éléments sont plus importants encore.

Dans le cadre d'un fichier numérique, une complication s'ajoute. Les signaux sont échantillonnés. C'est à dire que l'on ne travaille pas sur un signal continu comme sur les courbes précédentes. On ne dispose que de points prélevés à intervalle fixe.

Sur le graphique ci-dessus, on a relevé les points de la courbe à intervalle $T_e$, e pour échantillonnage. On utilise aussi souvent la fréquence d'échantillonnage $F_e = \frac{1}{T_e}$ qui s'exprime en Hz.

On se doute que l'on doit prendre assez de points, c'est à dire que $T_e$ doit être assez petit, ou $F_e$ assez grand. Le théorème de Shanon nous permet de chiffrer ce “assez”. On doit choisir $F_e$ au moins deux fois plus grand que la fréquence dans le spectre du signal. Par exemple, on sait que l'humain est capable d'entendre des sons jusqu'environ $20\,kHz$ – au-delà, on parle d'ultrasons – et donc on choisira $F_e \geqslant 40\,kHz$. Les CD audio sont à $44,1\,kHz$.

On fournit à l'algorithme fft un tableau de flottants signal.

• signal représente le signal échantillonné.
• La fréquence d'échantillonnage Fe est connue. On sait donc que la période d'échantillonnage est Te = 1/Fe.
• signal compte n = len(signal) items. On en déduit donc que cela représente une durée T = n * Te.

Il n'est pas nécessaire que fft connaisse T, Te et Fe. Ces grandeurs ne sont utiles que pour interpréter le résultat. fft renvoie un tableau que nous appellerons spectre.

• spectre est un tableau de nombres complexes.
• Nous pourrons adapter ce résultat pour ne conserver que le module de ces nombres.

>>> complex(3,2)
3+2j
>>> z = complex(3,2)
>>> abs(z)
3.6055512754639896

• Le tableau contient n items, autant que signal.
• L'item de rang i correspondant à la fréquence i * Fe/n = i * 1/T.

• Pour une raison liée au théorème de Shanon, le tableau spectre présente une symétrie de sorte qu'à l'indice n-i on a le complexe conjugué de l'indice i.
• Pour i != 0, si on veut interpréter spectre[i] en termes de sinus ou cosinus, il faudra multiplier par 2.

La théorie nous dit que $\displaystyle raie[\ell] = \frac{1}{n} \left|\sum_{k = 0}^{n-1} signal[k] \cdot \exp\left(\frac{-2\,i\,\pi}{n}k\ell\right)\right|$ ou raie est l'amplitude d'une raie du spectre.

En Python, on pourrait poser :

from cmath import exp, pi
def raie_spectre(signal, l):
    n = len(signal)
    I = complex(0, 1) # le nombre spécial tel que I**2 = -1
    acc = 0
    for k, s in enumerate(signal):
        acc += s * exp(-2*pi*I/n * k * l)
    return abs(acc)/n

Mais c'est algorithme est beaucoup trop lent, d'où l'intérêt d'une version rapide – D'où le qualificatif de fast dans fast fourier transform.

On veut calculer $\displaystyle X_\ell = \sum_{k = 0}^{n-1} x_k \exp\left(-\frac{2i\,\pi}{n} k\ell\right)$.

On commence par poser $w(n) = \exp\left(-\frac{2i\,\pi}{n}\right)$ de sorte que le calcul devient $\displaystyle X_\ell = \sum_{k = 0}^{n-1} x_k w(n)^{k\ell}$, avec de plus $w(n)^n = 1$.

$X_k$ servirons à calculer les raies du spectre et les $x_k$ sont les éléments de signal.

Le calcul ci-dessus peut-être décrit par cette figure.

La liste des $x_k$ permet de calculer la liste des $X_k$, en tenant compte du facteur $w(n)$.

L'astuce du calcul utilisée par l'algo FFT consiste à calculer la FFT pour les $\frac{n}{2}$ valeurs d'indice paire, la FFT pour les $\frac{n}{2}$ valeurs d'indice impair, et de combiner les 2 résultats.

Cela suppose que $n$ est pair. Et pour itérer le principe, on demandera même que $n$ soit une puissance de 2.

Ce que signifie ce schéma :

• $FFT(n)$ consiste à calculer les $n$ valeurs $X_k$ pour le signal constitué des $n$ termes $x_k$.
• Pour cela on fait le calcul $FFT\left(\frac{n}{2}\right)$ sur les $x_k$ avec $k$ pair. Cela nous donne $\frac{n}{2}$ valeurs $A_k$.
• On fait aussi le calcul $FFT\left(\frac{n}{2}\right)$ sur les $x_k$ avec $k$ impair. Cela nous donne $frac{n}{2}$ valeurs $B_k$.
• La première moitié des $X_k$ est obtenue par $X_k = A_k + B_k \cdot w(n)^k$
• La deuxième moitié des $X_k$ est obtenue par $X_{k} = A_{k-\frac{n}{2}} + B_{k-\frac{n}{2}}\cdot w(n)^k$

Bien sûr, $FFT\left(\frac{n}{2}\right)$ sera obtenu récursivement en découpant en 2 blocs $FFT\left(\frac{n}{4}\right)$… et on poursuit jusqu'à avoir des cellules avec un seul item en entrée.

Pour finir, quand on a tous les $X_k$, on peut calculer les éléments du spectre par $\left|\frac{X_k}{n}\right|.

Vous pouvez poursuivre par quelques expériences sur les signaux suivants – clic-droit dessus et enregistrer pour télécharger.

Pour charger le signal, je vous donne le script Python :

def lecture_wave(filename) -> list:
    """
    Ouvre le fichier filename
    Renvoie le tableau formé par les échantillons du canal de gauche
    """
    file = wave.open(filename, 'rb')
    bytes_per_data   = file.getsampwidth()
    channels_number = file.getnchannels()
    frames_number   = file.getnframes()
    fe             = file.getframerate()
    bytesList      = file.readframes(framesNumber)
    file.close()
    bytes_per_frame = bytes_per_data * channels_number
    # que ce soit mono ou stéréo, je ne récupère qu'un canal
    canal = []
    for i in range(0, frames_number*bytesPerFrame, bytes_per_frame):
        frame_bytes = [bytesList[i + j] for j in range(bytes_per_data)]
        value = convert_to_float(frame_bytes)
        canal.append(value)
    return canal, fe


def convert_to_float(frame_bytes) -> float:
    """
    frame_bytes: liste d'octets, du poids faible au poids fort
    Convertit en réel
    Le codage est sur n octets en CA2,
      la valeur max correspond au flottant 1
      la valeur min au flottant -1
    """
    n = len(frame_bytes)
    is_positif = (frame_bytes[-1] & 0x80 == 0)
    if is_positif:
        value = 0
    else:
        value = -1
    for i in range(n):
        if is_positif:
            value += frame_bytes[i]
        else:
            value -= (~frame_bytes[i])&0xFF
        if i < n-1:
            value /= 256
        else:
            value /= 128
    return value