Nucléo – Filtrer une bande fréquentielle plus efficacement

Comment filtrer une bande fréquentielle plus efficacement ?

NIVEAU 3

Objectifs

Mettre en oeuvre un filtre à réponse impulsionnelle infinie (IIR)
Caractériser un filtre à réponse impulsionnelle infinie (IIR)
Etudier la stabilité d’un filtre à réponse impulsionnelle infinie (IIR)

Pré-requis

Filtrer une bande fréquentielle / spectrale 3

Filtre IIR / Principe

Ce sont des filtres dits récursifs. La sortie est calculée en fonction d’échantillons d’entrée mais également des précédents échantillons de sortie, d’où le terme récursif. Ils ont une structure de ce type :

Contrairement à un filtre non-récursif (ou FIR), ces filtres dits récursifs calculent leurs échantillons de sortie à la fois en se basant sur une série de \(M\) échantillons d’entrée, mais en se basant également sur une série des \(N\) précédents échantillons de sortie.

La sortie est alors calculée de cette manière : \(s[k] = s_k = \sum_{i=0}^{M-1} \alpha_i \cdot e[k-i] + \sum_{p=1}^{N} \beta_p \cdot s[k-p]\)

L’ordre d’un tel filtre est donné par le maximum entre \(N\) et \(M\).

Il est a noter que ces filtres permettent généralement d’obtenir des ordres bien inférieurs aux filtres non-récursifs, pour un gabarit quasi-identique (i.e. des fréquences caractéristiques similaires pour une fréquence d’échantillonnage égale et des gains identiques dans les différentes zones utiles du gabarit). Ceci permet donc de gagner en temps de calcul. Cependant, ce type de filtre est plus difficile à mettre en oeuvre par leur rebouclage qui peut entrainer une instabilité.

Conception du filtre sous Matlab

Certaines étapes sont les mêmes que pour un filtre IIR.

Plusieurs logiciels existent pour concevoir des filtres numériques. Nous utiliserons Matlab, cette fois-ci sans sa boite à outils Filter Designer, car les filtres ainsi conçus n’ont pas la structure que permet d’implémenter les cartes Nucléo. En effet, ces dernières permettent l’implantation de structure de type Lattice.

Il est à noter que la démarche proposée ici se base sur les fréquences réduites \(\nu = f / F_e\) (où \(f\) est la fréquence considérée et \(F_e\) la fréquence d’échantillonnage souhaitée).

Nous allons voir ici une méthode permettant de définir le gabarit du filtre de manière non graphique et d’en afficher sa réponse en fréquence. Ce tutoriel est basé sur la page de Mathworks : IIR Filter Design.

Plusieurs types de filtres

Il existe 4 grands types de filtres IIR, dits classiques, concevables avec Matlab :

Butterworth : [b,a] = butter(n,Wn,options)
Chebyshev Type I : [b,a] = cheby1(n,Rp,Wn,options)
Chebyshev Type II : [b,a] = cheby2(n,Rs,Wn,options)
Elliptic : [b,a] = ellip(n,Rp,Rs,Wn,options)

Matlab propose alors 4 fonctions associées à la conception de ces filtres (voir ci-dessus). Chacune d’entre elles fournit les coefficients \(\beta\) (b) et \(\alpha\) (a) du filtre numérique associé. Par défaut les filtres réalisés sont de type passe-bas.
Dans tous les cas, il est nécessaire de préciser plusieurs paramètres :

n : l’ordre du filtre
Wn : un ensemble de fréquences caractéristiques (normalisées / réduites)
Rp : le taux d’ondulation autorisée dans la bande passante (en dB)
Rs : la valeur minimale d’atténuation dans la bande atténuée (en dB)

Remarque Il est à noter que l’on peut ajouter certaines options à l’ensemble de ces fonctions. Parmi ces options, il est possible de spéficier le type de filtre que l’on souhaite : ‘high’ pour un filtre passe-haut (à partir d’un modèle passe-bas) et ‘stop’ pour un filtre coupe-bande (à partir d’un modèle passe-bande).

Gabarit et conception

Pour pouvoir calculer la réponse en fréquence (par exemple) et les coefficients du filtre qui seront utilisés pour l’implémentation sur la carte Nucléo, il est donc indispensable de connaître l’ordre du filtre (noté \(n\)) et les fréquences caractéristiques de coupure (notées \(W_n\)).

Pour cela, Matlab propose d’autres fonctions qui permettent à partir du gabarit du filtre de récupérer ces informations.

Butterworth : [n,Wn] = buttord(Wp,Ws,Rp,Rs)
Chebyshev Type I : [n,Wn] = cheb1ord(Wp,Ws,Rp,Rs)
Chebyshev Type II : [n,Wn] = cheb2ord(Wp,Ws,Rp,Rs)
Elliptic : [n,Wn] = ellipord(Wp,Ws,Rp,Rs)

Filtre passe-bas

Par exemple, on pourra utiliser le code Matlab suivant pour générer un filtre passe-bas de fréquence d’échantillonnage de 10kHz, de bande passante 0-3.5kHz (variation de 1dB autorisée) et de bande atténuée 4-5kHz (à 60 dB en dessous du gain dans la bande passante) :

Fs = 10e3;
Fpass = 3.5e3;   Wpass = Fpass/Fs;
Fstop = 4e3;     Wstop = Fstop/Fs;
Apass = 1;
Astop = 60;

[n1, Wn1] = buttord(Wpass, Wstop, Apass, Astop)
[n2, Wn2] = cheb1ord(Wpass, Wstop, Apass, Astop) 
[n3, Wn3] = cheb2ord(Wpass, Wstop, Apass, Astop) 
[n4, Wn4] = ellipord(Wpass, Wstop, Apass, Astop)

Ce code fournit alors l’ordre \(n\) de chacun des types de filtres, ainsi que la fréquence caractéristique dans \(W_n\).

Filtre passe-bande

Pour concevoir un filtre passe-bande de fréquence d’échantillonnage de 10kHz, de bande-passante 2-3.5kHz (variation de 3 dB) et de bandes atténuées 0-1.5kHz et 4-5kHz (à 80 dB en dessous du gain dans la bande passante)

on pourra utiliser le code Matlab suivant pour générer un filtre passe-bas de fréquence d’échantillonnage de 10kHz, de bande passante 0-3.5kHz (variation de 1dB autorisée) et de bande atténuée 4-5kHz (à 80 dB en dessous du gain dans la bande passante) :

Fs = 10e3;
Fpass = [2e3 3.5e3];   Wpass = Fpass/Fs;
Fstop = [1.5e3 4e3];   Wstop = Fstop/Fs;
Apass = 1;
Astop = 80;

[n, Wn] = buttord(Wpass, Wstop, Apass, Astop);

Filtres passe-haut et coupe-bande

Pour la conception de filtres passe-haut et coupe-bande, ce sont les mêmes fonctions que précédemment. Seules les fréquences \(F_{pass}\) et \(F_{stop}\) sont inversées pour coller au gabarit.

C’est lors de la génération des coefficients que l’on choisit définitivement le type de filtre que l’on souhaite.

Génération des coefficients

Il est maintenant possible de générer les différents coefficients du filtre afin de pouvoir l’implémenter dans un calculateur et d’en obtenir la réponse en fréquence.

Pour cela, on va s’appuyer sur les fonctions suivantes :

Butterworth : [b,a] = butter(n,Wn,options)
Chebyshev Type I : [b,a] = cheby1(n,Rp,Wn,options)
Chebyshev Type II : [b,a] = cheby2(n,Rs,Wn,options)
Elliptic : [b,a] = ellip(n,Rp,Rs,Wn,options)

qui permettent de générer les coefficients \(\alpha\) et \(\beta\) du filtre à partir des données obtenues par les fonctions précédentes.

Par exemple, pour obtenir les coefficients d’un filtre elliptique, on pourra utiliser la fonction suivante :

[b,a] = ellip(n,Apass,Astop,Wn);

Les coefficients du filtre sont stockés dans les variables \(b\) et \(a\).

La structure implémentable sur une carte Nucléo est de type Lattice. Il reste donc une dernière étape avant de pouvoir récupérer les coefficients, les traduire en structure Lattice via la fonction tf2latc, qui renvoie les coefficients transformés.

[k v] = tf2latc(b,a);

Pour pouvoir par la suite les implémenter sur une carte Nucléo, il faut pouvoir exporter ces variables dans un format lisible. On va ici utiliser le format CSV (compatible avec Excel) et qui sépare les différentes valeurs par des virgules (,). On peut utiliser la fonction suivante :

csvwrite('k_coeff.csv',k')
csvwrite('v_coeff.csv',v')

Les coefficients obtenus s’appellent \(k\) et \(v\) pour une structure de type Lattice.

Réponse en fréquence et réponse impulsionnelle

A partir des coefficients précédents, on peut obtenir la réponse en fréquence à l’aide de la fonction suivante (où \(b\) et \(a[§latex] sont les coefficients du filtre et [latex]h\) la réponse en fréquence complexe aux points de pulsation réduite \(w\)) :

[h, w] = freqz(b, a);

ou la réponse impulsionnelle grâce à la fonction suivante (où \(b\) et \(a\) sont les coefficients du filtre et \(h\) la réponse impulsionnelle aux instants \(t\)) :

[h,t] = impz(b,a)

Dans les deux cas, il est possible de préciser le nombre de points sur lequel on souhaite réaliser ces opérations.

Ainsi pour obtenir la réponse en fréquence du filtre elliptique passe-bas précédent, on peut utiliser les instructions suivantes :

[h, w] = freqz(b, a, 2048);
plot(w/pi, 20*log10(abs(h)));

Le troisième argument de la fonction freqz permet de préciser le nombre de points sur lequel calculer la fonction de transfert. On obtient alors la figure suivante :

On peut faire de même avec la réponse impulsionnelle :

[h, t] = impz(b, a, 200);
plot(t, h);

On obtient alors la figure suivante (le troisième paramètre de la fonction impz correspond au nombre d’échantillons – à la période \(T_E = 1/F_E\)) :

Exemple d’un passe-haut

Pour un filtre passe-haut, de fréquence d’échantillonnage de 10kHz, de fréquence de bande-passante de 4kHz (avec 1dB de différence de gain maximale dans cette bande), de fréquence de bande atténuée de 3kHz (avec -80dB dans cette bande), on obtient :

Fs = 10e3;
Fpass = 4e3;   Wpass = Fpass/Fs;
Fstop = 3e3; Wstop = Fstop/Fs;
Apass = 1;
Astop = 80;

[n, Wn] = cheb2ord(Wpass, Wstop, Apass, Astop); figure; 
[b, a] = cheby2(n,Astop,Wn, 'high');
[h, w] = freqz(b, a); plot(w/pi, 20*log10(abs(h)));
[k v] = tf2latc(b,a); 
csvwrite('k_coeff.csv',k') 
csvwrite('v_coeff.csv',v')

Les coefficients obtenus s’appellent \(k\) et \(v\) pour une structure de type Lattice.

On obtient alors la réponse en fréquence suivante (pour un filtre de type Chebyshev 2) :

Implémentation sur Nucleo via MBED

Il est recommandé de reprendre le mode Suiveur, comme décrit dans le tutoriel Filtrer une bande fréquentielle / spectrale 3.

Code de base – Mode Suiveur

Une première étape à valider avant d’intégrer la partie filtrage numérique est la mise en oeuvre des deux modules de conversion : analogique-numérique (ADC) et numérique-analogique (DAC). Pour cela, il peut être intéressant de revenir sur les tutoriels Faire une action à intervalle régulier, Récupérer un signal analogique et Générer une tension analogique.

Le principe est celui représenté dans la figure suivante.

#include "mbed.h"
#define TE          0.00003
#define SAMPLE_RATE  1/TE

void convert(void);

/* E/S */
AnalogIn    mesure(A0);         // A0 on Arduino Header
AnalogOut   analog_out(PA_5);

/* Timer a repetition */
Ticker tik;

float in, out;

int main() {
    /* Association de la fonction convert à un timer */
    tik.attach(&convert, TE);
    
    while(1) {
    }
}

void convert(void){
    in = mesure.read();    // Lecture de l'entree analogique
    out = in;
    analog_out.write(out); // Ecriture de la sortie analogique
}

Dans cet exemple, le signal de sortie recopie le signal d’entrée. Cela permet de vérifier le bon fonctionnement de ces deux étages de conversion, avant même l’ajout d’un calcul intermédiaire.

La relation entre l’entrée et la sortie est alors du type : \(s[n] = e[n]\).

Dans le code précédent, on peut également noter que l’acquisition du signal d’entrée (nommé mesure) et de la restitution de ce même signal vers la sortie (signal nommé analog_out) se fait à intervalle régulier (ici 30 us – soit une fréquence d’échantillonnage de 33kHz) grâce à l’utilisation d’un module de gestion du temps par interruption, le Ticker tik.

Intégration des coefficients du filtre

Il est temps à présent d’intégrer l’élément central de ce traitement numérique, l’étape de filtrage numérique.

Pour cela, une bibliothèque, nommée mbed-dsp, est à votre disposition ici. Elle contient tous les objets en lien avec les filtres IIR et les méthodes permettant la mise en place de la structure matérielle associée à un filtre de type IIR.

Téléchargez cette bibliothèque dans un répertoire de votre ordinateur et importez-la dans votre projet (Cliquez droit sur votre projet, puis Import Library / From Import Wizard, puis Upload, sélectionnez le fichier zip contenant la bibliothèque, puis Import / Import).

Cette bibliothèque contient en particulier les objets et méthodes suivants :

arm_iir_lattice_instance_f32 : type permettant de construire un objet de type filtre IIR, dont les coefficients sont des flottants sur 32 bits
arm_iir_lattice_init_f32 : fonction permettant d’initialiser la structure matérielle d’un filtre IIR (de type Lattice), dont les coefficients sont des flottants sur 32 bits
arm_iir_lattice_f32 : fonction permettant de calculer la nouvelle sortie d’un filtre IIR, dont les coefficients sont des flottants sur 32 bits

Il faut, dans un premier temps, ajouter à votre projet précédent (contenant le code de base pour l’acquisition et la restitution des données analogiques) le tableau des coefficients et le nombre de coefficients, en dehors de toute fonction (voir paragraphe sur le suiveur).

Coefficients du filtre

Il faut ensuite créer les variables permettant de stocker les coefficients du filtre. Pour cela, il faut connaîtrer le nombre de coefficients. Pour cela, il faut regarder la taille des variables contenant les coefficients sous Matlab.

On voit cela dans l’espace de travail (ou Workspace) de Matlab.

Dans le cas du filtre de Chebyshev 2 de type passe-haut, on se retrouve avec des variables de type double et de taille 13 (pour \(a\) et \(b\)). C’est la même chose pour les coefficients Lattice \(k\) et \(v\) :

Remarques On peut remarquer que l’un des tableaux est plus petit d’une case que l’autre.

Il va donc falloir créer deux tableaux de float32_t (équivalent à des nombres réels simple précision) de taille spécifiée précédemment, dans votre projet MBED.

#define NB_POINTS 13
float32_t k_coeff[NB_POINTS-1];
float32_t v_coeff[NB_POINTS];

Remarque Les coefficients des filtres sont théoriquement constants au cours du temps. Cependant, les fonctions utilisées par la suite pour calculer les échantillons ont besoin de variables et non de constantes pour travailler, même si à terme ces tableaux ne seront pas modifiés.

Il faut ensuite remplir ces tableaux avec les valeurs des coefficients. Pour cela, il faut ouvrir les fichiers *.csv exportés précédemment sous Matlab avec un éditeur de texte.

Il faut ensuite copier l’ensemble des coefficients et les coller dans le tableau déclaré précédemment sous MBED.

float32_t k_coeff[NB_POINTS-1]; = {-0.7775,0.98328,-0.84677,0.91293,-0.85179,0.8205,-0.73222,0.5711,-0.32475,0.11381,-0.021763,0.0017756};  
float32_t v_coeff[NB_POINTS]; = {-0.00062406,-0.00071921,-0.00067741,0.012544,0.038442,0.068748,0.08427,0.074132,0.048655,0.025367,0.010842,0.0035784,0.00075118};

Structure de stockage

Il faut ensuite initialiser, en dehors de toute fonction également, les différentes structures qui accueilleront les résultats intermédiaires des calculs. Pour cela, il faut utiliser la suite d’instruction suivante :

float32_t iir_state[NB_POINTS+1];
arm_iir_lattice_instance_f32 my_iir;

La première ligne permet de créer un tableau où seront stockées les derniers échantillons acquis (nécessaires pour le calcul de la sortie du filtre). La taille de ce tableau doit être la taille du tableau \(v\) incrémentée de 1. La seconde ligne permet de créer un objet de type filtre IIR, dont les coefficients sont des flottants sur 32 bits.

Par la suite, il faut faire appel à la fonction qui va initialiser la structure matérielle, basée sur le tableau précédent et les coefficients du filtre. Pour cela, il faut faire appel à la fonction suivante, une seule fois dans le main (par exemple) :

arm_iir_lattice_init_f32(&my_iir, NB_POINTS, k_coeff, v_coeff, iir_state, 1);

Les paramètres à lui fournir sont : l’objet de type filtre IIR (ici my_iir), le nombre de coefficients (ici NB_POINTS), le tableau contenant les coefficients (ici k_coeff et v_coeff), le tableau permettant de stocker tous les échantillons nécessaires (ici iir_state) et le nombre total d’entrées à utiliser (ici une seule entrée, car un seul filtre calculé).

Enfin, il faut faire appel régulièrement à la fonction qui exécutera le traitement numérique des données par l’intermédiaire de la fonction suivante :

arm_iir_lattice_f32(&my_iir, &in, &out, 1);

Les paramètres à lui fournir sont : l’objet de type filtre IIR (ici my_iir), le nouvel échantillon (ici in), la valeur de sortie (ici out) et le nombre total d’entrées à utiliser (ici une seule entrée, car un seul filtre calculé).

Exemple complet d’un filtre passe-haut

Voici un exemple complet d’un filtre FIR de type passe-haut (voir gabarit dans la section Exemple d’un passe-haut de cet article).

Filtre PASSE-HAUT – Cheby2 – FS = 10kHz, Fp = 4kHz, Fstop = 3.5kHz

#include "mbed.h"
#include "dsp.h"

/* E/S */
AnalogIn    mesure(A0);         // A0 on Arduino Header
DigitalOut  clk_test(PA_9);     // D8
AnalogOut   sortie(PA_5);   // D13

#define TE          0.0001

// Filtre PASSE-HAUT - Cheby2 - FS = 10kHz, Fp = 4kHz, Fstop = 3.5kHz
#define NB_POINTS 13
float32_t k_coeff[NB_POINTS-1] = {-0.7775,0.98328,-0.84677,0.91293,-0.85179,
                                    0.8205,-0.73222,0.5711,-0.32475,
                                    0.11381,-0.021763,0.0017756};  
float32_t v_coeff[NB_POINTS] = {-0.00062406,-0.00071921,-0.00067741,0.012544,
                                    0.038442,0.068748,0.08427,0.074132,0.048655,
                                    0.025367,0.010842,0.0035784,0.00075118};

float32_t iir_state[NB_POINTS+1];
arm_iir_lattice_instance_f32 my_iir;

/* Conversion routine */
void convert(void);

/* Timer a repetition */
Ticker tik;

float in, out;

int main() {
    arm_iir_lattice_init_f32(&my_iir, NB_POINTS, k_coeff, v_coeff, iir_state, 1);
    
    tik.attach(&convert, TE);
    
    while(1) {
    }
}

/*
 * ROUTINE DE CONVERSION ET CALCUL
 */
void convert(void){
    
    in = mesure.read()-0.5;   
    clk_test = 1;  
    arm_iir_lattice_f32(&my_iir, &in, &out, 1);
    clk_test = 0; 
    sortie.write(out+0.5);
    
}

MInE Prototyper Prototyper avec Nucleo et MBED

Nucléo – Filtrer une bande fréquentielle plus efficacement