Sauver un signal audio numérique dans un fichier wave

1. Sauver un signal audio numérique dans un fichier wave On trouve dans le fichier ‘la3.wav’ les trois blocs de données : RIFF, fmt et data. Les données sont disposées comme suit : on lit l’octet de poids faible d’abord, puis les autres en allant vers l’octet de poids fort (format little endian ou ‘petit-boutien’ ) • 52 49 46 46, signifie RIFF • 24 53 07 00, lire 00 07 53 24, soit 7*65536+83*256+36=480036 bytes • 66 6D 74 signifie fmt en ASCII • format PCM (code 01 00) • monophonie (01 00) une voie • 40 1F 00 00, soit fe=8 kHz • 40 1F 00 00 est le byte rate en octets par seconde • 01 00 indique qu’un échantillon est codé sur un octet • 08 00, lire 00 08, un échantillon occupe 8 bits • L’en tête du fichier la3.wav précise donc : 8kHz, mono, 16 bits/échantillon, format PCM, taille 480 044 octets, bit rate 8000 octets par seconde, ... • Comment doubler fe ?

2. Quantification d’un signal discret avec Scilab signal discret : x=[x(nTe)=xn, n=0.. N-1] xn codé sur B bits, intervalle -1=< xn < 1 le pas de quantification Q= 2/2B taille de x : N*B en bits Signal binaire : xbinn= partieEntière(xn/Q) -2B-1 =< xbinn < 2B-1, par exemple B=8, -128=<xbinn<128 Signal quantifié (ou numérique) xquantn=xbinn*Q, avec -1 =< xquantn < 1 Erreur de quantification : en = xn – xquantn Rapport signal sur bruit ou SNR, s’exprime en dB :SNR= 20*log10(écartType(x)/écartType(e))Exemple : SNR =72dB sur ligne téléphonique grand public Page 2

3. Composition fréquentielle ou spectre et TFD Lire f0= 440 Hz a0=0.75 (~0.8) fe = N =32 NTe = Df =250Hz spectre/N <0.4 M=256 pts tracés Lire : f0= a0= fe = N = NTe = Df = spectre/N =

4. Sous échantillonnage d’un signal sinusoïdal pur Sous échantillonnage M=8 • ssech=s(1:8:length(s)); • N devient N/8= 1024 • Te devient 8*Te • fe devient fe/8 =1000Hz • a0/2 = 512/1024, a0=1 • f0 = 440Hz • taux compression : 8

5. Énoncés de la contrainte de Shannon Énoncés de la contrainte de Shannon • D’après ‘A Mathematical Theory of Communication’, july 1948, in Bell System Technical Journal, par Claude Elwood Shannon (1916-2001) • Pour échantillonner correctement un signal s(t), il faut respecter la contrainte de Shannon : • Contrainte de Shannon simplifiée : la fréquence d'échantillonnage doit être égale au moins au double de la fréquence maximale du spectre du signal : • s'il existe fmax telle que S(f >fmax)=0, • alors fe >=2*fmax • Contrainte de Shannon générale : la fréquence d'échantillonnage doit être égale au moins au double de la largeur du spectre du signal : • s'il existe fmin et fmax telle que S(f >0) =0 pour f<fmin et pour f>fmax, • alors fe >=2*(fmax-fmin) • Si la contrainte n’est pas respectée, les échantillons ne permettent pas de reconstituer le signal s(t) ! • Conséquence : seuls les signaux ‘à bande limitée’ (c’est-à-dire dont le spectre est nul au-dessus d’une fréquence fmax) peuvent être échantillonnés correctement, d’où le filtre dit ‘anti-aliasing’ des cartes sons qui limite le spectre du signal à l’intervalle [0, fe /2[ avant l’échantillonnage Filtre reconstructeur et formule de Shannon • Voici la formule de Shannon qui reconstruit s(t) à partir des échantillons s(nTe) (seulement si la contrainte de Shannon est respectée !) : • En terme de filtrage, la formule de Shannon applique au signal discret le filtre reconstructeur de Shannon, pour retrouver s(t) : • ce filtre multiplie par Te les composantes fréquentielles entre –fe/2 et fe/2, et multiplie par 0 toutes les autres composantes du spectre pour supprimer • voici la réponse fréquentielle de ce filtre tracée entre -fe et fe : • Si la contrainte de Shannon n‘est pas respectée pour l’échantillonnage, la formule est impuissante, s(t) est perdu, les échantillons sont inutiles ! Te f 0 -fe/2 fe/2

6. 8 0 440 Sur-échantillonnage d’un signal discret dans un rapport M signifie insertion de M-1 échantillons nuls entre deux échantillons du signal // avec Scilab ou Matlab sse=zeros(size(s)); sse(1:M:length(s))=ssech; • N/8 redevient N échantillons • 8*Te redevient Te • fe/8 redevient fe =8000Hz • a0/2 reste 512/1024, a0=1 • f0 = 440Hz • Effets du sur-échantillonnage sur le spectre : • S(f) ne change pas, car on ajoute des échantillons nuls • fe étant multipliée par 8,on voit M=8 périodes de S(f) entre 0 et fe Comment retrouver le spectre du signal de départ?

7. spectre (R échantillons) R/4 3*R/16 R/8 f 0 fe f 0 fe f 0 fe H1(f) 4 0 f 0 fe Compresser et décompresser le signal de spectre suivant dans un facteur M=4

8. Créer et appliquer un filtre de réponse fréquentielle donnée avec Scilab // filtre passe bande 1000Hz-2000Hz // gain 4, R=64, fe=8000Hz R=64; fe=8000; n=0:R-1; fr=n*fe/R; H=4*[zeros(1,R/8),ones(1,1+R/8), ... zeros(1,-1+R/2),ones(1,1+R/8), ... zeros(1,-1+R/8)]; plot2d3(fr,H) xgrid xtitle(['H2,avec R=',string(R)], ... 'fréquence (Hz)’, ‘H’) //calcul des coefficients du filtre h=fftshift(real(ifft(H))); plot2d3(n/fe,h) xtitle('coefficients du filtre',... 'temps (s)',... 'h=fftshift(real(ifft(H)))') xgrid(); // filtrage [y,fe]=wavread('piano.wav'); disp(fe) // fe=8000 sound(y,fe) yf= convol(h,y); wavwrite(yf,fe,'pianofilt.wav') sound(yf,fe) //Spectrogrammes (Goldwave) // énergie Ey=(y*y')/2 // énergie y = 163.96 Eyf=(yf*yf')/2 // énergie yf =89.62 La réponse fréquentielle du filtre est définie dans le vecteur H, les coefficients du filtre sont calculés dans le vecteur h, on filtre ‘piano.wav’, on compare spectrogrammes et énergies avant et après filtrage Définition de l’énergie du signal x de taille X échantillons

9. Utiliser un banc de filtres pour analyser automatiquement le spectre d’un signal audio Lit le signal audio dans e et fe function [s, E, Esignal, fe]=bancfiltres(M, R, fichier, play) //fichier ‘bancfiltres.sce’ //utilisation [s,e,es,fe]=bancfiltres(8,128,'piano.wav',0); [e,fe]=wavread(fichier); N= R/(4*M); H=[ones(1,N-1),0.9,0.5,0.1,zeros(1,R-2*N-3),... 0.1,0.5,0.9,ones(1,N-2)]; h=fftshift(real(ifft(H))); n=0:R-1; for j=0:M-1 bande(j+1,:)=2*cos((2*j+1)*n*%pi/(2*M)).*h; end for j=0:M-1 sfiltre=convol(e,bande(j+1,:)); s(j+1,:)=sfiltre(1:length(e)); wavwrite(s(j+1,:),fe,['s'+string(j+1)+'.wav']); end Esignal= e*e'/2; E=diag(s*s')/2; disp(['sum(E):',string(sum(E))]) bar([0:M-1]*fe/(2*M),100*E/Esignal) xtitle(['Analyse de ',fichier],'frequence (Hz)'... ,'energie (% energie totale)') xgrid(); if play then sound(sum(s,1),fe); end endfunction Crée filtre générateur h Crée banc de M filtres De taille R partir du filtre h Filtre e dans sfiltre Ramène longueur sfiltre à celle de y et sauve dans des fichiers wave Calcule et affiche un diagramme barre des énergies en % de l’énergie de e Joue sum(s,1) si play!=0 Définition de l’énergie d’un signal x de taille X échantillons

10. Structure de CODEC* utilisant un banc de filtres *CODEC : coder decoder, cf. vocoder, Chicago, 1939 La structure du CODEC inclut cinq étages : Banc de filtres de réponses fréquentielles B1,B2,… BM et de longueur R avec : (B1+B2+… BM)*X=XN*M échantillons, C=1/M Sous-échantillonnage de rapport M, autorisé parce que la largeur de spectre des signaux x1, x2, … xM vaut fe/M (et donc la condition de Shannon générale est vérifiée)M*N/M échantillons, C=1 Etage de compression: objectif, atteindre C > 1, mais la compression est destructive, elle modifie le signal compressé Sur-échantillonnage (intercale M-1 échantillons nuls) Banc de filtres interpolateurs, obtenu en multipliant par M les réponses fréquentielles des filtres B1, B2, … BM Synthèse additive, le signal décompressé est noté xrec Algorithme compression Étage sur-é- chantillonneur Étage sous-é- chantillonneur banc de M filtres banc de M filtres Étage de compression xse1 x1 xd1 B1 xd2 xse2 x2 B2 x xrec … … xM xdM xM xseM BM

11. Compresser avec un banc de filtres en diminuant la longueur du codage binaire des échantillons Y a t’il des bits inutiles dans les signaux suivants codés sur 8 bits issus d’un banc de filtres ? Quel est le taux de compression qui en résulte ? -0.8<s1<0.6 xmax= m= bits ? u = ? -0.25<s2<0.2 xmax= m = bits ? u= ? -0.10<s2<0.10 xmax= m= bits ? u= ? -0.03<s2<0.04 xmax= m= bits ? u= ? C= ?

12. Luminance et chrominance • L’œil est plus sensible à la luminance (noir et blanc, clair et sombre) qu’à la chrominance (couleur), cf. ci-dessous tiré de http://semsci.u-strasbg.fr/efflum.htm • D’où le codage YUV (ou Y Cb Cr) des images TV à partir des informations R, G (vert) et B Y est la luminance, maximum de lisibilité pour l’œil humain, traduction de l’image couleur en niveaux de gris CB et CR sont les chrominances bleue et rouge, pour reconstruire les informations de couleur

13. y Im(x,y) H s(t) t D y découpé en M lignes x 0 L x découpé enN colonnes En résumé, 1/ la durée devient la dimension, longueur ou largeur, donc le temps devient espace, et 2/ on passe à deux dimensions. Du signal audio à … l’image numérique

14. Détection des contours par un filtre d’image // filtrer l'image lena.jpg // à l’aide d’un filtre laplacien lena=imread('lena.png'); h = fspecial('laplacian'); imf = imfilter(lena,h); imshow(imf);

15. Utilisation de la DCT pour compresser une image (cet exemple est tiré de l’aide de Scilab : - - > help dct) la DCT 2D est appliquée à l’image A, dont la composition fréquentielle (fréquences spatiales, verticale et horizontale) est calculée. Les composantes d’amplitude inférieure à 1 sont négligées, La ligne ‘size(find(d<>0),’*’)’ trouve 165 composantes non nulles sur 1680 Le taux de compression vaut un peu plus de 10, C= 1680/165 L’image décompressée à droite est obtenue en appliquant la dct inverse, soit A1=dct(d,1);

16. Schéma de principe de la compression JPEG (tiré de Wikipédia) Division de l’image en blocs de 8x8 pixels appelés ‘macroblocs’ Séparation de chaque bloc en plans Y (luminance), Cr et Cb (chrominances rouge et bleu), ces deux plans étant sous échantillonnés d’un rapport 2 suivant la hauteur et suivant la largeur, Transformation DCT (Discrete Cosine Transform) de chaque bloc : on obtient 8x8 coefficients de Fourier qui définissent la composition fréquentielle du bloc Quantification des coefficients : les plus faibles en valeur absolue sont annulés ou codés sur un nombre de bits plus faible (le pas de quantification est augmenté) : compression avec pertes. Compression des coefficients restants : codage RLE (Run Length Encoding), codage de Huffman ou VLC (Variable Length Coding)

17. Utilisation de la DCT dans le principe de compression JPEG L’image à compresser est découpée en blocs de 8x8 pixels, auxquels on applique la DCT, qui calcule 8x8 coefficients pour chaque bloc selon la formule suivante : N=8, pixel(x,y), avec x=0..7, et y=0.. 7 est un bloc de 64 pixels DCT(i,j), i=0..7, j=0..7 est le tableau des 64 coefficients DCT du bloc C(i) vaut 1 pour i non nul, et sqrt(2) pour i = 0 (de même pour C(j) et j) Le tableau DCT(i,j), i=0..7, j=0..7 contient le spectre du bloc de pixels, les fréquences spatiales normalisées varient entre 0 et fx/2=0.5 horizontalement et verticalement. 64 intensités donnent 64 coefficients DCT, taux de compression : C=1 La DCT inverse reconstitue le bloc de pixels à partir des coefficients DCT(i,j) i=0..7, j=0..7

18. La DCT décompose chaque bloc de 8x8 pixels en une somme pondérée des 64 images élémentaires ci-dessous : Origine des fréquences j=7 j=1 j=0 j i=0 i=1 DCT(i,j) forts DCT(i,j) faibles i=7 i Tiré du cours de Pierre Nerzic cité page 1 La DCT calcule les 64 coefficients de pondération de la somme des images élémentaires, ou encore la composition fréquentielle spatiale. La DCT inverse reconstitue le bloc de pixels en faisant la somme pondérée des images élémentaires. Une image élémentaire contient une fréquence horizontale et une fréquence verticale. Les coefficients de pondération sont les DCT(i,j), i=0..7, j=0 .. 7. La correspondance entre les indices i et j et la fréquence spatiale normalisée est : fi = i/(2*N), et fj= j/(2*N). • le bloc élémentaire associé au coefficient DCT(i,j), i=0..7, j=0..7, avec la fréquence horizontale j/16, et la fréquence verticale i/16 • Les DCT(i,j) décroissent quand i et j augmentent dans la plupart des cas (les composantes des blocs de pixels ont des fréquences spatiales basses)

19. Application à l’image ‘cameraman.png’ pour [Q,C]=quantzone(sim, 6), C= 3.047 [Q,C]=quantzone(sim, 6) : on conserve les 21 premiers coefficients DCT des fréquences spatiales les plus bassesbasses, soit i+j<6, On annule les autres

20. Représenter une fréquence spatiale i(x)= 0.5+ 0.5*cos(2*%pi*f*x) x varie de 0 à L Définition : N pixels entre 0 et L Période échantillonnage : L/N Résolution horizontale: fx = N/L Pixellisation : x= k*L/N, k= 0 … N-1 ik= 0.5*(1+cos(2*%pi*f*k*L/N) Normalisation de fe : L/N= 1 fx=1 pixel par unité de longueur  x= 0 .. N-1 ik=0.5*(1+cos(2*%pi*f*x)

21. Comparaison de DCT et de FFT s= [sn, n=0 .. N-1], de taille N, avec sn=s(n/fe) S=fft(s)=[Sk, k=0..N-1], de taille N, avec Sk=s(kfe/N) D= dct(s)= [Dm, m= 0 .. N-1], avec Dm=D(fm) S et D ont la même taille que s, soit N L’exponentielle complexe devient un cosinus, D est un vecteur réel, à la différence des Sk qui sont complexes Les fréquences des valeurs calculées sont différentes Pour fft, Sk=S(fk) Pour dct, Dm=D(fm) DCT est donc une variante de FFT, qui calcule un spectre réel et non complexe pour des fréquences entre 0 et fe/2 et non pas 0 et fe