Traiter le signal audio numrique avec Matlab ou

Traiter le signal audio numérique avec Matlab ou Scilab Jean-Paul Stromboni, Polytech'Nice Sophia, Dpt Sciences Informatiques, SI 3 Durée 50 minutes, avec Matlab/Scilab, un vidéo projecteur, et des hauts parleurs Après ce chapitre, vous devrez savoir comment : • utiliser l’environnement de MATLAB : fenêtre de commande, éditeur de scripts, appel de fonction … et savoir passer à Scilab • utiliser les langages de scripts de MATLAB et de Scilab • exécuter des scripts et des fonctions MATLAB (et Scilab) • synthétiser des signaux audio numériques composés d’harmoniques et dotés d’enveloppes, reproduire le timbre d’instruments de musique • tracer des chronogrammes, des spectrogrammes, … Le TD n° 3 utilise Matlab ou Scilab pour : Synthétiser des signaux audio, représenter les chronogrammes, lire et écrire au format Wave, calculer et afficher le spectre, etc … 2011 -12, S. S. I. I. , SI 3, séance 3, Traiter le signal audio numérique avec Matlab ou Scilab, page 1

Utiliser MATLAB (ou Scilab) est un objectif du cours S. S. I. I. • MATLAB (pour MATrix LABoratory) de Mathworks, est un logiciel réputé dans le domaine du calcul scientifique. • Le département S. I. possède 25 licences MATLAB sur le réseau local de l’école à utiliser en travaux dirigés (cf. procédure d’installation) • Pour utiliser hors de l’école, on suggère d’installer également Scilab (à télécharger sur http: //www. scilab. org) outil libre, gratuit et multiplateforme très semblable à MATLAB (il est intéressant de comparer). • Pour apprendre MATLAB, on analyse des exemples tirés de : http: //www-gmm. insa-toulouse. fr/~guillaum/AM/ avec l’accord de Philippe Guillaume, Professeur à l'INSA de Toulouse, et auteur de l’ouvrage : 'Musique et Acoustique : de l’instrument à l’ordinateur', collection Hermès, éditeur Lavoisier. 2011 -12, S. S. I. I. , SI 3, séance 3, Traiter le signal audio numérique avec Matlab ou Scilab, page 2

Se repérer dans l’environnement MATLAB (celui de Scilab plus simple contient la plupart des outils de calcul et de tracé) dossier travail fenêtre plot Publish in html fenêtre edit fenêtre de commande. m file 2011 -12, S. S. I. I. , SI 3, séance 3, Traiter le signal audio numérique avec Matlab ou Scilab, page 3

Utiliser la fenêtre de commande de Matlab (traduire en Scilab) >> % en Matlab, le prompt est ’>>’ >> N=3 % avec ou sans caractère '; ' ? N = 3 >> N=3; >> n >> Message=[‘S. I. ‘, num 2 str(N), ’. ’]; >> disp(Message) >> help disp % aide succincte % il y a un éditeur de ligne de commande % on a droit aux commandes de shell >> pwd, ls, dir, cd % lancer une application MSDOS !notepad % exécuter le script MATLAB ’sinus. m’ >> sinus % sinus. m doit être dans le PATH >> format long % 10 chiffres décimaux % effacer la fenêtre Command Window >> clc % noter l’instruction eval >> eval([‘la’, num 2 str(3), ’=440’]) // En SCILAB, le prompt est ’ ’ // est un commentaire N=3; N n Message=[’date’, string(28), ’septembre’] disp(Message) help disp // il y a un éditeur de ligne commande pwd, dir, cd // lancer une application Windows ou Unix unix(‘notepad. exe’) exécuter un script ou une fonction Scilab exec(’sinus. sce’) // fichiers. sce et. sci sinus // %pi format(20) %pi clc execstr(’la 3= 440; ’); la 3 2011 -12, S. S. I. I. , SI 3, séance 3, Traiter le signal audio numérique avec Matlab ou Scilab, page 4

Analyse du script tiré du fichier ‘sinus. m’ et traduction en Scilab % arpège de sons sinusoïdaux ( Matlab) clear all; close all; Fe = 22050; h = 1/Fe; f 0 = 220; T = 1. 5; N = 13; fr = f 0*[1: 1: N]; am = 1; x = []; for k = 1: N % création d'une enveloppe tr = T*[0. 02. 98 1]; yr = [0 1 1 0]; env = envelop(tr, yr, Fe); % création du signal sinusoïdal th = 0: h: T; s = sin(2*pi*fr(k)*th); % application de l’enveloppe y = s. *env; % concaténation x = [x, am*y]; T = T*. 8; am = am*. 8; end plot(x) wavplay(x, Fe); wavwrite(x, Fe, '. /sinus. wav'); // arpège de sons sinusoïdaux ( Scilab) clear all; close all; Fe = 22050; h = 1/Fe; f 0 = 220; T = 1. 5; N = 13; fr = f 0*(1: N); am = 1; exec('envelop. sce'); x = []; for k = 1: N tr = T*[0. 02. 98 1]; yr = [0 1 1 0]; env = envelop(tr, yr, Fe); th = 0: h: T; y = sin(2*%pi*fr(k)*th). *env; x = [x, am*y]; T = T*. 8; am = am*. 8; end plot 2 d(x) sound(x, Fe); wavwrite(x, Fe, '. /scilabsinus. wav'); 2011 -12, S. S. I. I. , SI 3, séance 3, Traiter le signal audio numérique avec Matlab ou Scilab, page 5

Analyser la fonction Matlab ‘cloche’ tirée du fichier ‘cloche. m’ function s = cloche(f 1, T, Fe) % imitation d'une cloche % f 1 = fréquence du principal % Fe = fréquence d'échantillonnage % T = durée du son %--------------------h = 1/Fe; th = 0: h: T; f = f 1*[0. 5 1 1. 188 1. 530 2. 0000 2. 470 2. 607 2. 650 2. 991 3. 367 4. 137 4. 487 4. 829 5. 385 5. 863 6. 709 8. 077 8. 547 9. 017 9. 530. . . 11. 026 12. 393]; a = [350 950 500 150 700 100 250 370 1000 180 300 150 300 100 50 20 10 35 5 15]; s = synthad(a, f, 0*f, T, Fe); t = T*[0. 001. 4. 6. 9 1]; a = [0. 6 1. 4. 2. 1 0]; env = envelop(t, a, Fe); s = s. *env; End % optionnel Ligne d’en tête Commentaire accessible dans l’aide . . . vecteur des fréquences harmoniques vecteur des amplitudes composantes fréquentielles allure de la courbe d'enveloppe a(t) a 1 0. 6 0 0. 4 0. 6 t 2011 -12, S. S. I. I. , SI 3, séance 3, Traiter le signal audio numérique avec Matlab ou Scilab, page 6

Analyser la fonction synthad tirée du fichier synthad. m function s = synthad(a, f, p, T, Fe) % synthèse additive % cette fonction crée un son de duree T, % compose des partiels f(n), d'amplitude a(n) % et de phase a l'origine p(n). % Fe est la fréquence d'échantillonnage %-------------------% création du vecteur temps discret dt = 1/Fe; t = 0: dt: T; n = length(t); % création du son, une boucle ajoute une a une % les composantes fréquentielles s = zeros(1, n); K = length(f); for k = 1: K s = s+a(k)*sin(2*pi*f(k)*t+p(k)); end % normalisation pour que les valeurs % restent dans l'intervalle [-0. 99] s =. 99*s/max(abs(s)); cumul des harmoniques trouvés dans les vecteurs a : amplitude, f : fréquence et p : phase maximum de s ramené à 0. 99*max(abs(s)) 2011 -12, S. S. I. I. , SI 3, séance 3, Traiter le signal audio numérique avec Matlab ou Scilab, page 7

Analyser la fonction envelop utilisée dans sinus. m et cloche. m function env = envelop(t, a, Fe) % enveloppe parametree par t et a = env(t), % t contient une liste d'instants t_k % a contient la liste des amplitudes a_k aux instants t_k % env est le son echantillonne a la frequence Fe, % affine par morceaux, tel que env(t_k) = a_k %----------------------------lt = length(t); T = t(lt); h = 1/Fe; th = 0: h: T; % test validite de t if t(1) >= T error('t incompatible dans envelop'); end % test compatibilite t et a if lt ~= length(a) error('t et a de longueurs différentes'); end % au cas ou t ne serait pas strictement croissant : for k = 2: lt-1 if (t(k) <= t(k-1)) | (t(k) >= t(lt)) t(k) = (t(k-1)+t(lt))/2; end n = length(th); env = zeros(1, n); ni = lt-1; c = zeros(1, ni+1); b = c; h 2 = 0; for k = 1: ni h 1 = h 2+1; h 2 = 1+floor(t(k+1)/h); cb = [t(k) 1; t(k+1) 1][a(k) ; a(k+1)]; c = cb(1); b = cb(2); env(h 1: h 2) = c*th(h 1: h 2)+b; end env =. 99*env/max(env); a a(k+1) a(k) t(k+1) t a= c*t+b pour t(k) < t <t(k+1) Que valent c et b ? Résoudre a(k)=c*t(k)+b a(k+1)=c*t(k+1)+b AB calcule la solution x de A*x = B on s'en sert ici pour trouver les coefficients directeurs c et b de l'enveloppe entre t(k) et t(k+1) 2011 -12, S. S. I. I. , SI 3, séance 3, Traiter le signal audio numérique avec Matlab ou Scilab, page 8

Utiliser la fonction ‘cloche’ dans le script ‘gammes’ Que réalise le script gammes ? Comment procède t’on pour créer la gamme ? % tiré du fichier gammes. m clear all % clavier azerty : notes = ['a', 'z', 'e', 'r', 't', 'y', 'u', 'i', 'o', 'p', … 'q', 's', 'd']; Fe = 22050; f 0 = 440; % la 3 est la première note temp = 2. ^((0: 12)/12); % fr = f 0*temp; % fréquence des notes de la 3 à la 4 T = 1. 5; % durée des notes for k = 1: 13 % on crée des notes de flute note = notes(k); eval([note '= cloche(fr(k), T, Fe); ']); %ou flute end % et on joue : disp('pour jouer, rentrez une note parmi : '); disp('entrer a z e r t y u i o p q s ou d, suivi de enter'); x = 0; % taper x pour terminer note = [1, 1]; % length(note)= ? while length(note) >1 note = input('note suivante? (x pour terminer)'); if length(note) == 1 disp('termine'); break end soundsc(note, Fe); % met l'amplitude à l'échelle end efface tout vecteur ligne de caractères … continuateur de ligne ‘a’ note pour k=1 Pour k=1, eval exécute : a = cloche(fr(1), T, fe); soundsc ramène le signal entre -1 et le joue 2011 -12, S. S. I. I. , SI 3, séance 3, Traiter le signal audio numérique avec Matlab ou Scilab, page 9

Et si vous préférez la flute, comment utiliser le fichier ‘flute. m’? function s = flute(f 1, T, Fe) % son flute %--------------------------h = 1/Fe; th = 0: h: T; nt = length(th); a = [1000 50 80 10 5 2. 1 1]; nh = length(a); N = 1: nh; % synthèse additive du spectre f = N*f 1; s = synthad(a, f, 0*f, T, Fe); % enveloppe t = T*[0. 1. 2. 9 1]; a = [0. 8 1. 8 0]; env = envelop(t, a, Fe); s = s. *env; % ajout de souffle br = bruit(T, Fe); fn = f 1/Fe*2; wn = (fn*[. 8 1]). ^(2); b = fir 1(50, wn); br = filter(b, 1, br); t = T*[0. 05. 8 1]; a = [0 1. 5 0]; env = envelop(t, a, Fe); br = br. *env; s = s+br/7; s = s/(max(abs(s))); 2011 -12, S. S. I. I. , SI 3, séance 3, Traiter le signal audio numérique avec Matlab ou Scilab, page 10

Traduction de gammes. m en Scilab dans gammes. sce //-------------------------------// on fabrique une gamme chromatique au clavier //------------------------------clear all getf("envelop. sci"); getf("synthad. sci"); getf("cloche. sci"); // clavier azerty : notes = ['a', 'z' 'e' 'r' 't' 'y' 'u' 'i' 'o' 'p' 'q' 's' 'd']; Fe = 22050; f 0 = 440; // la 3 est la première note temp = 2. ^((0: 12)/12); // fr = f 0*temp; // fréquence des notes de la 3 à la 4 T = 1. 5; // durée des notes for k = 1: 13 // on crée des notes de flute note = notes(k); execstr( strcat([note, '= cloche(fr(k), T, Fe); '])); // ou flute, ou orgue end // et on joue : disp('pour jouer, rentrez une note parmi : '); disp('entrer a z e r t y u i o p q s ou d, suivi de enter'); x = 0; // taper x pour terminer note = a; // length(note)= ? while length(note) >1 note = input('note suivante ? '); if length(note) == 1 disp('termine'); break end sound(note, Fe); // met l'amplitude du son à l'échelle end La fonction getf obsolète est remplacée par exec("envelop. sci"); // par exemple 2011 -12, S. S. I. I. , SI 3, séance 3, Traiter le signal audio numérique avec Matlab ou Scilab, page 11

Utiliser le type ‘cell’ de Matlab, exemple de ‘creegammes. m’ %% créer des notes et des airs avec les cellules gamme={'do', 'dod', 'red', 'mi', 'fad', 'sold', 'lad', 'si'}; dt=power(2, 1/12); for g=1: 5 frla=110*2^(g-1); for n=1: length(gamme) eval([gamme{n}, num 2 str(g), '=frla*dt^(n-10); ']) end %%créer une mélodie jouer(la 3, 1); jouer(si 3, 0. 75); jouer(dod 4, . 3, 0. 5); jouer(mi 4, . 3, . 5); jouer(re 4, . 3, . 5); jouer(fad 4, . 3, 1); jouer(mi 4, . 6, 1); jouer(la 4, 1); jouer(sold 4, 1); jouer(la 4, . 8, 1); 2011 -12, S. S. I. I. , SI 3, séance 3, Traiter le signal audio numérique avec Matlab ou Scilab, page 12

Voici la fonction jouer utilisée par le script creegammes (jouer. m) function note=jouer(fr, Dur, amp, Fs) % jouer possède 4 arguments, fr est la fréquence de la note % dur sa durée en seconde, ampl son amplitude initiale, et Fs % la fréquence d'échantillonnage. L'enveloppe est linéaire. % Par défaut, fr= 440 Hz, Dur= 1 s, amp= 1, et Fs= 8000 Hz. % En l'absence d'argument de sortie, la note créée est jouée. f=440; D=1; a=1; fe=8000; switch nargin case 1 f=fr; case 2 f=fr; D=Dur; case 3 f=fr; D=Dur; a=amp; case 4 f=fr; D=Dur; a=amp; fe=Fs; end % construire la note t=[0: 1/fe: D]; note=a*sin(2*pi*f*t). *(1 -t/D); if nargout==0 wavplay(note, fe) end Traduction Scilab de jouer. m dans jouer. sce : Function note=jouer(fr, Dur, amp, Fs) // fr est la fréquence de la note, Dur est sa durée en seconde, // ampl est son amplitude, Fs la fréquence d'échantillonnage. // enveloppe linéaire, fr=440 Hz, Dur=1 s, amp=1, Fs=8 k. Hz nbin=argn(2); //nbout=argn(1); //marche pas, toujours égal à 1 ? ? ? f=440; D=1; a=1; fe=8000; select nbin case 1 then f=fr; case 2 then f=fr; D=Dur; case 3 then f=fr; D=Dur; a=amp; case 4 then f=fr; D=Dur; a=amp; fe=Fs; end t=[0: 1/fe: D]; note=a*sin(2*%pi*f*t). *(1 -t/D); // tester nbout pour reproduire nargout endfunction 2011 -12, S. S. I. I. , SI 3, séance 3, Traiter le signal audio numérique avec Matlab ou Scilab, page 13

calculer et afficher le spectre avec Matlab %% paramétrage N=1024 Fe=44100 D=2; f 1=880; %% son de cloche s=cloche(f 1, D, Fe); wavplay(s, Fe) %% spectre fr=[0: N-1]*Fe/N; sp=abs(fft(s, N))/N; plot(fr(1: N/2), sp(1: N/2)) grid xlabel('fréquence Hz') ylabel('spectre') title(['spectre', num 2 str(N), . . . ' points']) 2011 -12, S. S. I. I. , SI 3, séance 3, Traiter le signal audio numérique avec Matlab ou Scilab, page 14

Afficher le spectrogramme avec Matlab %% spectrogramme d'un arpège s=cloche(f 1, D, Fe); s=[s, cloche(5*f 1/8, D, Fe), . . . cloche(3*f 1/4, D, Fe)]; wavplay(s, Fe) spectrogram(s, N, 0, N, Fe, 'yaxis') colorbar % cf. spectrogramme ci-contre %% idem pour un accord acc= cloche(f 1, D, Fe)+. . . cloche(5*f 1/4, D, Fe)+. . . cloche(3*f 1/2, D, Fe); acc= 0. 99*acc/max(abs(acc)); wavplay(acc, Fe) spectrogram(acc, N, 0, N, Fe, 'yaxis') colorbar % non tracé 2011 -12, S. S. I. I. , SI 3, séance 3, Traiter le signal audio numérique avec Matlab ou Scilab, page 15

Traduire le fichier creegammes. m en Scilab (creegammes. sce) //créer des notes et des airs avec les cellules gamme=['do', 'dod', 'red', 'mi', 'fad', 'sold', 'lad', 'si']; [nl, nc]=size(gamme); dt=2^(1/12); for g=1: 5, frla=110*2^(g-1); for n=1: nc, execstr([gamme(n)+string(g)+'=frla*dt^(n-10); ']) end //créer une mélodie exec('jouer. sce'); s=jouer(la 3, 1); s=[s, jouer(si 3, 0. 75)]; s=[s, jouer(dod 4, . 3, 0. 5)]; s=[s, jouer(mi 4, . 3, . 5)]; s=[s, jouer(re 4, . 3, . 5)]; s=[s, jouer(fad 4, . 3, 1)]; s=[s, jouer(mi 4, 1)]; s=[s, jouer(la 4, 1)]; s=[s, jouer(sold 4, 1)]; s=[s, jouer(la 4, 1)]; sound(s, 8000); savewave('majoie. wav', s, 8000) 2011 -12, S. S. I. I. , SI 3, séance 3, Traiter le signal audio numérique avec Matlab ou Scilab, page 16

Traduction en Scilab des deux functions ‘synthad’ et ‘envelop’ function s = synthad(a, f, p, T, Fe) // synthese additive // cette fonction cree un son de duree T, // compose des partiels f(n), d'amplitude a(n) // et de phase a l'origine p(n). // Fe est la frequence d'echantillonnage //----------------------// creation du vecteur temps discret dt = 1/Fe; t = 0: dt: T; n = length(t); // creation du son, boucle pour ajouter une a une // les composantes frequentielles s = zeros(1, n); K = length(f); for k = 1: K s = s+a(k)*sin(2*%pi*f(k)*t+p(k)); end // normalisation pour que les valeurs soient // toutes dans l'intervalle [-0. 99] s =. 99*s/max(abs(s)); endfunction [env] = envelop(t, a, Fe) lt = length(t); T = t(lt); h = 1/Fe; th = 0: h: T; if t(1) >= T // test de validite de t error('t incompatible dans envelop'); end if lt ~= length(a) // test de compatibilité de t et a error('t et a de longueur différente dans envelop'); end // au cas où t ne serait pas strictement croissant : for k = 2: lt-1 if (t(k) <= t(k-1)) | (t(k) >= t(lt)) t(k) = (t(k-1)+t(lt))/2; end n = length(th); env = zeros(1, n); ni = lt-1; c = zeros(1, ni+1); b = c; h 2 = 0; for k = 1: ni h 1 = h 2+1; h 2 = 1+floor(t(k+1)/h); cb = [t(k) 1; t(k+1) 1][a(k) ; a(k+1)]; c = cb(1); b = cb(2); env(h 1: h 2) = c*th(h 1: h 2)+b; end env =. 99*env/max(env); endfunction 2011 -12, S. S. I. I. , SI 3, séance 3, Traiter le signal audio numérique avec Matlab ou Scilab, page 17

Traduire pour Scilab le script Matlab dans ‘cloche. m’ function s = cloche(f 1, T, Fe) // imitation d'une cloche // f 1 = fréquence fondamentale // Fe = fréquence d‘échantillonnage // T = durée du son //----------------------h = 1/Fe; th = 0: h: T; f = f 1*[0. 5 1 1. 188 1. 530 2. 0000 2. 470 2. 607 2. 650 2. 991. . . 3. 367 4. 137 4. 487 4. 829 5. 385 5. 863 6. 709 8. 077. . . 8. 547 9. 017 9. 530 11. 026 12. 393]; a = [350 950 500 150 700 100 250 370 1000 180 300. . . 100 150 300 100 50 20 10. . . 35 5 15]; s = synthad(a, f, 0*f, T, Fe); t = T*[0. 001. 4. 6. 9 1]; a = [0. 6 1. 4. 2. 1 0]; env = envelop(t, a, Fe); s = s. *env; endfunction 2011 -12, S. S. I. I. , SI 3, séance 3, Traiter le signal audio numérique avec Matlab ou Scilab, page 18

Analyser le script ci-dessous, et répondre aux questions posées function [sdet, fe]=note(fr, dure, a, fe) % la fonction permet de créer une note sinusoïdale d’amplitude a % de fréquence fr et de durée dure avec la fréquence d’échantillonnage fe % par défaut, fr=440 Hz, si a n'est pas précisé, a=0. 75, % dure=1 si non précisé, si fe n'est pas précisé, fe=44100 Hz % une enveloppe exponentielle est attachée à la note créée % la note est jouée s'il n'y a pas d'argument de sortie % Noter : ce commentaire apparaît dans l’aide : help note if nargin==0, fr=440; dure=1; a=1; fe=44100; end if nargin ==1, dure=1; a=1; fe=44100; end if nargin ==2, a=1; fe=44100; end if nargin ==3, fe=44100; end temps=[0: 1/fe : dure]; env=exp(-temps/(dure/3)); sdet= a*cos(2*pi*temps*fr). *env; if nargout==0, wavplay(sdet, fe), stem(temps, sdet) xlabel('temps (s)'), ylabel('note'), grid end % optionnel note. m Quel doit être le nom du fichier qui le contient ? On saisit >> note; % dans la fenêtre de commande, quel est le résultat ? On saisit s=note; % Quel est le résultat ? Comment utiliser note pour créer un arpège puis pour créer un accord ? arpege= [note(440), note(5*440/4), note(3*440/2), note(880, 2, . 5)]; accord= note(440)+note(5*440/4)+note(3*440/2); 2011 -12, S. S. I. I. , SI 3, séance 3, Traiter le signal audio numérique avec Matlab ou Scilab, page 19

L’éditeur de scripts de MATLAB • • Taper edit, ou File/New/m-file, ou utiliser la barre d’outils. ‘%’ est le commentaire ‘%%’ permet de découper le script en cellules exécutables une par une Publish in html permet de faire un compte rendu en html des résultats du Execute Cell script (images comprises) F 9 evaluate selection F 5 run L’exécution du script ead. m depuis command window : >> ead seulement si le dossier contenant ead. m est dans le Path, ou dans le current directory File/Set Path, modifie le path : Dock Publish html – suivi de Add Folder – puis Save, et Close • En Scilab, faire exec(‘fichier. Script’). – Fichiers. sce, ou. sci, ou. txt, ou … – Les fichiers scripts doivent être dans le répertoire de travail, ou il faut préciser leur chemin 2011 -12, S. S. I. I. , SI 3, séance 3, Traiter le signal audio numérique avec Matlab ou Scilab, page 20

Ajouter une fonction nouvelle dans MATLAB Matlab contient déjà des fonctions, telles que plot(), wavread() … l'utilisateur peut en ajouter, qui doivent respecter les mêmes règles : • Si la fonction est nommée 'notepure', le fichier script doit se nommer 'notepure. m' • notepure. m doit se trouver dans le Path !! • On appelle ‘notepure’ depuis un script ou après le prompt de commande avec : – notepure(‘la 3. wav’, 0. 5, 440, 1. 5); – s = notepure(‘si. wav’, 0. 5, 440, 1. 5); • Le nombre d'arguments donnés à l'appel d'une fonction peut varier, c'est Matlab qui gère : – plot(t, s), – plot(s), – ou encore h=plot(t, s 1, t, s 2) • Les deux variables 'nargin' et 'nargout' transmettent à la fonction appelée le nombre d'arguments d’entrée et de sortie spécifiés lors de l'appel. • Dans Scilab, il faut exécuter l’instruction exec(‘mon. Fichier. Script’) pour pouvoir appeler les fonctions définies dans mon. Fichier. Script • la dernière ligne ‘endfunction’ est obligatoire 2011 -12, S. S. I. I. , SI 3, séance 3, Traiter le signal audio numérique avec Matlab ou Scilab, page 21