Comment fonctionne un filtre passe-bas par programmation?

Je travaille sur un simple filtre passe-bas pour une mesure <100 Hz dans mon application. Mais jusqu'à présent, je lutte avec la théorie derrière tout cela. C'est cool que je l'ai fait fonctionner, mais je l'apprécierais vraiment si je savais comment / pourquoi cela fonctionne.

J'ai trouvé le code suivant:

void getLPCoefficientsButterworth2Pole(const int samplerate, const double cutoff, double* const ax, double* const by)
{
    double PI = M_PI;
    double sqrt2 = sqrt(2);

    double QcRaw  = (2 * PI * cutoff) / samplerate; // Find cutoff frequency in [0..PI]
    double QcWarp = tan(QcRaw); // Warp cutoff frequency

    double gain = 1 / ( 1 + sqrt2 / QcWarp + 2 / ( QcWarp * QcWarp ) );

    by[2] = ( 1 - sqrt2 / QcWarp + 2 / ( QcWarp * QcWarp ) ) * gain;
    by[1] = ( 2 - 2 * 2 / ( QcWarp * QcWarp ) ) * gain;
    by[0] = 1;

    ax[0] = 1 * gain;
    ax[1] = 2 * gain;
    ax[2] = 1 * gain;
}

Pour calculer les coefficients. Ensuite, dans les échantillons audio, je les «passe-bas» de cette façon:

        xv[2] = xv[1];
        xv[1] = xv[0];

        xv[0] = pData[j];
        yv[2] = yv[1];
        yv[1] = yv[0];

        yv[0] = (ax[0] * xv[0] + ax[1] * xv[1] + ax[2] * xv[2]
                   - by[1] * yv[0]
                   - by[2] * yv[1]);

        pData[j] = yv[0];

Pour obtenir une conception passe-bas.

Je me pose quelques questions:

1. Je reçois les échantillons audio dans un simple tableau flottant *. Quel est ce nombre flottant? La seule chose que je vois est un nombre, comment est-ce un morceau de son?
2. Le code utilise des calculs antérieurs (trois d'entre eux) dans le nouveau calcul par échantillon. Est-ce à dire que les 2 premiers échantillons de données ne sont pas filtrés correctement? (ce n'est pas important car ce ne sont que 2 échantillons, mais je me demande simplement)
3. En essayant de tout apprendre, j'ai trouvé quelques formules pour le filtre Butterworth (2e pôle). Comment ces formules sont-elles reflétées dans ce code? Aucune des formules que j'ai trouvées n'a ces calculs que vous pouvez voir dans la fonction 'getLPCoefficientsButterworth2Pole ()'.

1. Le nombre float * tableau est un pointeur sur le tableau. Il s'agit d'un nombre unique qui contient l'adresse du premier élément du tableau de valeurs flottantes.

2. Habituellement, la condition initiale (c'est-à-dire les éléments «passés» initiaux de x et y) est 0, mais si leurs valeurs ne sont pas égales à 0 ce n'est pas un gros problème non plus, car après un certain temps, les conditions initiales n'ont aucun effet sur la sortie signal pour tout filtre stable. Et votre filtre est évidemment stable.

3. Une fonction de transfert passe-bas du deuxième ordre avec caractéristique de Butterworth et fréquence de coupure $\omega_a=2*\pi*f_a$ (avec $f_a$ en Hertz) est donné par

$$H(s) = \frac{\omega_a^2}{s^2+\sqrt{2}\omega_a s + \omega_a^2}\tag{1}$$ Il s'agit d'un résultat standard que vous pouvez facilement trouver sur le Web. Afin d'obtenir un filtre à temps discret, vous pouvez appliquer une soi-disant transformation bilinéaire:

$$s = 2f_s\frac{z-1}{z+1}\tag{2}$$

où $f_s$est la fréquence d'échantillonnage. Cela est nécessaire car l'axe de fréquence des signaux analogiques ($0\le f\le\infty$) doit être mappé sur la plage de fréquences autorisée pour les signaux à temps discret ($0\le f\le f_s/2$). Comme cette transformation déforme les fréquences, nous devons calculer la fréquence de coupure analogique souhaitée à partir de la fréquence de coupure donnée$f_d$ dans le domaine du temps discret:

$$\omega_a = 2f_s \tan(\frac{\omega_d}{2})\text{ with }\omega_d = 2\pi f_d/f_s$$

Si vous insérez (2) dans (1) vous obtenez

$$H(z) = k\frac{z^2+b_1z+b_2}{z^2+a_1z+a_2}\tag{3}$$

avec

$$k=\frac{\alpha^2}{1+\sqrt{2}\alpha+\alpha^2}\text{ (gain)}$$ $$a_1 = \frac{2(\alpha^2-1)}{1+\sqrt{2}\alpha+\alpha^2},\quad a_2=\frac{1-\sqrt{2}\alpha+\alpha^2}{1+\sqrt{2}\alpha+\alpha^2},\quad b_1=2,\quad b_2=1$$

où j'ai utilisé $\alpha=\tan(\frac{\omega_d}{2})$. La fonction de transfert (3) correspond à l'équation du domaine temporel

$$y(n) = kx(n) + kb_1x(n-1) + kb_2x(n-2) - a_1y(n-1) - a_2y(n-2)$$

Comme vous pouvez le voir, il y a une certaine ressemblance avec votre code. Cependant, il existe également des différences et vous devez vérifier la réponse en fréquence de votre filtre. Je pense que les équations ci-dessus sont correctes, mais c'est à vous de les vérifier et de décider quelle version est la bonne.

Vous avez demandé comment fonctionne un filtre passe-bas et avez mentionné que le filtre utilise les valeurs passées de vos données. Il s'agit d'une discussion non technique de ce qui se passe dans un filtre passe-bas.

Le filtre passe-bas prend des vues différentes (décalées dans le temps) de votre signal, les met à l'échelle et les ajoute ensemble. Vous pouvez imaginer dessiner votre signal 3 fois, l'un étant courant, le second étant décalé d'un temps d'échantillonnage, le troisième étant décalé de 2 temps d'échantillonnage.

Aux basses fréquences, toutes les vues sont très similaires (le décalage d'un seul échantillon change à peine où vous êtes sur le signal à tout instant). Dans ce cas, les trois versions s'ajouteront de manière constructive (ou du moins non destructive), de sorte que le signal passe à travers le filtre.

Passant maintenant à des fréquences plus élevées, chaque version décalée du signal devient plus distincte à un instant donné (point d'échantillonnage) et peut même en fait s'inverser en signe. À ces fréquences plus élevées, les trois versions de votre signal ont tendance à s'annuler (ajoutées de manière destructive) afin que le signal soit atténué.

Différents types de filtres organisent ce type d'interférence constructive / destructrice sur les bandes de fréquences appropriées pour créer des filtres passe-bas, passe-bande ou passe-haut.

Cela dépend de la façon dont vous le souhaitez. Pour moi, j'ai mis en œuvre en C. Donc, ce que j'ai fait, c'est que j'ai généré des coefficients de filtre en matlabutilisant firls. Ensuite, j'ai stocké ces coefficients de filtre dans un tableau, puis je passe ces coefficients de filtre à la fonction de filtre dans C. Il matlabest assez facile de générer les coefficients de filtre, puis d'utiliser l'opération de convolution ou d'utiliser un filtre spécifique. Mais Cvous devez également implémenter le filtre. Pour ma part, je dois le faire pour les DSP, donc implémenté dans Cet obtenu les résultats en matlabutilisant la fonction mex. Pour calculer les coefficients, utilisez cette commande dans matlab:

n=17                  %filter order
f=[0 0.167 0.333 1]   %Frequency band edges
a=[1 1 0 0]           %Desired amplitudes
fir= firls (n, f,a )