Comment supprimer une bande fréquentielle ?

NIVEAU 3

Objectifs

  • Mettre en oeuvre un filtre à réponse impulsionnelle finie (FIR)
  • Caractériser un filtre à réponse impulsionnelle finie (FIR)

Pré-requis

Filtre numérique et gabarit

Avant toute chose, il faut définir le gabarit du filtre, c’est à dire ses différentes spécificités : fréquences caractéristiques, atténuation autorisée…

Voici les 4 gabarits associés aux 4 grands types de filtres possibles :

Dans le cas de filtres numériques, on parle souvent de fréquence réduite \(\nu\). Il s’agit de la fréquence considérée (\(f\)) ramenée à la fréquence d’échantillonnage (\(F_e\) ou \(F_s\)). Ainsi : \(\nu = f / F_e\). Elle est souvent appelée \(W\) sous Matlab.

  • \(G_0\) : gain (en dB) dans la bande passante (souvent 0 dB)
  • \(F_S\) ou \(F_e\) : fréquence d’échantillonnage (en Hz)
  • \(W = f / F_S\) : fréquence réduite (ou normalisée)
  • \(R_{pass}\) : variation maximale de gain possible dans la bande passante (en dB)
  • \(A_{stop}\) : valeur maximale de gain (par rapport au gain dans la bande passante) autorisée dans la bande atténuée (en dB)
  • \(W_{px}\) : fréquences caractéristiques de la bande passante
  • \(W_{sx}\) : fréquences caractéristiques de la bande atténuée

Filtre FIR / Principe

Ce sont des filtres dits non-récursifs. La sortie est calculée uniquement avec une série d’échantillons d’entrée. Ils ont une structure de ce type :

Chacun des nouveaux échantillons de sortie est calculée en fonction de \(M\) échantillons d’entrée, auxquels on associe des coefficients permettant de donner la réponse impulsionnelle souhaitée au filtre (i.e. sa réponse en fréquence… tout est qu’une question de référentiel de travail – temporel ou fréquentiel).

Ainsi, les échantillons de sortie sont régis par l’équation suivante : \(s[k] = s_k = \sum_{i=0}^{M-1} \alpha_i \cdot e[k-i]\)

L’ordre d’un tel filtre est \(M\).

La transformée de Fourier d’un tel filtre est alors donnée par la relation suivante : \(H(\nu) = \sum_{p=0}^{M-1} \alpha_i \cdot e^{-j 2 \pi p \nu}\) où \(\nu\) est la fréquence réduite égale à \(\nu = f/F_e\) où \(F_e\) est la fréquence d’échantillonnage.

Conception du filtre sous Matlab

Plusieurs logiciels existent pour concevoir des filtres numériques. Nous utiliserons Matlab et sa boite à outils Filter Designer.

Matlab – Filter Designer

Lancez Matlab. Puis, dans l’onglet APPS, sélectionnez Filter Designer (dans la rubrique Signal Processing and Communications). Une nouvelle fenêtre s’ouvre alors :

Si vous n’avez pas ce type de fenêtre, cliquez sur la dernière icone en bas à gauche pour sélectionner l’outil de conception de filtre.

Cet outil vous permet de renseigner le gabarit du filtre souhaité :

  • Méthode de conception (Design Method) : IIR ou FIR, avec des choix particuliers pour chacun d’entre eux
  • Type de réponse : passe-bas (Lowpass), passe-haut (Highpass), passe-bande (Bandpass), coupe-bande (Bandstop)
  • Fréquence d’échantillonnage (notée Fs) – Cette fréquence est à adapter en fonction des capacités du microcontroleur qui effectuera le calcul des nouveaux échantillons du filtre.
  • Fréquences caractéristiques : coupure (notée Fc), bande-passante (notée Fpass), bande atténuée (notée Fstop) en Hz (par défaut)
  • Atténuation : dans la bande-passante (notée Apass), dans la bande atténuée (notée Astop) en dB (par défaut)

Il faudra ici choisir un filtre de type FIR.

Une fois les différents paramètres saisis, vous pouvez cliquer sur le bouton Design Filter.

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 :

Dans la fenêtre du haut, on peut voir la réponse fréquentielle du filtre ainsi obtenu. On peut également noter l’ordre du filtre (ici 52).

Il est encore possible de modifier certaines valeurs et de relancer la conception du filtre (bouton Design Filter). Cette réponse fréquentielle sera alors modifiée.

Génération des coefficients

L’étape précédente permet de créer le filtre souhaité et d’en vérifier son comportement fréquentiel attendu (via l’affichage de la réponse fréquentielle théorique du filtre conçu).

Il faut à présent obtenir de cette réponse les différents coefficients du filtre afin de pouvoir l’implémenter dans un calculateur.

Pour cela, il suffit d’aller dans le menu Targets et de sélectionner Generate C Headers….

Une nouvelle fenêtre s’ouvre permettant de sélectionner :

  • le nom du tableau qui contiendra les coefficients du filtre (noté Numerator)
  • le nom de la constante du nombre de coefficients du filtre (noté Numerator length)
  • le type de coefficients que l’on souhaite

Pour l’implémentation sur une carte Nucleo, il est préférable de passer dans un format de type : Single precision float.

En validant, l’application vous propose de sauvegarder un fichier header (extension .h) que vous pouvez renommer et stocker où vous le souhaitez. Ce fichier est de la forme suivante :

/*
 * Filter Coefficients (C Source) generated by the Filter Design and Analysis Tool
 * Generated by MATLAB(R) 9.1 and the Signal Processing Toolbox 7.3.
 * Generated on: 13-Oct-2017 14:23:25
 */

/*
 * Discrete-Time FIR Filter (real)
 * -------------------------------
 * Filter Structure  : Direct-Form FIR
 * Filter Length     : 53
 * Stable            : Yes
 * Linear Phase      : Yes (Type 1)
 */

/* General type conversion for MATLAB generated C-code  */
#include "tmwtypes.h"
/* 
 * Expected path to tmwtypes.h 
 * /usr/local/MATLAB/R2016b/extern/include/tmwtypes.h 
 */
/*
 * Warning - Filter coefficients were truncated to fit specified data type.  
 *   The resulting response may not match generated theoretical response.
 *   Use the Filter Design & Analysis Tool to design accurate
 *   single-precision filter coefficients.
 */
const int BL = 53;
const real32_T B[53] = {
  -1.923207674e-05,8.679173334e-05,-0.0001035222012,-9.687029524e-05,0.0004889828269,
  -0.0006440563011,-8.671574244e-19, 0.001367349061,-0.002238002373, 0.000988851185,
    0.00250719348,  -0.0055896868, 0.004243299831, 0.002845831681, -0.01105924882,
    0.01178205013,-3.774590571e-17,  -0.0181965474,  0.02652300335, -0.01072031818,
   -0.02558417059,  0.05551689863, -0.04286225513, -0.03121924214,   0.1481076032,
    -0.2561196983,   0.2999954224,  -0.2561196983,   0.1481076032, -0.03121924214,
   -0.04286225513,  0.05551689863, -0.02558417059, -0.01072031818,  0.02652300335,
    -0.0181965474,-3.774590571e-17,  0.01178205013, -0.01105924882, 0.002845831681,
   0.004243299831,  -0.0055896868,  0.00250719348, 0.000988851185,-0.002238002373,
   0.001367349061,-8.671574244e-19,-0.0006440563011,0.0004889828269,-9.687029524e-05,
  -0.0001035222012,8.679173334e-05,-1.923207674e-05
};

Ce code contient un tableau avec les coefficients du filtre, nommé (par défaut) B dans cet exemple, et le nombre de coefficients associés, nommé BL dans cet exemple.

Implémentation sur Nucleo via MBED

Principe

Ce type de système numérique, ici un filtre numérique, nécessite de travailler sur des échantillons numériques du signal d’entrée et restitue en sortie un échantillon numérique.

Afin de pouvoir réaliser cette opération de filtrage sur des signaux réels (et donc analogiques), il est indispensable que ce traitement numérique de l’information soit accompagné de modules de conversions analogique-numérique (ADC – Analog-Digital Converter) et numérique-analogique (DAC – Digital-Analog Converter).
Il faudra donc mettre en oeuvre la chaîne d’acquisition représentée dans la figure suivante pour pouvoir réaliser ce système numérique.

Les cartes Nucleo L476RG intègrent déjà ces deux composants (voir les tutoriels Récupérer un signal analogique et Générer une tension analogique).

Le programme de traitement numérique devra alors réaliser les tâches suivantes (voir également figure suivante) :

  • Echantilloner le signal d’entrée, à intervalle régulier
  • Effectuer le calcul de la sortie à partir du nouvel échantillon, des précédents (stockés en mémoire) et des coefficients du filtre
  • Convertir la valeur de sortie de ce module de traitement en valeur analogique

Selon le type de filtre à réaliser, en particulier de son ordre, le temps de calcul de la phase de traitement peut différer d’un filtre à l’autre. Il faudra alors garantir que la période d’échantillonnage est suffisamment longue pour que le système ait le temps de réaliser toutes les étapes précédentes. On pourra alors s’intéresser au temps d’exécution de ces calculs (voir le tutoriel Caractériser un traitement numérique).

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 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 FIR et les méthodes permettant la mise en place de la structure matérielle associée à un filtre de type FIR.

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_fir_instance_f32 : type permettant de construire un objet de type filtre FIR, dont les coefficients sont des flottants sur 32 bits
  • arm_fir_init_f32 : fonction permettant d’initialiser la structure matérielle d’un filtre FIR, dont les coefficients sont des flottants sur 32 bits
  • arm_fir_f32 : fonction permettant de calculer la nouvelle sortie d’un filtre FIR, 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.

/*
 * Coefficients du filtre
 */
const int BL = 53;
float32_t B[53] = {
  -1.923207674e-05,8.679173334e-05,-0.0001035222012,-9.687029524e-05,0.0004889828269,
  -0.0006440563011,-8.671574244e-19, 0.001367349061,-0.002238002373, 0.000988851185,
    0.00250719348,  -0.0055896868, 0.004243299831, 0.002845831681, -0.01105924882,
    0.01178205013,-3.774590571e-17,  -0.0181965474,  0.02652300335, -0.01072031818,
   -0.02558417059,  0.05551689863, -0.04286225513, -0.03121924214,   0.1481076032,
    -0.2561196983,   0.2999954224,  -0.2561196983,   0.1481076032, -0.03121924214,
   -0.04286225513,  0.05551689863, -0.02558417059, -0.01072031818,  0.02652300335,
    -0.0181965474,-3.774590571e-17,  0.01178205013, -0.01105924882, 0.002845831681,
   0.004243299831,  -0.0055896868,  0.00250719348, 0.000988851185,-0.002238002373,
   0.001367349061,-8.671574244e-19,-0.0006440563011,0.0004889828269,-9.687029524e-05,
  -0.0001035222012,8.679173334e-05,-1.923207674e-05
};

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

float32_t   state[BL];
arm_fir_instance_f32 fir; 

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 seconde ligne permet de créer un objet de type filtre FIR, 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_fir_init_f32(&fir, BL, B, state, 1); 

Les paramètres à lui fournir sont : l’objet de type filtre FIR (ici fir), le nombre de coefficients (ici BL), le tableau contenant les coefficients (ici B), le tableau permettant de stocker tous les échantillons nécessaires (ici 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_fir_f32(&fir, &in, &out, 1);

Les paramètres à lui fournir sont : l’objet de type filtre FIR (ici fir), 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é).

Filtrage numérique

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).

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

#define TE          0.00003
#define SAMPLE_RATE  1/TE
 
const int BL = 53;
float32_t B[53] = {
  -1.923207674e-05,8.679173334e-05,-0.0001035222012,-9.687029524e-05,0.0004889828269,
  -0.0006440563011,-8.671574244e-19, 0.001367349061,-0.002238002373, 0.000988851185,
    0.00250719348,  -0.0055896868, 0.004243299831, 0.002845831681, -0.01105924882,
    0.01178205013,-3.774590571e-17,  -0.0181965474,  0.02652300335, -0.01072031818,
   -0.02558417059,  0.05551689863, -0.04286225513, -0.03121924214,   0.1481076032,
    -0.2561196983,   0.2999954224,  -0.2561196983,   0.1481076032, -0.03121924214,
   -0.04286225513,  0.05551689863, -0.02558417059, -0.01072031818,  0.02652300335,
    -0.0181965474,-3.774590571e-17,  0.01178205013, -0.01105924882, 0.002845831681,
   0.004243299831,  -0.0055896868,  0.00250719348, 0.000988851185,-0.002238002373,
   0.001367349061,-8.671574244e-19,-0.0006440563011,0.0004889828269,-9.687029524e-05,
  -0.0001035222012,8.679173334e-05,-1.923207674e-05
};

float32_t   state[BL];
arm_fir_instance_f32 fir;   

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() {
    arm_fir_init_f32(&fir, BL, B, state, 1);
    
    tik.attach(&convert, TE);
    
    while(1) {
    }
}

void convert(void){
    in = mesure.read();
    arm_fir_f32(&fir, &in, &out, 1);
    
    analog_out.write(out);
}

Tutoriel lié

MInE Prototyper Prototyper avec Nucleo et MBED

Nucléo – Supprimer une bande fréquentielle