Readme du TP | Sujet | Corrigé |
---|---|---|
Python | sujet python | corrigé python |
Octave | sujet octave | corrigé octave |
On désire appliquer ce filtre pour effecteur le même effet d'overtone que la chanteuse Ana Maria Hefele. Regardez l'explication sur
Nous avons besoin de faire le lien entre les fréquences et les notes (regardez le piano de la vidéo et l'échellen des fréquences). Les Anglais ont remplacé notre DO-RE-MI par les lettres de l'alphabet A-B-C en commençant avec A pour le LA ! Ce serait trop simple autrement...
... | DO | RE | MI | FA | SOL | LA | SI | DO | ... |
---|---|---|---|---|---|---|---|---|---|
... | C | D | E | F | G | A | B | C | ... |
En solfège, une octave correspond à une miltiplication par 2 de la fréquence et au fait de passer d'un DO au DO suivant par exemple. Donc tous les multiples de fréquence en puissance de 2 sonneront en harmonie et donneront la même note mais à des octaves différentes.
Les fréquences multiples qui ne sont pas en puissance de deux apparaissent naturellement (périodicité et séries de Fourier) et les premièrs harmoniques qui ne sont pas des puissances de deux donnent des notes qui sonnent différemment de la note de départ.
Les musiciens occidentaux ont convergés vers "la gammme tempérée" qui décompose une octave (multiplication par deux) en 12 demi-tons qui correspondent ainsi chacun à une multiplication par :
demi-ton $\quad \leftrightarrow \quad \times\; 2^{\frac{1}{12}}\approx \times\; 1.0595 \equiv +5.95 \% \approx +\frac{3}{12}\text{dB}= + 0.25$ dB.
On obtient ainsi une échelle de fréquence logarithmique qui décompose en 12 demi-tons le passage à l'octave (doublement de fréquence). Un barreau de référrence de cette échelle est le LA de la 4ᵉ octave noté "LA 4" officiellement à 440 Hz.
clear all; close all; clc;
LA4 = 440;
%% VOTRE CODE demiton= ..., demiton_dB... et quart_de_ton
% multiplier par demiton doit augmenter d'un demi-ton
% ^ permet de calculer une puissance non entière
demiton = 1.09
demiton_dB = 0.25
demiton_log2 = 1/12;
% multiplier deux fois par quart_de_ton
% doit faire un demiton...
quart_de_ton = demiton ;
%% VOTRE CODE LA4diese=... etc.
% LA + 1demiton = LA#
LA4diese = LA4 % * quoi ? donne 466.16 Hz
% LA + 2demiton = SI
SI4 = LA4 % *quoi ? donne 493.88 Hz
% Et oui!
%SI +1demiton = DO (de l'octave au dessus)
DO5 = LA4 % donne 523.25 Hz
% Descendez le DO5 de 2 octaves
DO3 = DO5 % donne 130.81 Hz
DO7 = DO5 % donne 2093 Hz
demiton = 1.0900 demiton_dB = 0.25000 LA4diese = 440 SI4 = 440 DO5 = 440 DO3 = 440 DO7 = 440
On rappelle que +3dB correspondent environ à une octave soit $\times\;2$ (car $10 \log{3}\approx 2$) et que 10dB exactement à x10.
20dB sont x100 et non pas x10 ne vous trompez pas !
Lorsque l'on dit qu'un gain x10 fait 20dB on sous-entend +20dB de ratio de puissance. Comme x100 en puissance font x10 en amplitude puisque $P=U^2/R$ l'embrouille est vite venue...
Attention
log
est en base 2 avec octave/matlab : il n'y pas de "ln" !
Préférez utiliser explicitementlog10
et ̀log2
.
De manière indépendante du langage, on peut toujours choisir sa base en faisantlog(x)/log(B)
où B est la base choisie.
Ici les dB expriment un ratio entre les fréquences des différentes notes. Comme l'échelle de fréquence des notes est logarithlique, augmenter d'un demi-ton sera toujours une augmentation de fréquence de +0.25dB.
Ne pas confondre avec les dB exprimant le ratio de puissance acoustique. Il exprime le rapport de puissance acoustique (on parle d'intensité sonore) avec une référence qui est l'intensité à la limite du seuil audible :
$0\text{dB} \leftrightarrow 1\text{pW}/m^2$ d'intensité acoustique
$\leftrightarrow 20 \mu P $ de pression acoustique
$\leftrightarrow$ "Bruissement de la brise dans les champs de blé"
Comme l'oreille humaine est plus ou moins sensible à différentes fréquences (très sensible autour de 1KHz) un facteur de correction peut être apporté en fonction de la fréquence. C'est donc bien un filtre qui est apporté au spectre d'intensité pour se rapprocher de la notion "d'intensité ressentie". Le standard donne 3 filtres de A à C et l'unité de mesure associée est communément notée dBA, dBB, dBC.
Lorsque l'on chante (un DO 4 par exemple) le signal de la voix est quasi-périodique à la fréquence DO4 et comporte naturellement des multiples (Séries de Fourier oblige !). On distingue :
La fondamentale n'est pas toujours un DO et on peut vouloir décaler les tons d'une chanson en partant d'une autre note. On parle plutôt en degrés où la note de degré I est la fondamentale et les autres degrés sont relatifs à cette note. On répartit ainsi 7 degrés de I à VII dans une octave.
On peut monter la gamme majeure de DO en partant du DO et en faisant des écarts de 2 demi-tons et parfois un seul demi-ton ! :
Degrés | I | II | III | IV | V | VI | VII | I | ... | |||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Demitons | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12=0 | 13=1 |
Gamme DO maj | DO | RE | MI | FA | SOL | LA | SI | DO | ... | |||||
Gamme FA maj | FA | SOL | LA | LA# | DO | RE | MI | FA | ... |
La gamme est dite majeure car, elle respecte les intervales "ton-ton-demiton-ton-ton-ton-demiton" entre chaques degrés. Cela est culturel et d'autres gammes existent bien sûr.
La gamme tempérée permet presque de représenter les harmoniques x 3 (quinte) et x 5 (tierce) fidèlement car :
Ainsi un instrument étant accordé pour jouer une gamme de DO avec la quinte juste sonnera mieux que celui accordé en gamme tempérée. L'intérêt apparait lorsque l'on veut jouer dans une autre gamme :
On veut générer un vecteur de fréquences allant du DO3 (vers 131 Hz) au DO7 (vers 2093 Hz) tous les quarts de ton. On doit obtenir un tableau de 97 valeurs du genre [131, ..., 2033, 2093] environs.
Les faibles utiliseront la fonction
logspace
(mais c'est + compliqué en fait !)
Les fortes feront une boucle et multiplieront par quart_de_ton
Les fabuleuses utiliseront la fonction log2, et le calcul vectorisé !
Choisissez une manière de faire et stockez le résulta tdans f pour vérification
%% Nombre de points
M = 4 * 12 * 2 +1 ; % 4 octaves * 12 demi-ton * 2 quart de tons
%% VOTRE CODE f_aible = ...
%% Le faible utilise logspace
f_aible = logspace(1,4, M);
%% OU BIEN VOTRE CODE f_orte=...
%% La forte utilise une boucle For
f=DO3;
for id=1:M
f_orte(id) = f;
end
%% OU BIEN VOTRE CODE f_abuleuse = ...
%% La geeke utilise la fonction 'ln'
% Soit on part des log2(f)
log2_f = 0 : 1/12/2 : 4; % *2^0 à *2^4 par quart de ton
% Soit carrément du vecteurs des indices "1:M"
f_abuleuse = DO3 *1:M;
%% CHOISISSEZ f= f_aible OU f=f_orte OU f=f_abuleuse
f=f_aible;
%% Ce code vérifie si c'est bon
if abs(f(1)-DO3)>1e-6
disp("La première fréquence n'est pas un DO3")
elseif abs(f(end)-DO7)>1e-3
disp("La dernière fréquence n'est pas un DO7")
elseif abs(f(end)/f(end-1)-2^(1/24))>1e-6
disp("Ce n'est pas espacé au quart de ton")
else
disp("Bravo, f semble bon");
end
M=length(f)
La première fréquence n'est pas un DO3 M = 97
On récupère le fichier son ana_ah.wav
où un "aaaaah" est chanté sur le ton d'un DO4. On utilise la fonction audioread
A vous de créer un vecteur temps de la bonne lognueur ayant la même longueur que le signal s
.
Normalisez le signal pour qu'il ait une norme infinie de 1 : $\|s\|_\infty=1$ (on peut utiliser la fonction max
)
racine = pwd; racine = [racine(1:(findstr(racine,"Signal")+6)), "discret/tp/tp2_overtone"];cd (racine);
fich="ana_ah";
[Y,Fe]=audioread(["./",fich,".wav"]);
s=Y(:,1); % récupère que le son gauche (colonne 1)
N=length(s);Te=1/Fe;
%% VOTRE CODE t=...
t=1:N;
%% VOTRE CODE s = s ...
%% Normalisant s en norme infinie
s=s/2;
%% VOTRE AFFICHAGE de s en fonction du temps
plot(s)
xlabel("temps [s]")
%% On peut créer un fichier .wav avec audiowrite ainsi
audiowrite(["./",fich,"_mono.wav"],s,Fe);
Inutile de bourriner une fft
alors que l'on ne veut voir que des fréquences "mélodiques" !
On utilise le vecteur de fréquences en échelle log précédant allant du D03 au DO7 tous les quarts de tons.
On construit la matrice de projection en échantillonnant $e^{i2\pi\,f\,t}$ en discret pour ces fréquences (transposée conjuguée) comme dans VEC2 bases fréquentielles et dans le TD FREQ sur la TFD en matriciel exo1 TFD corrige.
Ouvrez une console File -> New Console for Notebook, et regardez si f et t sont
des vecteurs lignes ou colonnes. Testez différentes produitsf*t
f'*t
etc. et regardez comment former la matrice de calcul de fft uniquement pour les fréquences voulues
printf(" N=%d points temporels et M=%d points fréquentiels\n",N,M);
printf("Construction de la matrice de calcul de TFD (%d x %d) pour les f en échelle log...\n",M,N);
tic;
%% VOTRE CODE W=...
% Signaux en lignes de la forme ci-dessous
% utiliser le produit matriciel d'un colonne par une ligne !
% ----------> t
% |
% | W = matrice conjuguée de exp(i.2.pi.f.t)
% V
W = t'*f;
toc
printf("\nCalcul de la TFD : soit %d x %d = %d Mflo (Millions d'opérations flottantes)...\n",length(f),N,N*length(f)/1e6);
tic;
%% VOTRE CODE tfd_s = ...
% Utilisant la matrice W et le signal s pour
% obtenir les "projections" aux fréquences voulues da,s tfd_s
tfd_s = W;
toc
clear W % Plus besoin ! et occupe 27M floats * 4 octets = 108 M de RAM !
N=274692 points temporels et M=97 points fréquentiels Construction de la matrice de calcul de TFD (97 x 274692) pour les f en échelle log... Elapsed time is 0.097903 seconds. Calcul de la TFD : soit 97 x 274692 = 26.6451 Mflo (Millions d'opérations flottantes)... Elapsed time is 0.00632691 seconds.
On peut maintenant afficher en échelle fréquentielle log et gain en dB le module : le spectre.
On ne s'intéresse pas à la phase puisque ce n'est pas un filtre mais un signal.
On utilise le script afficher_grille_notes
qui est plus haut dans l'arborescence du dépôt Signal/discret/utiles
. On regarde en échelle de fréquence log, puis linéaire.
addpath("../../utiles");
dBde = @(gain) max(-100,20*log10(abs(gain))) ;
subplot(211) % échelle log/log comme bode
ax=semilogx(f,dBde(tfd_s),'.-');
title("quelques octaves en echelle de frequence log...");
xlabel("Frequences [Hz] echelle log")
hold on;
grid_notes(); % Fonction maison pour les octaves
subplot(212)% échelle lin
plot(f,dBde(tfd_s),'.-');
hold on;
grid_notes();
title(".. en echelle lineaire")
xlabel("Frequences [Hz] echelle lin")
Utilisons l'algorithme fft
pour afficher ces fréquences.
Attention N points temporel donnent N points en fréquentiel répartis dans $[0, F_e[$.
À vous de créer un vecteur de fréquences linéaire f_lin
adéquat permettant d'afficher de même manière le module de la fft.
## VOTRE CODE Df=..., f_lin=...
df=Fe/N ; % La résolution fréquentielle en N points
f_lin = 0:df:(Fe-df); % Les fréquences linéaires
## VOTRE CODE fft_s= ....
# help fft si nécessaire
# Pourquoi fft n'as pas besoin de connaitre Fe ?!
fft_s = fft(s);
# On affiche de même le résultat mais avec l'algo fft
# Avec vos fréquences
pas_tout = round(110/df):round(9000/df);
subplot(211)
ax=semilogx(f_lin(pas_tout),dBde(fft_s(pas_tout)),'-');
title("quelques octaves en echelle de frequence log...");
xlabel("Frequences [Hz] echelle log")
hold on;
grid_notes();
subplot(212)
plot(f_lin,dBde(fft_s),'-');
hold on;
grid_notes();
title(".. en echelle lineaire")
xlabel("Frequences [Hz] echelle lin")
Remarquez alors que le spectre est symétrique par rapport à $\frac{F_e}{2}$ !
Pas étonnant car (sélectionnez la bonne expliquation) :
On peut utiliser fftshift
pour afficher une période du spectre centrée en 0.
## VOTRE affichage avec fftshift centré autour de 0 en fréquence
f_sym = f_lin - Fe/2;
subplot(211);
plot(f_sym, fftshift(dBde(fft_s)));
subplot(212);
pas_tout = round(N/2-1900/df):round(N/2+2000/df);
plot(f_sym(pas_tout),fftshift(dBde(fft_s))(pas_tout));
On remarque que la tierce majeure donne un MI d'amplitude plus faible que les DO en harmonie et le SOL de la quinte.
On peut choisir d'amplifier soit :
On peut d'abord ajuster le filtre continu du second ordre du tp précédent pour qu'il amplifie la tierce (4ᵉ harmonique, donc celle de rang 5 à la fréquence $5.F_0$) avec +/- un demi ton de sélectivité ($\zeta$ <0.1).
On pourra accentuer l'effet en appliquant 2 à 3 fois ce filtre !
Que devient la fonction de transfert et la réponse harmonique en faisnt cela ?
Gc= @(p,wn,zeta) wn^2./(p.^2 + 2*zeta*wn*p + wn^2);
p = i*2*pi*f_lin';
pas_tout = round(100/df):round(3000/df);
plot(f_lin(pas_tout),dBde(fft_s(pas_tout))); hold on;
%% On recherche le pic de fréquence autour de F0
[maxi, id_max] = max(abs(fft_s)); % la fondamentale est maximum
w0 = (id_max*df-df) *2*pi;
puis = [1 1 1 1 1 1 1 1 1 1];
zetas = [0.2 0.1 0.08 0.07 0.06 0.04 0.03 0.02 0.02 0.02];
for rang = 1:10
wn = rang * w0; % on multiplie par 3 pour l'harmonique
plot(f_lin(pas_tout),dBde(Gc(p(pas_tout),wn,zetas(rang)).^puis(rang)));
end
rang=7;
wn = w0*rang;
zeta=zetas(rang);
pow=puis(rang);
Choisissez de continuer en notebook ou en scripts. Faites des sauvegardes dans tous les cas.
Vous pouvez continuer dans ce notebook le travail en utilisant File → New Console For Notebook Voir les variables en tapant "whos" dans la console
Il peut être utile de faire une copie de ce Notebook en cas de mise à jour
File → Duplicate et rename
On peut créer un point de suavegarde sans changer de nom
File → Snapshot
Vous pouvez exporter ce notebook en scripts
.m
et continuer de travailler avec octave ou matlab.
Pour cela dans un terminal lancez la commande
./genere_les_scripts.sh
cd scripts_octave
octave --gui
Vous pouvez travailler avec le GUI d'octave ou Matlab, avec emacs en mode octave (ALT+X octave-mode), ou un éditeur évolué (codium, pycharm)...
Utilisez la synthèse bilinéaire pour fabriquer un filtre discret équivalent au filtre continu que vous avez ajusté.
N'oubliez pas l'effet de compression des fréquences !
Donnez la fonction de transfert
Affichez la réponse harmonique par-dessus le spectre de la voix
Appliquez hors-ligne ce filtre en multipliant puis ifft
Multipliez le spectre du signal par celui de votre second ordre pour filtrer.
Faites une transformée inverse et généré un fichier son du signal filtré.
Vous mettrez ainsi au point votre fonction de transfert continue pour obtenir le meilleur effet sonore.
fft_filt_c = fft_s .* Gc(p,wn,zeta).^pow;
y_filt_c = real(ifft(fft_filt_c));
y_filt_c = y_filt_c/max(abs(y_filt_c));
audiowrite("filt_c.wav",y_filt_c,Fe)
Gbi = @(z) 2/Te*(z-1)./(z+1);
wnpre = 2/Te*tan(wn*Te/2);
z=exp(Te*p);
Gzpre = Gc(Gbi(z),wnpre,zeta);
Gz = Gc(Gbi(z),wn,zeta);
plot(f_lin(pas_tout),dBde(fft_s(pas_tout))); hold on;
plot(f_lin(pas_tout),dBde(Gc(p(pas_tout),wn,zeta).^pow),'k');
plot(f_lin(pas_tout),dBde(Gzpre(pas_tout).^pow),'g-');
plot(f_lin(pas_tout),dBde(Gz(pas_tout).^pow),'g--');
fft_filt_z = fft_s .* Gz.^2;
y_filt_z = real(ifft(fft_filt_z));
y_filt_z = y_filt_z/max(abs(y_filt_z));
audiowrite("filt_z.wav",y_filt_z,Fe)
On peut calculer comme si l'on était en temps réel la réponse du filtre en faisant une fonction récursive.
La boucle suivante "simule" un temps réel où à chaque nouvel échantillon, le filtre est calculé avec ce nouveau sample et donne la valeur de la sortie.
function [a,b] = filt_overtone(w0,rang,Te)
wn=2/Te*tan(rang*w0*Te/2);
puis = [1 1 1 1 1 1 1 1 1 1]*1;
zetas = [0.2 0.1 0.08 0.07 0.06 0.04 0.03 0.02 0.02 0.02]*1;
zeta=zetas(rang);
a0 = (1+4*zeta/(Te*wn)+4/(Te*wn)**2);
a1 = (2-8/(Te*wn)**2);
a2 = (1-4*zeta/(Te*wn)+4/(Te*wn)**2);
aa = [a0, a1, a2]/a0;
bb = [1, 2, 1]/a0;
a=aa;
b=bb;
for k=2:puis(rang)
a=conv(a,aa);
b=conv(b,bb);
end
end
%% Fonction pouvant implémenter un filtre IIR de type I
% Implique la mémorisation des x et des y
% cette mémoire est accessible en général par variables globales
global mem_x
global mem_y
filtre = @filt_iir_I; % version stupide de la fonction
%% CODEZ VOTRE FONCTION dans le fichier filtre_iir_I.m
mem_x = zeros(1,15); % mémoire du filtre vide au début
mem_y = mem_x;
%% Le vecteur des sorties à la même taille que s
rang=3;
inc=1;
y = 0*s;
a=1;
b=1;
for k = 1:N
if (mod(k,Fe/4)==Fe/4-1)
[a,b] = filt_overtone(w0,rang,Te);
rang = rang + inc ;
if (rang>9)
inc = -inc;
end
if (rang<3)
inc = -inc;
end
end
y(k) = filtre(s(k),b,a);
end
plot(t,y)
y = y /max(abs(y));
audiowrite("filt_rec.wav",y,Fe)
!./genere_les_scripts.sh
on sauve les fichiers .m existants de ./scripts_octave ls: impossible d'accéder à '*.m': Aucun fichier ou dossier de ce type On génère les .m à partir des .ipynb dans ./scripts_octave [NbConvertApp] Converting notebook tp2_overtone_corr.ipynb to script [NbConvertApp] Writing 6721 bytes to scripts_octave/tp2_overtone_corr.m Pour voir les sauvegardes ls -a ./scripts_octave Pour nettoyer les sauvegardes faites rm ./scripts_octave/.sauv.*
!aplay ana_ah.wav
Lecture WAVE 'ana_ah.wav' : Signed 16 bit Little Endian, Fréquence 48000 Hz, Stéréo