1 / 21

Traiter le signal audio numérique avec Matlab ou Scilab

Jean-Paul Stromboni, Polytech'Nice Sophia, Dpt Sciences Informatiques, SI3 Durée 50 minutes, avec Matlab/Scilab, un vidéo projecteur, et des hauts parleurs. Après ce chapitre, vous devrez savoir comment :. Traiter le signal audio numérique avec Matlab ou Scilab.

shanae
Download Presentation

Traiter le signal audio numérique avec Matlab ou Scilab

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Jean-Paul Stromboni, Polytech'Nice Sophia, Dpt Sciences Informatiques, SI3 Durée 50 minutes, avec Matlab/Scilab, un vidéo projecteur, et des hauts parleurs Après ce chapitre, vous devrez savoir comment : Traiter le signal audio numérique avec Matlab ou Scilab 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 spectres, 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 …

  2. 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 multi-plateforme 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 dePhilippe Guillaume, Professeur à l'INSA de Toulouse, et auteur de l’ouvrage : 'Musique et Acoustique : de l’instrument à l’ordinateur', collection Hermès, éditeur Lavoisier.

  3. 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

  4. Utiliser la fenêtre de commande de Matlab (traduire en Scilab)  // En SCILAB, le prompt est ’’  // est un commentaire • N=3  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(’la3= 440;’); • la3 >> % en Matlab, le prompt est ’>>’ >> N=3 % avec ou sans caractère ';' ? N = 3 >> N=3; >> n >> Message=[‘S.I.‘,num2str(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’,num2str(3),’=440’])

  5. Analyse du script tiré du fichier ‘sinus.m’ et traduction en Scilab // arpège de sons sinusoïdaux (Scilab) clear all; close all; Fe = 22050; h = 1/Fe; f0 = 220; T = 1.5; N = 13; fr = f0*(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 plot2d(x) sound(x,Fe); wavwrite(x,Fe,'./scilabsinus.wav'); % arpège de sons sinusoïdaux (Matlab) clear all; close all; Fe = 22050; h = 1/Fe; f0 = 220; T = 1.5; N = 13; fr = f0*[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 (Matlab) clear all; close all; Fe = 22050; h = 1/Fe; f0 = 220; T = 1.5; N = 13; fr = f0*[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');

  6. Analyser la fonction Matlab ‘cloche’ tirée du fichier ‘cloche.m’ function s = cloche(f1,T,Fe) % s = cloche(f1,T,Fe) % imitation d'une cloche % f1 = fréquence du principal % Fe = fréquence d'échantillonnage % T = durée du son %----------------------------------------- h = 1/Fe; th = 0:h:T; f = f1*[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 100 50 20 10 ... 35 5 15]; s = synthad(a,f,0*f,T,Fe); t = T*[0 .001 .01 .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 des composantes fréquentielles allure de la courbe d'enveloppe a(t) a 1 0.6 0 0.4 0.6 t

  7. Analyser la fonction synthad tirée du fichier synthad.m function s = synthad(a,f,p,T,Fe) % 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 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))

  8. a a(k+1) a(k) t t(k) t(k+1) 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 end n = length(th); env = zeros(1,n); ni = lt-1; c = zeros(1,ni+1); b = c; h2 = 0; for k = 1:ni h1 = h2+1; h2 = 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(h1:h2) = c*th(h1:h2)+b; end env = .99*env/max(env); 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 A\B 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)

  9. 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; f0 = 440; % la3 est la première note temp = 2.^((0:12)/12); % fr = f0*temp; % fréquence des notes de la3 à la4 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 1 et le joue

  10. Et si vous préférez la flute, comment utiliser le fichier ‘flute.m’? function s = flute(f1,T,Fe) % s = flute(f1,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*f1; 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 = f1/Fe*2; wn = (fn*[.8 1]).^(2); b = fir1(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)));

  11. 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; f0 = 440; // la3 est la première note temp = 2.^((0:12)/12); // fr = f0*temp; // fréquence des notes de la3 à la4 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 exec("envelop.sci"); // par exemple

  12. Utiliser le type ‘cell’ de Matlab, exemple de ‘creegammes.m’ %% créer des notes et des airs avec les cellules gamme={'do','dod','re','red','mi','fa','fad','sol','sold','la','lad','si'}; dt=power(2,1/12); for g=1:5 frla=110*2^(g-1); for n=1:length(gamme) eval([gamme{n},num2str(g),'=frla*dt^(n-10);']) end end %%créer une mélodie jouer(la3,.3,1); jouer(si3,.3,0.75); jouer(dod4,.3,0.5); jouer(mi4,.3,.5); jouer(re4,.3,.5); jouer(re4,.3,.5); jouer(fad4,.3,1); jouer(mi4,.3,1); jouer(mi4,.6,1); jouer(la4,.4,1); jouer(sold4,.4,1); jouer(la4,.8,1);

  13. 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= 440Hz, Dur= 1s, amp= 1, et Fs= 8000Hz. % 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=440Hz, Dur=1s, amp=1, Fs=8kHz 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

  14. %% paramétrage N=1024 Fe=44100 D=2; f1=880; %% son de cloche s=cloche(f1,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', num2str(N), ... ' points']) calculer et afficher le spectre avec Matlab

  15. Afficher le spectrogramme avec Matlab %% spectrogramme d'un arpège s=cloche(f1,D,Fe); s=[s,cloche(5*f1/8,D,Fe),... cloche(3*f1/4,D,Fe)]; wavplay(s,Fe) spectrogram(s,N,0,N,Fe,'yaxis') colorbar % cf. spectrogramme ci-contre  %% idem pour un accord acc=cloche(f1,D,Fe)+ ... cloche(5*f1/4,D,Fe)+ ... cloche(3*f1/2,D,Fe); acc= 0.99*acc/max(abs(acc)); wavplay(acc,Fe) spectrogram(acc,N,0,N,Fe,'yaxis') colorbar % non tracé

  16. Traduire le fichier creegammes.m en Scilab (creegammes.sce) //créer des notes et des airs avec les cellules gamme=['do','dod','re','red','mi','fa','fad','sol','sold','la','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 end //créer une mélodie exec('jouer.sce'); s=jouer(la3,.3,1); s=[s,jouer(si3,.3,0.75)]; s=[s,jouer(dod4,.3,0.5)]; s=[s,jouer(mi4,.3,.5)]; s=[s,jouer(re4,.3,.5)]; s=[s,jouer(re4,.3,.5)]; s=[s,jouer(fad4,.3,1)]; s=[s,jouer(mi4,.3,1)]; s=[s,jouer(mi4,.4,1)]; s=[s,jouer(la4,.4,1)]; s=[s,jouer(sold4,.4,1)]; s=[s,jouer(la4,.4,1)]; sound(s,8000); savewave('majoie.wav',s,8000)

  17. Traduction en Scilab des deux functions ‘synthad’ et ‘envelop’ function s = synthad(a,f,p,T,Fe) // 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 0.99] s = .99*s/max(abs(s)); endfunction function [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 end n = length(th); env = zeros(1,n); ni = lt-1; c = zeros(1,ni+1); b = c; h2 = 0; for k = 1:ni h1 = h2+1; h2 = 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(h1:h2) = c*th(h1:h2)+b; end env = .99*env/max(env); endfunction

  18. Traduire pour Scilab le script Matlab dans ‘cloche.m’ function s = cloche(f1,T,Fe) // s = cloche(f1,T,Fe) // imitation d'une cloche // f1 = fréquence fondamentale // Fe = fréquence d‘échantillonnage // T = durée du son //--------------------------------------------- h = 1/Fe; th = 0:h:T; f = f1*[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 100 50 20 10 ... 35 5 15]; s = synthad(a,f,0*f,T,Fe); t = T*[0 .001 .01 .4 .6 .9 1]; a = [0 .6 1 .4 .2 .1 0]; env = envelop(t,a,Fe); s = s.*env; endfunction

  19. 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=440Hz, si a n'est pas précisé, a=0.75, % dure=1 si non précisé, si fe n'est pas précisé, fe=44100Hz % 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 end% optionnel 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 ? note.m 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);

  20. 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 script (images comprises) F9 evaluate selection F5 run L’exécution du script ead.m depuis command window : >> eadseulement si le dossier contenant ead.m est dans le Path, ou dans le current directory File/Set Path, modifie le path : suivi de Add Folder puis Save, et Close En Scilab, faire exec(‘fichierScript’). 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 L’éditeur de scripts de MATLAB Dock Publish html Execute Cell

  21. 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(‘la3.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,s1,t,s2) 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(‘monFichierScript’) pour pouvoir appeler les fonctions définies dans monFichierScript la dernière ligne ‘endfunction’ est obligatoire 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 :

More Related