Calculs audio numérique pour les nuls

[alka]

Suite a discussion sur le sujet Corrections Filtrages semi automatisé viewtopic.php?f=2&t=148, je me lance.

Ce sujet pourra par la suite servir à toutes les questions matheuses expliquées pour que tout le monde puisse comprendre

Introduction sur les nombres complexes

préliminaire : il faut que le langage de programmation supporte les nb imaginaires.
Il faut au moins pouvoir les déclarer et utiliser les opérateurs +, *, ^ et module et argument.
(j'appelle ça des imaginaires comme les anglo-saxons. C'est tout de suite moins compliqué que des complexes)


Quelques rappels sur les nombres imaginaires (rappels limités a ce qui est utile dans le contexte.)

Un nombre imaginaire s'écrit sous forme algrébrique C = a + ib
(ou a + jb quand on est physicien qui veut pas utiliser i qui représente le courant. Je noterai i.)

On peut se représenter géométriquement C comme un point de coordonnées (a,b) dans un plan.

C peut aussi s'écrire sous forme angulaire (r,φ) avec r = module et φ = argument , cad le rayon et l'angle du point C.
φ en radians

C'est la forme angulaire (appelée aussi polaire) qu'on préfère en audio et on va vite comprendre pourquoi.

On peut passer des coordonnées polaires (r,φ) aux coordonnées algébriques a+ib par les formules trigonométriques habituelles.



note: φ est dans l'intervalle ]Pi, Pi] cad ]-180°, +180°]. D'où l'emploi de atan2 pour son calcul.


Application à l'audio

En audio, on manipule fréquences, amplitudes et phase.

Par exemple, dans un fichier de mesure on voit pour chaque fréquence f l'amplitude G(dB) et la phase φ (degrés).
On lui fait correspondre la valeur imaginaire H de coordonnées polaires (A,θ) tel que :
A = 10^(G/20)
θ = Pi/180 * φ (pour passer de degrés en radians)

On peut repasser d'une valeur complexe H à son Amplitude et Phase par la conversion inverse :
Amplitude (dB) = 20 log (Mod(H)) c'est log10
Phase (degrés) = 180/Pi (Arg(H)) pour le passage de radians en degrés

Maintenant qu'on est rassuré sur le fait de pouvoir toujours convertir facilement de l'un a l'autre, on peut commencer.

[alka]

Fonctions de transfert complexe

Tout filtre a une fonction de transfert qu'on peut exprimer en imaginaire. Elle s'applique à un nb imaginaire et son résultat est imaginaire.
Les Fonctions de transfert complexes sont souvent notées H(z) dans la littérature.

Au lieu d'avoir deux fonctions pour un filtre T(f) et P(f) qui donnent respectivement l'amplitude et la phase résultante pour la fréquence f ,
on aura H(z) qui, pour l'imaginaire z, donne en résultat une valeur imaginaire dont on pourra tirer Amplitude et Phase.

Mais qui est z ?

Dans toutes les fonctions de transfert exprimées en complexe, on considère qu'elles s'appliquent a un complexe unitaire. Cad dont le module vaut 1. Géométriquement, il est placé sur le cercle de rayon 1 autour du point d'origine O.
On peut l'écrire en polaires z = (1,w) ou en algébrique z = cos(w)+ i sin(w)

Le w en question est défini avec la fonction de transfert.
Pour les fonctions de transfert de biquads par exemple, par définition, w = 2.Pi.f/Fs avec Fs la fréquence de sampling


Quelques propriétés des nombres imaginaires unitaires :

De façon générale le conjugué de C = a + ib c'est C barre = a - ib (noté avec un trait au dessus de C, une barre)
En polaires, le conjugué de C(r,w) est Cb(r,-w).

Dans le cas particulier d'un imaginaire unitaire, le conjugué est aussi son inverse
Pour z = cos(w) + i sin(w) le conjugué c'est zb = cos(w) - i sin(w) = 1/z = z^-1


Réponse complexe d'un biquad
Connaissant la fonction de transfert imaginaire d'un biquad H(z), on calcule amplitude et phase de la réponse du filtre à cette fréquence

pour chaque freq f on calcule :
w = 2 * Pi * f /Fs
zb = complexe [cos(w) - i sin(w)] (note : zb est l'inverse du classique z = cos(w) + i sin(w) )

Hz = (b0 + b1*zb + b2*zb^2) / (1 + a1*zb1 + a2*zb^2)

Hz est la réponse complexe, calculée avec les opérateurs d'addition, multiplication et de puissance complexes.

Si on veut l'Amplitude et Phase de H(z) à la fréquence f :
Amplitude (dB) = 20 log (Module(Hz))
Phase (degres) = 180/pi * Argument(Hz)

[alka]

Opérations sur les fonctions de transfert

Maintenant qu'on sait que tout peut etre exprimé en nb imaginaires, on va pas se gener.

La réponse mesurée en fréquence d'un HP formée de valeurs d' Amplitude(dB) et Phase(degrés) sera transformée, pour chaque fréquence, en valeur imaginaire Hm :
Hm= complex(r,w)
r = 10^(G/20)
w= pi/180*Phi pour avoir w en radians
ce qui s'écrit aussi en algébrique Hm = r*cos(w) + i r*sin(w)

Chaque filtre (low pass, high pass, peq, délai, etc...) a une fonction de transfert imaginaire aussi. Elle s'applique a un imaginaire et donne en résultat un imaginaire.

Et maintenant l'opération qui fait la magie et la simplicité en nb imaginaires :
La réponse de fonctions de transfert en série est la multiplication des fonctions de transfert
multiplication complexe bien sur.

Concrètement:
si je cascade un LP de fonction de transfert H1 avec trois PEQ de fonction de transfert H2, H3 et H4 et que j'applique le tout à la mesure d'un HP convertie en complexes Hm , l'opération pour trouver Résultat se résume à faire Resultat = Hm * H1 * H2 * H3 * H4

L'amplitude et phase résultante à une fréquence donnée se déduisent par l'opération qu'on connait bien sur le module et l'argument de Resultat pour la fréquence f :
Amplitude (dB) = 20 log (Mod(Resultat))
Phase (degrés) = 180/pi (Arg(Resultat))

Les autres opérations spéciales utiles :
- inverser la phase c'est faire * -1
- appliquer un gain G(dB) c'est multiplier par 10^(G/20)
- appliquer un délai de T (secondes) c'est multiplier par le nb imaginaire (cos(phi) - i sin(phi)) avec Phi = 2.Pi.f. T

edit: correction formule pour le sens du déphasage dû à un délai. merci thierry38.

[alka]

Les filtres classiques et Biquads

reste une chose, c'est quoi les fonctions de transfert des filtres classiques ?

on pense aux filtres calssiques : Low pass, high pass, (Butt, Bessel, LR) et aussi aux autres : PEQ, band pass, low shelve, high shelve, etc

La plupart de ces filtres ont été étudiés et on trouve un peu partout leur fonction de transfert imaginaire
wikipedia par exemple donne celle d'un low pass du 2eme ordre : en fonction de sa fréquence f0 et de son Q



la façon d'écrire H(jw) c'est pour dire que c'est la réponse complexe pour la pulsation w

Donc
soit on connait toutes les fonctions de transfert des filtres dont on a besoin et on les implémente,
soit, on passe par les biquads. C'est le choix que j'ai fait.

Pourquoi passer par les biquads ?
dabord parce que j'ai eu du mal a trouver les fonctions de transfert complexes de tous les filtres qui m'intéressent et aussi parce que je voulais pouvoir vérifier les coefficients (bénéfice collatéral)

Passer par les biquads comme étape intermédiaire est intéressant parce que :
- on peut exprimer n'importe quel filtre sous forme de biquad ou de cascade de biquads et
- la fonction de transfert d'un biquad est parfaitement connue.



le code vba correspondant :

Code : Tout sélectionner

    ' pour chaque fréquence Freq :
     
    w = 2 * PI * Freq / Fs   

    Zb = Complexe(Cos(w), -Sin(w))   ' zb est l'inverse de z
    Zb2 = ImPower(Zb, 2)          ' zb2 est le carré de zb
   
    ' H = (b0 + b1*zb + b2*zb2) / (1 + a1*zb + a2*zb2)
   
    Hnum = ImSum(b0, ImProduct(b1, Zb), ImProduct(b2, Zb2))   
    Hdenom = ImSum(1, ImProduct(a1, Zb), ImProduct(a2, Zb2))
   
    H = ImDiv(Hnum, Hdenom)

   Amplitude = 20 log(module(H))
   Phase = 180/Pi(Arg(H))

Suffit de savoir déterminer les coefficients des biquads.
Grace au succès des DSP ça a été fait pour toutes sortes de filtres classiques et moins classiques.

Le cookbook de Robert Bristow Jones donne les clés pour calculer les coefficients de biquads pour les filtres de base. La moitié du web s'est inspirée de ses travaux ! http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt

Ensuite, pour les filtres LP et HP spéciaux (Butt, Be, LR de différents ordres) on peut les faire par composition des filtres de base et en adaptant Fc et Q selon les tables calculées par nos ancêtres. http://www.ti.com/lit/an/sloa049b/sloa049b.pdf
par exemple :
- un BUTT12 est un filtre du second ordre avec Q=0,707. Nécessite un biquad.
- un BUTT24 c'est le chainage de deux BUTT12, le premier avec Q=0,5412 et le second avec Q=1,3017. Nécessite deux biquads.

Voici dans un pdf comment combiner les filtres de base du second ordre pour obtenir les filtres classiques (BUT, Bessel, LR) dans des ordres plus élevés. download/file.php?id=1534

En sachant faire des biquads, on sait tout faire ou presque

Voilà le code en VBA de mon excel pour le calcul des coeffs des filtres de base. Se lit comme du pseudo code :

Code : Tout sélectionner

' intermediate variables
w0 = 2 * PI * Fc / Fs

Select Case Filter
    Case "LPF6"
    ' Low Pass 1rst order
    '
    K = Tan(w0 / 2)
    alpha = 1 + K
    b0 = K / alpha
    b1 = b0
    b2 = 0
    a0 = 1
    a1 = -((1 - K) / alpha)
    a2 = 0
    
    Case "HPF6"
    ' High Pass 1rst order
    '
    K = Tan(w0 / 2)
    alpha = 1 + K
    b0 = 1 / alpha
    b1 = -b0
    b2 = 0
    a0 = 1
    a1 = -((1 - K) / alpha)
    a2 = 0
    

    Case "LPF12"
    ' Low Pass 2nd order
    ' LPF6: H(s) = 1 / (s^2 + s/Q + 1)
    
    alpha = Sin(w0) / (2 * Q)
    b0 = (1 - Cos(w0)) / 2
    b1 = 1 - Cos(w0)
    b2 = (1 - Cos(w0)) / 2
    a0 = 1 + alpha
    a1 = -2 * Cos(w0)
    a2 = 1 - alpha


    Case "HPF12"
    ' High Pass 2nd order
    ' HPF6: H(s) = s^2 / (s^2 + s/Q + 1)
    alpha = Sin(w0) / (2 * Q)
    b0 = (1 + Cos(w0)) / 2
    b1 = -(1 + Cos(w0))
    b2 = (1 + Cos(w0)) / 2
    a0 = 1 + alpha
    a1 = -2 * Cos(w0)
    a2 = 1 - alpha
    
    Case "PEQ"
    ' Peaking EQ
    '
    alpha = Sin(w0) / (2 * Q)
    A = Sqr(10 ^ (Gain / 20))
    
    b0 = 1 + alpha * A
    b1 = -2 * Cos(w0)
    b2 = 1 - alpha * A
    a0 = 1 + alpha / A
    a1 = -2 * Cos(w0)
    a2 = 1 - alpha / A
    
    Case "APF"
    ' AllPass filter.
    '
    alpha = Sin(w0) / (2 * Q)
    
    b0 = 1 - alpha
    b1 = -2 * Cos(w0)
    b2 = 1 + alpha
    a0 = 1 + alpha
    a1 = -2 * Cos(w0)
    a2 = 1 - alpha

    Case "LSH"
    ' LowShelve filter.
    '
    alpha = Sin(w0) / (2 * Q)
    A = Sqr(10 ^ (Gain / 20))
    
    b0 = A * ((A + 1) - (A - 1) * Cos(w0) + 2 * Sqr(A) * alpha)
    b1 = 2 * A * ((A - 1) - (A + 1) * Cos(w0))
    b2 = A * ((A + 1) - (A - 1) * Cos(w0) - 2 * Sqr(A) * alpha)
    a0 = (A + 1) + (A - 1) * Cos(w0) + 2 * Sqr(A) * alpha
    a1 = -2 * ((A - 1) + (A + 1) * Cos(w0))
    a2 = (A + 1) + (A - 1) * Cos(w0) - 2 * Sqr(A) * alpha
    
    Case "HSH"
    ' HighShelve filter.
    '
    alpha = Sin(w0) / (2 * Q)
    A = Sqr(10 ^ (Gain / 20))
    
    
    b0 = A * ((A + 1) + (A - 1) * Cos(w0) + 2 * Sqr(A) * alpha)
    b1 = -2 * A * ((A - 1) + (A + 1) * Cos(w0))
    b2 = A * ((A + 1) + (A - 1) * Cos(w0) - 2 * Sqr(A) * alpha)
    a0 = (A + 1) - (A - 1) * Cos(w0) + 2 * Sqr(A) * alpha
    a1 = 2 * ((A - 1) - (A + 1) * Cos(w0))
    a2 = (A + 1) - (A - 1) * Cos(w0) - 2 * Sqr(A) * alpha
    
    Case Else
    ' error
    MsgBox ("Erreur fonction BQcoeffs. Filtre " & Filter & " inconnu.")
    Exit Function
    
End Select

coeffs(1) = b0 / a0
coeffs(2) = b1 / a0
coeffs(3) = b2 / a0
coeffs(4) = a1 / a0
coeffs(5) = a2 / a0


[alka]

maintenant je suis fatigué et vais faire la sieste

retour ce soir pour les questions si je n'ai anesthésié l'auditoire et peut etre un exemple concret

[thierryvalk]

Super grand merci, j'attend avec impatience l'exemple concret commenté.

[thierry38]

Hello,

- appliquer un délai de t (secondes) c'est multiplier par le nb imaginaire (cos(phi)+ i sin(phi)) avec Phi = 2.Pi.f. t

J'aurai dit cos - i sin pour le retard (l'angle est décalé "dessous" dans le sens inverse des aiguilles d'une montre).
ou e(-s.t) dans Laplace.avec s=jw et t=1/Fs *n(pour le domaine discret/numérique)

le + étant les Diracs en avance.
D'ou Z^+1,Z^+2,Z^+xxx

[alka]

Retard ou avance est relatif.
sin(-w) = -sin(w)

dans ma logique, un délai c'est un retard, donc un déphasage plus grand. Un w positif doit faire grandir la phase, donc + dans la formule.
Si on veut appliquer un délai qui correspond a une avance, il faut appliquer un délai négatif.

La suite arrive, le temps de faire l'exemple concret.

[thierry38]

La convention (logique) est de considérer le - comme un retard.(tel un condensateur et son impédance -j/cw)

D'ou z-1,z-2,z-3... qui sont des échantillons retardés.(un retard Laplacien s'écrit e(-s.t)).Toujours avec s=jw et t=1/Fs
La valeur du coeff devant n'est pas représentée (enfin si,c'est =1 (Dirac) pour simplifier)

avec z=e(-s.t)----->cos (w.t) - j sin (wt).

Z+1,z+2,z+x sont des échantillons en avance,le FIR en phase maximum par exemple.

Ou plus simplement,une fonction quelconque y=f(t)
une représentation graphique de y=f(t-Tau) montrera un décalage sur la droite de Tau (évènement passé).

[alka]

ah ces matheux.... trop de gros mots en quelques phrases

raisonnons physiquement et simple : un micro et un HP
mesure 1 : admettons la phase mesurée au micro à 1kHz est de 45°
mesure 2 : on éloigne le HP de 1cm. Donc on ajoute un délai de 0,03ms qui est un retard pour le sens commun.
Ce délai a 1000Hz correspond a un déphasage de 10°
La nouvelle phase mesurée sera de 35° ou de 55° ?
Yaka essayer pour voir. Moi je dis 55°

[alka]

Exemple concret d'application des fonction de transfert complexe

voilà une illustration concrete des traitements. dans excel bien entendu

imaginons le HP qui a la réponse suivante : (pure invention)


on a la mesure dans un fichier qui contient Freq, Amplitude(dB) et Phase(degres) de REW par exemple

colonnes A à F concernent le HP
colonne A : Freq contient les frequences mesurées (pas beaucoup, c'est un exemple. c'est du 1/2 octave !)
colonne B : Amp(dB) amplitude G mesurée
colonne C : Phase(deg) phase P mesurée

je calcule :
colonne D : le gain A = 10^(G/20)
colonne E : la phase en radians Omega(rad)=Pi/180*P
j'en déduis la colonne F qui est la réponse du HP en nb complexe:
HPc(jw) = A*cos(omega) + i A*sin(omega)


On veut lui appliquer deux traitements pour l'exemple :
Filtre 1 HighPass BUTT12 à 200Hz
Filtre 2 PEQ à 1000Hz, Q=2 G=+10dB

En valeur intermédiaire, on calcule les coeffs des biquads de ces deux filtres (cases orange du tableau)

A chacun de ces biquads correspond une réponse complexe qui figure respectivement dans les colonnes G et H ( H1c et H2c)

La réponse finale est HRc = HPc*H1c*H2c dans la colonne I
Colonnes J et K contiennent respectivement RAmp et Rphase de la réponse finale calculées avec les formules habituelles:
RAmp = 20*log(mod(HRc))
RPhase = 180/Pi*Arg(HRc)

et voilà les courbes.
bleu = HP original. rouge avec les filtres appliqués



une partie des lignes de données est masquée pour ne pas encombrer

[alka]

l'excel

ExempleTV.xlsm

(26.86 Kio) Téléchargé 913 fois

ya la macro qui calcule la fonction de transfert d'un biquad dedans.

edit : oops, mauvaise version. je viens de remettre la bonne

[thierryvalk]

la réponse finale est HRc = HPc*H1c*H2c dans la colonne I

Pourquoi on multiplie ?
Actuellement on additionne amplitude et phase de chaque réponse en fréquence de chaque biquad.


sinon, pour résumer :

On a les impulses qui sont du temporel.
Via une FFT peut passer en fréquentiel.
En sortie de FFT on peut garder une table en imaginaire.
Si l’on désire appliquer des corrections IIR, on calcule pour chaque fréquence une fonction de transfert en imaginaire ce qui nous fera une nouvelle table.
On multiple chaque nombre imaginaire des 2 tables et l’on obtient une réponse en fréquence et phase imaginaire que l’on peut remettre en amplitude et phase.

Juste ?

[alka]

la réponse finale est HRc = HPc*H1c*H2c dans la colonne I

Pourquoi on multiplie ?
Actuellement on additionne amplitude et phase de chaque réponse en fréquence de chaque biquad.

on multiplie les valeurs complexes .
La réponse en amplitude c'est après avoir appliqué le log() et on se rappelle que log(a*b) = log(a) + log(b)
la multiplication complexe s'occupe de tout. Reste plus qu'a prendre amplitude et phase a la fin et ça roule.

démonstration:
si on a F1 et F2 de réponse complexe respectivement (A1,w1) et (A2,w2)
la mise en série des deux a pour réponse R la multiplication complexe des réponses. R=F1*F2 cad R (A1*A2, w1+w2)
Pour l'amplitude et phase de R :
GR (dB) = 20log(A1*A2) = 20log(A1) + 20log(A2) = GF1 + GF2
PhaseR(°) = 180/Pi*(w1+w2) = PhaseF1 + PhaseF2

Pour ton résumé : suis pas sur, faut que je re-regarde les fft.

ps: a la relecture, je réalise que je réemploie le mot complexe. Zut j'avais dit que j'utiliserai nombres imaginaires...

[thierryvalk]

Sans remettre en cause ces nombres qui permettent beaucoup de chose, vu d'un côté purement informatique c'est plus relatif.
Par exemple, tu as utilisé les biquads pour le calcul de filtres.
Un biquad c'est des coefficients, en informatique on peut les réunir dans un structure et donc on ne va plus manipuler des a1, b1... mais biquadX et biquadX.a1 pour accéder a un coefficient particulier.
On crée quelques fonctions de traitement de ces biquads et l'on en vient a oublie leurs lourdeurs.
Au final, un nombre complexe c'est une structure comprenant un nombre (partie) réel et un second imaginaire et l'on utilise bien des fonctions pour les calculs.

[alka]

disons que les langages de programmation bien faits optimisent les opérateurs complexes et ça ne ralentit pas vraiment.
D'un coté tu calcules des logs a répétition et de l'autre des multiplications et un seul log. Va savoir quel chemin prend le plus de temps.

J'utilise les biquads comme intermédiaire qui n'est en effet pas optimal; On pourrait prendre les fonctions de transfert de chaque type de filtre directement sans les convertir en biquads, ca enlèverait du temps de calcul et augmente la complexité du code.
Au final, c'est pas bien important. L'opération se fait devant l'utilisateur et a la demande. Si ça prend 200ms ou 300ms au lieu de 100 on s'en fiche un peu.

Pour ton résumé

On a les impulses qui sont du temporel. oui
Via une FFT peut passer en fréquentiel. oui
En sortie de FFT on peut garder une table en imaginaire. oui la sortie de la FFT c'est un tabelau de nombres complexes. tu gardes bien la moitié du tableau qui contient les freq positives.
Si on désire appliquer des corrections IIR, on calcule pour chaque fréquence une fonction de transfert en imaginaire ce qui nous fera une nouvelle table. je pense qu'on dit la meme chose mais pas sur pour le vocabulaire. Si H1, H2,... Hn sont les fonctions de transfert complexes des différents filtres, on fabrique la table de Resultat qui contient des complexes Resultat = FFT*H1*H2*...*Hn. On multiplie bien chaque valeur de la table, exactement comme je fais dans excel. On pourrait dire que la colonne HPc est la table FFT plutot que le résultat d'une mesure de HP.
On multiple chaque nombre imaginaire des 2 tables et l’on obtient une réponse en fréquence et phase imaginaire que l’on peut remettre en amplitude et phase.oui tu peux le faire en deux étapes si tu préfères: Dabord une table Corrections= H1*H2*...*Hn et ensuite Resultat = FFT * Corrections

[thierryvalk]

Je commence a sortir du brouillard et j'espère avoir compris.
Mais je pense bien que oui.
Merci pour ces explications structurées, car en général les explications sur les maths se font avec des maths pures et dures tout aussi indigestes et incompréhensibles que le problème de base.

[alka]

oui c'est comme en politique: Tu prend une question simple de bon sens et tu la soumets a un politique.
Après quelques minutes, tu comprend même plus la question que t'as posée


Filtre a phase linéaire en IIR ?

changeons de sujet : en sachant faire une FFT inverse et une convolution en c#, je pense bien qu'on pourrait tenter de rendre les filtres a phase 0 ou linéaire directement dans dstudio

FFT on sait faire

algorithme de invFFT c'est faisable a partir de l'algorithme de la FFT.
http://www.dsprelated.com/showarticle/800.php

la convolution peut se faire par définition a partir de FFT et invFFT.

FFT existe aussi en excel. Excel est vraiment lent et ça ne marcherait que pour des tables de 1024 ou moins. A essayer un jour de pluie
je met le lien pour ma mémoire
http://www.dsprelated.com/showthread/co ... 8342-1.php