====== Algorithme FFT - Fast Fourier Transform ====== {{ .:fourier.jpg?nolink&200|}} 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 [[https://fr.wikipedia.org/wiki/Traitement_du_signal|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. ===== Idée d'une superposition ===== ==== Corde vibrante ==== 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. {{ .:cordes_fourier.jpg?direct&400 |}} 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$... ==== Signal ==== 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. ==== Fréquence pure ==== 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. {{ youtube>bopZQdsJb_E }} {{ .:sin440.png?direct&400 |}} ==== Superposition ==== 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 : {{ .:son_composite.png?direct&400 |}} ==== Percevoir les 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é. ==== La transformation de Fourier ==== 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$$ {{ .:spectre.svg |}} 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. Sans rapport avec notre projet mais intéressant : Dans le [[https://www.seuil.com/ouvrage/traite-des-objets-musicaux-pierre-schaeffer/9782020026086|Traité des objets musicaux]], Pierre Schaeffer montre que l'amorce du son compte plus que le spectre. C'est le « ffffff » caractéristique du souffle au début d'une note de flûte qui fait que l'on identifie immédiatement le son d'une flûte. Même si le « ffffff » ne dure qu'une fraction de seconde, notre oreille le perçoit et lui donne une grande importance. Cette situation est typique : nous sommes capable d'enregistrer des sons avec une grande précision, nous pouvons visualiser des signaux, faire des calculs très savants, mais tout cela ne mène pas loin si on n'a pas identifié d'abord quels étaient les éléments importants du signal. ==== Échantillonnage ==== 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. {{ .:dotted_graphe.png?direct&400 |}} 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 [[https://fr.wikipedia.org/wiki/Th%C3%A9or%C3%A8me_d%27%C3%A9chantillonnage|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$. ===== La FFT ===== ==== Entrées / sorties ==== 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. Python possède un type ''complex''. >>> complex(3,2) 3+2j >>> z = complex(3,2) >>> abs(z) 3.6055512754639896 Remarquez que l'on utilise ''j'' pour le nombre tel que $j^2 = -1$ car ''i'' entraînerait trop de confusions. * Le tableau contient ''n'' items, autant que ''signal''. * L'item de rang ''i'' correspondant à la fréquence ''i * Fe/n = i * 1/T''. Les fréquences apparaissant dans ''spectre'' sont donc de la forme $f = i\frac{F_e}{n}$. Dans notre exemple, nous avions un signal avec la fréquence de base -- //fondamental// -- à $440\,Hz$ et les autres -- //harmoniques// -- à $880\,Hz$, $1320\,Hz$ et $1760\,Hz$. En fonction du choix de $F_e$, rien ne garantit que ces 4 fréquences soient de la forme $i\frac{F_e}{n}$ ! Que va-t-il se passer ? ''spectre'' va-t-il rater les fréquences de notre signal au motif qu'on ne regarde pas précisément les bonnes fréquences ? Non, bien sûr ! Par exemple, si ''Fe = 4000'' et que l'on prélevé un son de ''T = 0.01'' secondes, on obtiendra un tableau ''signal'' de ''n = T * Fe = 40''. Le pas de fréquence est ''Fe/n = 1/T = 100'' et dans le tableau ''spectre'', les indices ''0, 1, 2, ...'' correspondent aux fréquences ''0, 100, 200, ...''. Notre fréquence qui est à $440\,Hz$ se répartira donc sur les indices voisins, surtout sur les plus proches. * 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 FFT est calculée en utilisant une version complexe des fonctions trigonométriques : $\exp(i\cdots) = \cos(\cdots) + i\cdot \sin(\cdots)$ ce qui entraîne $\cos(\cdots) = \frac{1}{2}\exp(i\cdots) + \frac{1}{2}\exp(-i\cdots)$. Le spectre obtenu par ''fft'' nous donne les intensités des $\exp(i\theta)$. Ainsi, si on a par exemple $3\cos(\cdots) = \frac{3}{2} \exp(i\cdots) + \frac{3}{2} \exp(-i\cdots)$, alors dans ''spectre'' on aura une valeur de ''3/2'' au rang ''i'' correspondant et au rang ''n-i''. ===== L'algorithme ===== ==== Version normale, trop lente ==== 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 **À faire :** implémentez cette fonction. Produisez le tableau ''signal'' pour le cas de $$f(t) = \sin(2\pi\cdot 440 t) + 0,3 \sin(2\pi\cdot 880 t) + 0,5 \sin(2\pi\cdot 1320 t) + 0,1 \sin(2\pi\cdot 1760 t)$$ Choisissez ''T = 0.1'' et ''Fe = 4000''. Représentez graphiquement ce signal avec ''matplotlib.pyplot''. Faites les calcule de toutes les raies du spectre pour obtenir le tableau ''spectre''. Représentez graphiquement le contenu de ''spectre'' avec ''matplotlib.pyplot''. Vous devriez voir apparaître nettement le spectre avec les raies. 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//. ==== Version rapide, fft ==== 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. {{ .:fft_cell.png?direct&400 |}} 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. {{ .:fft_2cells.png?direct&400 |}} 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. {{ .:fft_cell_0.png?direct&400 |}} Pour finir, quand on a tous les $X_k$, on peut calculer les éléments du spectre par $\left|\frac{X_k}{n}\right|. **À faire :** implémentez l'algorithme FFT. Testez sur le signal $f$ déjà rencontré. Faites un comparatif de temps d'exécution avec la version simple. ===== Signaux réels ===== Vous pouvez poursuivre par quelques expériences sur les signaux suivants -- //clic-droit dessus et enregistrer pour télécharger.// {{ .:a.wav |a}} {{ .:e.wav |e}} {{ .:i.wav |i}} {{ .:o.wav |o}} {{ .:u.wav |u}} {{ .:aeiaeo.wav |aeiaeo}} 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