Calculs audio numérique pour les nuls

[alka]

Tu peux détailler un peu ton test avec ton vogelpick ?

Le principe est le suivant :
une IR convoluée avec sa symétrique garde la meme réponse en fréquence a une magnitude doublée après FFT mais sa phase devient linéaire.

c'est sans doute cette propriété qu'exploite Rephase et d'autres pour linéariser la phase. Si en plus l'impulse initiale est centrée sur 0, ça fait un filtre zéro phase.

On a vu un peu plus haut que la convolution d'IR dans le domaine temporel correspond a la multiplication des réponses fréquentielles.

mon intuition me disait que si y(t) et la symétrique de x(t) alors il y a bien des chances que Y (jw) soit le conjugué de X (jw)

Pour tester j'ai repris l'excel dont on a parlé hier :
pour Filtre1 un filtre HighPass BUTT12 a 200Hz (celui qui a servi qui a servi quelques pages plus haut)
et réponse Filtre 2 = conjugué réponse du Filtre 1
le tout appliqué à un HP bidon qui a une réponse en frequence bleue et une phase nulle.
Résultat = Filtre1 * Filtre2 * HPc (multiplication complexe)
résultat en rouge
On voit que le BUTT12 s'applique comme prévu et que la phase au lieu de bouger ne bouge pas du tout.

ca prouve pas que j'ai raison, mais ça prouve que je n'ai peut etre pas tort

Ca me plait de faire des filtres a phase linéraire en IIR

Edit : En regardant de plus près, l'effet sur l'amplitude du BUTT12 est doublé. c'est pas tout a fait ça mais pas loin
ca s'explique tres bien : (r,w) * (r,-w) donne (r^2, 0)

[thierry38]

Il y a la méthode du reverse IIR mais...

un IIR--->réponse impulsionelle infinie.
càd qu'il faut que le signal par exemple à 50Hz ait fait tout sa période dans le biquad------>et après le même biquad avec les coeff inversés.
ça se fait en tamponnant le signal (par exemple 100ms de mise mémoire) avant de le passer dans les biquads inversés,mais au final ,la latence est quasi la même que du FIR.

En gros ça veut dire qu'en temps réel c'est impossible/merdique.
ça se fait en offline,càd qu'il faut passer ses fichiers musicaux à la moulinette,et ensuite les écouter...
De plus pour faire des pentes de 100dB/oct ou autres spécialités à la rePhase,c'est imbitable avec des cascades de biquads.

Comme de toute façon les DSP/CPU modernes (et les PC) traitent des longs FIR sans soucis,c'est la solution la plus élégante/souple.

mon intuition me disait que si y(t) et la symétrique de x(t) alors il y a bien des chances que Y (jw) soit le conjugué de X (jw)

conjugué,je ne sais pas si c'est le terme,
Souvent appelé retournement temporel.(la phase est croissante quand les basses fréquences précèdent les hautes fréquence
s)

édité une nouvelle fois,une énorme faute de gammaire---->il faut que...est ,.....corrigé par ait.cette honte

[alka]

c'est marrant d'essayer
et puis pondre le biquad qui donne la réponse conjuguée d'un autre biquad c'est pas gagné.


en même temps je ne me faisais aucune illusion : c'est pas alka en 2016 qui va trouver comment faire de la phase linéaire en IIR
mais je fais bien semblant

faut vraiment se mettre au FIR donc... je n'y échapperai pas

[alka]

j'ai relu internet et j'avais mal lu le premier coup : une IR convoluée avec sa symétrique (sa version retournée temporellement si on veut) a bien la phase linéaire mais sa réponse en magnitude est mise au carré.

exactement ce qui se passe dans ma manip. mon vogeplick a raison même si je vois pas a quoi ça peut me servir

promis je ferai le ménage dans ce fil

[louisr]

j'ai relu internet et j'avais mal lu le premier coup : une IR convoluée avec sa symétrique (sa version retournée temporellement si on veut) a bien la phase linéaire mais sa réponse en magnitude est mise au carré.

exactement ce qui se passe dans ma manip. mon vogeplick a raison même si je vois pas a quoi ça peut me servir

Oui pour la phase et le module c'est assez intuitif, mais pour le conjugué (ou retournement temporel), bien vu !
j'essayerai de montrer ça mathématiquement dans le train tout à l'heure

[thierry38]

et puis pondre le biquad qui donne la réponse conjuguée d'un autre biquad c'est pas gagné.

Simplement inverser les coefficents

a0 devient a2,et vice versa.

a0 <----> a2
b0 <----> b2

et en normalisant,diviser tous les coeffs par a2 au lieu de a0.(avec la forme classique,les a au dénominateur (feedback))

ci-dessous,mais avec des b au dénominateur (http://www.dspguide.com/ch19/4.htm)



Il n'y a plus de z-1,z-2...mais z+1,z+2...(càd des samples en avance)
il suffit de réécrire l'équation en multipliant par z-2,(=retarder le tout de 2 samples)

y[n]=a0.x[n-2] + a1.x[n-1] + a2.x[n] + b0.y[n-2] + b1.y[n-1] + b2.y[n].

[alka]

tiens, j'avais lu cette page et pas percuté qu'on pouvait appliquer le reverse recursion a un biquad. bien vu le coup du * z^-2
ca me parait trop simple du coup. j'essaie

[alka]

alors c'est plutot incroyable : ça marche !

en vert la réponse normale d'un LR24 a 1kHz
(pour mémoire, un LR24 c'est un BUTT12 * BUTT12)

en rouge la réponse d'un BUTT12 * conjugué_BUTT12.
Même résultat que le LR24 pour l'amplitude (rouge et vert sont confondus) mais déphasage nul tout le long




Ca ouvre des perspectives cette affaire. Vais paufiner mes jouets excel pour voir jusqu'ou on peut aller avec ces "filtres IIR reversal convertis en biquad sans retard"

[alka]

puisque tout semble possible en ce moment: est ce possible de calculer les coeffs d'un biquad dont la réponse en amplitude est unitaire et la réponse en phase est inchangée ?

sais pas si je suis clair : reformulons...
On a les 5 coeffs du biquad B qui donne la réponse complexe H (A,w) avec A = gain et w = phase
je voudrais en déduire les coeffs du biquad B' qui a pour réponse H'(1,w)
possible ?

[alka]

test dans le dspiy d'un "BUTT12 conjugué" selon le principe précédent : marche pas.

grosse oscillation : le biquad ainsi créé n'est probablement pas stable.

Dommage, c'est pas de cette façon qu'on crééra a moindre frais un filtrage non déphasant en IIR.

La méthode "reverse recursion" fonctionne en mémorisant des samples pour appliquer l'algorithme en décalé, et ça je vois pas comment faire dans sigmastudio. Piste abandonnée.

[thierryvalk]

grosse oscillation : le biquad ainsi créé n'est probablement pas stable.

J'ai déjà lu qu'il existe une méthode pour voir si un biquad est stable ou pas.
La question: avec nos filtres et surtout le fait de mettre les biquads les uns derrière les autres, peut-on arriver à avoir un système instable ?

[louisr]

grosse oscillation : le biquad ainsi créé n'est probablement pas stable.

Une fonction de transfert est stable si toutes les racines de son polynôme caractéristique sont à partie réelles strictement négative. Tu dois pouvoir calculer ça sur ton biquad

[alka]

je dirai biquad stable si poles dans le cercle unitaire. Pas sur que les zeros doivent l'etre aussi. Est ce compatible avec ta formulation Louis ? humm... demande réflexion

je viens de faire le calcul et effectivement un des poles du biquad_conjuqé est juste en dehors du cercle unitaire. vraiment peu. module = 1,05. Peut etre qu'en trichant on peut faire quelquechose. reste a savoir tricher intelligemment

@thierry : nos filtres sont composés de filtres qui sont stables. Leur mise en cascade qui est une multiplication ne peut que rester stable. On ne créé pas de nouveau pole autre que ceux qui existent déjà.

[thierryvalk]

Je souhaiterais que l'on contrôle ce que je fais.
J'ai des fonctions, qui fonctionnent et le tout est de faire un ensemble efficace.

Pour commencer une fonction qui calcule amplitude et phase pour un biquad à fréquence donnée.
fresponse est une structure comportant magnitude et phase.
biquad est une structure de biquad tout simplement.
Il faudrait donc s'arrêter à un moment et retourner non plus fresponse mais H(z)

Code : Tout sélectionner


private fresponse getResponseAt(double f, biquad s)
        {
            double a1 = s.a1;
            double a2 = s.a2;
            double b0 = s.b0;
            double b1 = s.b1;
            double b2 = s.b2;
            fresponse frp = new fresponse();
            f = 2 * f / fs;

            double omega = Math.PI * f;
            double numeratorReal = b0 + b1 * Math.Cos(omega) + b2 * Math.Cos(2 * omega);
            double numeratorImag = -(b1 * Math.Sin(omega) + b2 * Math.Sin(2 * omega));
            double denominatorReal = 1 + a1 * Math.Cos(omega) + a2 * Math.Cos(2 * omega);
            double denominatorImag = -(a1 * Math.Sin(omega) + a2 * Math.Sin(2 * omega));
            double magnitude = Math.Sqrt((numeratorReal * numeratorReal + numeratorImag * numeratorImag)
                / (denominatorReal * denominatorReal + denominatorImag * denominatorImag));
            
            double phase = Math.Atan2(numeratorImag, numeratorReal) - Math.Atan2(denominatorImag, denominatorReal);
            if (phase >= Math.PI) { phase -= 2 * Math.PI; }
            else if (phase <= -Math.PI) { phase += 2 * Math.PI; }

            if (magnitude != 0) { magnitude = 20 * Math.Log(magnitude) / Math.Log(10); }

            frp.phase = phase * 180 / Math.PI;
            frp.magnitude = magnitude;

            return frp;
        }
 

[louisr]

je dirai biquad stable si poles dans le cercle unitaire. Pas sur que les zeros doivent l'etre aussi. Est ce compatible avec ta formulation Louis ? humm... demande réflexion

je viens de faire le calcul et effectivement un des poles du biquad_conjuqé est juste en dehors du cercle unitaire. vraiment peu. module = 1,05. Peut etre qu'en trichant on peut faire quelquechose. reste a savoir tricher intelligemment

Exact, j'ai écrit un peu vite : je parlais bien pôle de la fonction de transfert, que j'ai la mauvaise habitude de mal nommer. Et également ma relation est vraie pour une transformé de Laplace, càd un temps continu. Pour une transformé en Z c'est bien les pôles dans le cercle unitaires !



Thierry je n'ai pas trop compris ta question

[thierryvalk]

La question, c'est ou qu'il est le nombre complexe dans cette fonction qui va me servir pour faire ma courbe de x biquad + autre courbe ?

La fonction qui calcule l'ensemble des points en f pour un biquad en les rajoutant a la courbe en cours est celle-ci.
On y passe donc pour chaque biquad en série.
Cette fois elle retourne un tableau de fresponse.

Code : Tout sélectionner

private fresponse[] filtre_courbe(fresponse[] fc, biquad s)
        {
        double x = 0;
        int i = 0;

        fresponse frp = new fresponse();

        while (x < 20000)
        {
            i++;
            x = Math.Pow(2, (double)i / 48);

            frp = getResponseAt(x, s);

            if (i == 1) { fc[i].magnitude = frp.magnitude; fc[i].phase = frp.phase; }
            else
            {
                fc[i].magnitude = fc[i].magnitude + frp.magnitude;
                fc[i].phase = fc[i].phase + frp.phase;

                if (fc[i].phase >= 180) { fc[i].phase -= 2 * 180; }
                    else if (fc[i].phase <= -180) { fc[i].phase += 2 * 180; }
            }
        }
        return fc;
        }
 

[alka]

mais pourquoi tout est si compliqué et illisible.... ce que tu fais pour calculer les réels et imaginaires des numérateurs et dénominateurs me parait suspect.

tu as les opérations + * ^ complexes dans c# ?

En complexes, ca s'écrit en 8 lignes (et encore, c'est détaillé )

Code : Tout sélectionner

 
    w = 2 * PI * Freq / Fs   

    Zb = Complexe(Cos(w), -Sin(w))   
    Zb2 = ImPower(Zb, 2)
   
    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))


[alka]

si tu dois multiplier des biquads entre eux, tu ne calcule pas leur amplitude et phase, tu fais directement la multiplicaition des réponses complexes H.

et tu prend Amplitude Phase tout a la fin si tu veux tracer le résultat.

Pour l'addition, pareil.