'Profil mélodique' du dernier mouvement de Herz und Mund und Tat und Leben (J.S. Bach, Jésus que ma joie demeure).Qui n'a jamais eu un air en tête sur lequel il est incapable de mettre un nom? Pour le non musicien, c'est d'autant plus frustrant qu'il est incapable de noter la mélodie, qu'il risque de l'oublier et ne peut la communiquer à d'autres qu'en la chantonnant, la sifflant, etc... Dans des circonstances qui ne sont jamais favorables. Le résultat est souvent désastreux...
Le
code Parsons propose un remède à ce problème. Inventé par
Denys Parsons en 1975, c'est une notation musicale accessible au non-initié. ...Enfin presque.
Le code correspondant à une mélodie est constitué d'une suite de quatre symboles :
- '
*', la première note de la mélodie
- '
U', une note plus 'haute' (aigüe) que la précédente (Up)
- '
D', une note plus 'basse' que la précédente (Down)
- '
R', une note de la même hauteur que la précédente (Repeat)
Ceci, indépendamment du rythme, de l'écart entre les notes, etc... Cela ne décrit donc pas la mélodie avec précision, ce code est une espèce de 'signature'. Il est impossible de jouer la mélodie en décodant le code Parsons. À la limite, on peut essayer de deviner. Par contre, on peut vérifier assez facilement que le code Parsons correspond bien à la mélodie.
Denys Parsons avait compilé plus de 10 000 mélodies dans un ouvrage appelé
Directory of Tunes and Musical Themes en 1975 (et ré-édité/amendé en 2010 chez Piatkus Books sous le nom
The Directory of Classical Themes. Il s'agit d'une espèce de dictionnaire de mélodies triées par code.
À partir du code, ce bouquin ou le site
Musipedia me permet de mettre un nom sur la mélodie. Donc, c'est génial. Sauf que, pour un non-initié de mon acabit, il est difficile de décider laquelle de deux notes qui se suivent est la plus grave ou même s'il s'agit de la même note. Sauf dans les cas les plus évidents. Et, encore... Par exemple, je suis incapable de trouver *UUDRUUDUUUDUU à partir de
Frère Jacques.
C'est là que l'informatique vient à mon secours. Si j'enregistre mon sifflement, sous Linux,
sndfile-spectrogram (de
sndfile-tools) me permet de visualiser la hauteur des notes :
On voit que, sifflé,
Jésus que ma joie demeure donne un beau
sonagramme bien net et relativement facile à décoder. On trouve assez facilement
uuudrudrududdduudud que la recherche par contour mélodique sur
http://musipedia.org trouve facilement.
On voit que le sifflement se passe du côté de 1000 Hz (une octave plus haut que le
La 440 Hz); il était donc inutile d'échantillonner à 22.05 Khz; échantillonné à 4 KHz, cela donne :
C'est bien, mais ce serait encore mieux d'écrire un petit programme qui génère le code à partir du sifflement (comme le fait le site si on active le micro en FLASH; mauvaise idée...). En fait, il 'suffit' de calculer la transformée de Fourier du signal audio capté par le micro. Reste à trouver le code pour accéder au micro et le code pour faire une FFT.
*
Parrot.c, un petit programme de démonstration d'accès à l'audio sous Linux
*
Discrete Fourier Transform à pied un tutorial sur les FFT comprenant un code simple.
Accès aux données audio sous Linux
Pour pouvoir analyser un son, il faut accéder à l'échantillonnage du signal. /dev/dsp nous donne accès aux données brutes du
PCM. Parrot.c montre comment initialiser le périphérique pour y lire 8000 échantillons de 8 bits par seconde. On en dérive un petit programme (minimum) très simple qui, suivant qu'il s'appelle 'r' (record) ou 'p' (play), copie les octets échantillonnés de /dev/dsp dans la sortie standard (utilisation : $ r > mon_enregistrement.raw) ou 'joue' les octets échantillonnés de l'entrée standard dans /dev/dsp (utilisation $ p < mon_enregistrement.raw). Il n'en faut pas plus.
#include <stdlib.h>
#include <stdio.h>
#include <fcntl.h>
#include <sys/ioctl.h>
#include <linux/soundcard.h>
#define SAMPLE_RATE 8000 /* the sampling rate */
#define SAMPLE_SIZE 8 /* sample size: 8 or 16 bits */
#define CHANNELS 1 /* 1 = mono 2 = stereo */
typedef unsigned char uchar;
uchar audio_buf[1024];
int main(int argc, char *argv[])
{
int fd, n;
int arg; /* argument for ioctl calls */
if ((fd = open("/dev/dsp", O_RDWR)) < 0)
{
perror("open of /dev/dsp failed");
exit(1);
}
arg = SAMPLE_SIZE; /* sample size */
if (ioctl(fd, SOUND_PCM_WRITE_BITS, &arg) < 0)
perror("SOUND_PCM_WRITE_BITS ioctl failed");
if (arg != SAMPLE_SIZE)
perror("unable to set sample size");
arg = CHANNELS; /* mono or stereo */
if (ioctl(fd, SOUND_PCM_WRITE_CHANNELS, &arg) < 0)
perror("SOUND_PCM_WRITE_CHANNELS ioctl failed");
if (arg != CHANNELS)
perror("unable to set number of channels");
arg = SAMPLE_RATE; /* sampling rate */
if (ioctl(fd, SOUND_PCM_WRITE_RATE, &arg) < 0)
perror("SOUND_PCM_WRITE_RATE ioctl failed");
if (argv[0][0] == 'p')
{ /* play stdin into /dev/dsp */
while ((n = read(0, audio_buf, sizeof(audio_buf))) > 0)
write(fd, audio_buf, n);
}
else if (argv[0][0] == 'r')
{ /* record /dev/dsp into stdout */
while ((n = read(fd, audio_buf, sizeof(audio_buf))) > 0)
write(1, audio_buf, n);
}
exit(0);
}
On pourrait aussi se servir des outils standards, l'enregistreur de son et le 'couteau suisse'
Sox (SOund eXchange) pour convertir les différents formats audio de/vers un 'raw' 8 bits à 8 KHz (ce qui est suffisant pour l'usage qu'on en fait). L'avantage du raw 8 bits, c'est que l'on n'a pas besoin de se préoccuper des en-tête, d'un autre codage,... Les outils font ça très bien; moi, tout ce que je veux faire, c'est analyser le signal audio pour voir si je peux détecter la différence de hauteurs (fréquences) entre deux notes sifflées dans le micro. (Si je veux communiquer les fichiers audio à quelqu'un qui n'a que des outils standard, j'utilise
LAME pour produire un .mp3)
Je peux m'assurer que tout cela fonctionne en faisant un petit test (j'enregistre 'Frère Jacques', je transforme mon fichier .raw en un .wav pour pouvoir en faire un sonagramme avec
sndfile-spectrogram puis, je 'joue' le morceau tel que je l'ai enregistré; l'enregistrement est interrompu par ^C (
-C))):
$ ./r > FJ.raw
$ sox -v 2.0 -r 8k -b 8 -e unsigned -c 1 FJ.raw FJ.wav
$ sndfile-spectrogram FJ.wav 640 480 FJ.png
$ ./p < FJ.raw
Il faut noter que la gestion du son sur Linux n'est pas des plus simple. ALSA a remplacé OSS et le contrôle se fait via divers programmes comme alsamixer, alsactl, gnome-audio-contrl,... (et, par exemple, il est difficile de savoir où se trouve la pré-amplification du micro dans une Ubuntu 8.04)
Traitement du signal
'On peut montrer que tout signal périodique est décomposable en une somme de sinusoïdes' ...ou quelque chose comme ça. Quelles sinusoïdes?, avec quelles intensités?, comment les calculer?
FFT à pied explique que
* les fréquences sont [1..N/2] * fréquence_échantillonnage/N
* l'intensité est la somme des produits de la mesure par le sinus/cosinus à 'cet endroit'
C'est assez simple, mais il existe une méthode plus efficace qui élimine des calculs redondants : la FFT. C'est cette méthode qui est utilisée dans la pratique. On passe d'un temps calcul en N² en Nlog(N); pour N=1024, on passe de 1M à 10K (on boucle 128 fois moins (mais la boucle est un peu plus complexe)).
J'aurais sans doute pu utiliser KissFFT ou un autre librairie plus ou moins généraliste et optimisée. Mais, pour faire simple, j'utilise le code court et auto-suffisant de Discrete Fourier Transform à pied :
//
// Listing 1.4: The Discrete Fast Fourier Transform (FFT):
//
#define M_PI 3.14159265358979323846
void smbFft(float *fftBuffer, long fftFrameSize, long sign)
/*
FFT routine, (C)1996 S.M.Bernsee. Sign = -1 is FFT, 1 is iFFT (inverse)
Fills fftBuffer[0...2*fftFrameSize-1] with the Fourier transform of the time
domain data in fftBuffer[0...2*fftFrameSize-1]. The FFT array takes and returns
the cosine and sine parts in an interleaved manner, ie.
fftBuffer[0] = cosPart[0], fftBuffer[1] = sinPart[0], asf. fftFrameSize must
be a power of 2. It expects a complex input signal (see footnote 2), ie. when
working with 'common' audio signals our input signal has to be passed as
{in[0],0.,in[1],0.,in[2],0.,...} asf. In that case, the transform of the
frequencies of interest is in fftBuffer[0...fftFrameSize].
*/
{
float wr, wi, arg, *p1, *p2, temp;
float tr, ti, ur, ui, *p1r, *p1i, *p2r, *p2i;
long i, bitm, j, le, le2, k, logN;
logN = (long)(log(fftFrameSize)/log(2.)+.5);
for (i = 2; i < 2*fftFrameSize-2; i += 2) {
for (bitm = 2, j = 0; bitm < 2*fftFrameSize; bitm <<= 1) {
if (i & bitm) j++;
j <<= 1;
}
if (i < j) {
p1 = fftBuffer+i; p2 = fftBuffer+j;
temp = *p1; *(p1++) = *p2;
*(p2++) = temp; temp = *p1;
*p1 = *p2; *p2 = temp;
}
}
for (k = 0, le = 2; k < logN; k++) {
le <<= 1;
le2 = le>>1;
ur = 1.0;
ui = 0.0;
arg = M_PI / (le2>>1);
wr = cos(arg);
wi = sign*sin(arg);
for (j = 0; j < le2; j += 2) {
p1r = fftBuffer+j; p1i = p1r+1;
p2r = p1r+le2; p2i = p2r+1;
for (i = j; i < 2*fftFrameSize; i += le) {
tr = *p2r * ur - *p2i * ui;
ti = *p2r * ui + *p2i * ur;
*p2r = *p1r - tr; *p2i = *p1i - ti;
*p1r += tr; *p1i += ti;
p1r += le; p1i += le;
p2r += le; p2i += le;
}
tr = ur*wr - ui*wi;
ui = ur*wi + ui*wr;
ur = tr;
}
}
}
Ce code, très simple et évident (en ce sens qu'il ne s'agit 'que' de quelques boucles imbriquées et d'affectations), ne semble pas optimisable de manière flagrante. Ce n'est pas très long, c'est vraisemblablement le cas d'école de la FFT, l'utilisation est évidente : 3 paramètres et cela calcule la FFT .
Premier test
Comme premier test, on peut traiter 1024 échantillons d'un sifflement (ou deux) et calculer la magnitude du signal aux fréquences calculées par la FFT.
unsigned char audio_buf[1024];
float fft_buf[2*1024]; /* alternance réel+imaginaire */
void bufcpy(float *fp, uchar *ucp, int len)
{ /* mise en place des mesures audio_buf dans le fft_buffer */
while (--len >= 0)
{
*fp++ = (float)*ucp++;
*fp++ = 0.0; /* i-part */
}
}
void print_mag(float *fp, int len)
{
float r_f, i_f;
while (--len >= 0)
{
r_f = *fp++;
i_f = *fp++;
fprintf(stdout, "%f\n", sqrt(r_f*r_f+i_f*i_f));
}
}
main()
{ /* ~pseudo-code */
...
read(fd, audio_buf, sizeof(audio_buf));
bufcpy(fft_buf, audio_buf, 1024);
smbFft(fft_buf, 1024, -1);
print_mag(fft_buf, 1024);
}
On obtient quelque chose comme :
On repère bien deux pics à certaines fréquences (mais l'abscisse est en 'indices de buffer'; pour connaître la fréquence correspondante, il faut multiplier par 8000/512). L'indice 0 a été supprimé (il correspond au niveau de base (?)) et on remarquera que l'information utile se trouve entre 0 et 512, le reste étant un calcul miroir (c'est produit par la méthode de calcul (FFT); de toute façon on ne calcule pas de sinus de fréquence supérieure à la moitié de la fréquence d'échantillonnage (Théorème d'échantillonnage de Nyquist-Shannon)
Il faut noter que cette échelle est linéaire alors que notre perception musicale ne l'est pas; l'oreille est sensible aux rapports entre les notes : une octave correspond à un doublement de fréquence. L'octave la plus grave du piano va de 32.70 à 65.41 Hz (environ 33 Hz pour 12 notes (les 7 blanches et 5 noires du piano - la gamme tempérée); la plus aigüe de 2093,00 Hz à 4186,01 Hz (2093 Hz pour 12 notes) (ref.). On a donc intérêt à siffler 'aigu', il y aura plus de Herz entre les notes, elles seront plus faciles à individualiser.
Contour mélodique
En effectuant une boucle et en n'imprimant que l'indice de la fréquence maximum, on obtient une espèce de sonagramme simplifié, un contour mélodique.
int get_max(float *fp, int len)
{
int i, imax, n;
float max, mean, mag2;
mean = max = 0.0; n = 0;
for (i = 1; i < len; i++)
{
mag2 = fp[2*i] * fp[2*i] + fp[2*i+1] * fp[2*i+1];
if (mag2 > max)
{
imax = i;
max = mag2;
}
mean += mag2; n += 1;
}
if (max < (3.0 * mean/n))
imax = 0; /* no signal */
return(imax);
}
void print_max(float *fp, int len)
{
fprintf(stdout, "%d\n", get_max(fp, len));
}
main()
{
... ~ pseudo-code
for (;;)
{
read(fd, audio_buf, sizeof(audio_buf));
bufcpy(fft_buf, audio_buf, 1024);
smbFft(fft_buf, 1024, -1);
print_max(fft_buf, 512); /**/
}
}
Par exemple, pour Frère Jacques, on obtient quelque chose comme
On remarque que la cinquième note n'est pas exactement au niveau de la quatrième (alors qu'elle devrait...). Il faudra donc être un peu tolérant (ou siffler plus juste).
Cela s'annonce donc bien. On peut noter que l'on peut jouer sur deux paramètres : la fréquence d'échantillonnage et taille du buffer. Dans la configuration actuelle avec 8000 échantillons par seconde et un buffer de 1024 échantillons, on ne calcule que 8 maxima par seconde avec une précision d'environ 8 Hz. Il faut donc maintenir la note suffisamment longtemps pour que 'ce soit clair'. Si on réduit la taille du buffer, la discrimination des fréquences sera moins fine. Si on augmente la fréquence d'échantillonnage (...à déterminer)
Le code Parsons
En essayant de décoder la succession des notes par
void print_parsons(float *fp, int len)
{
static int last_one, sound;
int imax;
imax = get_max(fp, len);
if (imax < 100)
{
sound = 0;
return;
}
if (sound)
return; /* already processed */
if (last_one == 0)
fprintf(stdout, "*");
else {
if (imax > (last_one + 10))
fprintf(stdout, "u");
else if (imax < (last_one - 10))
fprintf(stdout, "d");
else fprintf(stdout, "r");
}
last_one = imax;
sound = 1;
fflush(stdout);
}
On obtient le code Parsons au fur et à mesure.
Contour mélodique schématique
C'est curieusement l'algorithme le plus sophistiqué...
On peut voir le code ici.
Entendre et voir de quoi l'on parle...
Une courte présentation multi-média (assez primitive) pour fixer les idées.
Elle a été réalisée avec
Imagination 2.1.
