[algo] Comment calculer une moyenne glissante sur un microcontrôleur à faibles ressources

Pour commencer, je dois dire que l’idée pour cet article m’a été suggérée par Audran (@alf_arobase) lors de discussions à l’Electrolab sur son super projet que je ne dévoilerai pas sans son autorisation! Mais c’est un truc très bien! Bref, cet article est surtout une synthèse des trucs dont on a discuté hier soir, et que j’ai relativement mal expliqué. Avec du recul, ça va mieux.

Le problème

Supposons que l’on ait à mesurer en temps réel un signal analogique à l’aide d’un petit microcontrôleur, pour, par exemple, déterminer quand ce signal passe sous un seuil, et déclencher un signal. Ceci peut être fait de manière triviale, en effectuant une mesure avec l’ADC, en la comparant avec le seuil, puis en décidant l’état du signal de sortie.

Ceci, en réalité, n’est qu’une théorie. Dans la réalité, plusieurs effets vont rendre cette technique impraticable:

  • Le signal à mesurer peut être très bruité;
  • L’ADC peut manquer de précision et sortir des valeurs de mesure différentes alors que le signal reste constant.

En conséquence, quand le signal sera proche de la limite, il va varier rapidement entre les états “signal détecté” et “signal non détecté”. Comme les mesures se feront en général à un rythme rapide, les variations du signal de sortie seront également très rapides, et probablement assez aléatoires.

Que peut on faire? Idéalement, on voudrait se passer de toutes ces petites variations rapides, pour retenir la valeur moyenne du signal à un moment donné. On peut aussi appliquer un filtre passe-bas au signal de l’ADC, car il est équivalent de dire que seules les basses fréquences (les variations lentes) du signal de l’ADC nous intéressent.

Première option: y’a moyen de moyenner (le signal)

La première chose qu’on peut faire, c’est diminuer l’influence des imperfections de l’ADC, c’est le plus facile.
Au lieu de mettre à jour le signal de sortie à chaque fois qu’un échantillon ADC est lu, on peut lire 4 échantillons, ou N, en faire la moyenne, et utiliser cette moyenne pour prendre la décision sur le signal de sortie. Dans ce cas, chaque fois qu’on moyenne 2^n valeurs d’ADC, on gagne n bits de résolution.

Cette méthode s’appelle la décimation. Ça s’invente pas.

Exemple, si je moyenne 4 valeurs, la moyenne pourra prendre, après la virgule, les valeurs x.00, x.25, x.50, ou x.75. Quatre valeurs après la virgule peuvent être codées avec deux bits, c’est deux bits de gagnés sur la résolution, au prix d’une vitesse d’acquisition plus faible. Autre manière de voir: additionner 4 nombres de M bits nécessite M+2 bits pour coder le résultat.

En embarqué, surtout avec des microcontrôleurs 8 bits, on n’aime pas les virgules, c’est donc un bon exemple d’application des calculs en virgule fixe. Faire la moyenne de N valeurs revient à les additionner, et c’est tout. Le résultat aura un dénominateur de N. Si on veut une moyenne du même ordre de grandeur que les résultats de mesure initiaux, on peut alors utiliser une division entière. Si on est malin, on fait une moyenne sur 2^N valeurs, et la chère division se transforme en simple décalage à droite.

Deuxième option : calculer une moyenne mobile

La moyenne c’est bien, mais cela pose un problème. Pour avoir des valeurs bien filtrées, il faut faire une moyenne sur beaucoup de valeurs, par exemple 64 ou 128, ce qui va diviser le “débit virtuel” de l’ADC moyenné par autant, et peut devenir vraiment un gros problème. De plus, chaque échantillon décimé va représenter un groupe d’acquisitions, une “case”, mais sans lien continu entre elles.

La moyenne mobile va nous aider.

La moyenne mobile, on la rencontre très fréquemment sur les graphes boursiers. En effet, les cours de bourse varient très chaotiquement, et un moyennage permet d’y voir plus clairement des tendances. Mais attention: dans ce domaine, faire une moyenne selon la première méthode reviendrait à utiliser un cours journalier pour calculer un cours moyen par semaine. On comprend bien que pour calculer des tendances au jour le jour, on ne peut pas se contenter d’une valeur toute les semaines, même si on tente de l’interpoler.

L’astuce consiste alors à calculer chaque jour une moyenne sur les N jours précédents. Chaque jour, la nouvelle moyenne prend en compte une partie des valeurs qui ont déjà été utilisées la veille, et ajoute l’information du nouveau jour. Ainsi, chaque jour, chaque valeur prend en compte les valeurs précédentes, mais cela atténue les variations rapides, puisque les données de chaque jour interviennent plusieurs fois dans le calcul. Voici à quoi cela ressemble sur le cours de, disons, RedHat pour rester dans l’open source. Ici j’ai demandé à tracer la moyenne mobile sur 50 jours (un classique). J’ai aussi tracé le graphe des valeurs par semaine (donc décimées), je vous l’affiche d’abord:

Cours de RedHat sur 1 an, décimation par semaine
Cours de RedHat sur 1 an, décimation par semaine
Cours de RedHat sur 1 an, décimation par jour, et moyenne mobile sur 50 jours
Cours de RedHat sur 1 an, décimation par jour, et moyenne mobile sur 50 jours

Bon, en bref: la décimation est effectivement moins bruitée que le graphe quotidien (moins de variations rapides). Mais elle contient moins de points (évidemment, 52 semaines par an au lieu de 365 jours). Et si votre seuil de basculement était à 50 dollars, eh bien les basculements seraient toujours aussi nombreux que sans filtrage. La décimation n’aurait eu aucun effet sur ce point.

Par contre la moyenne mobile, en rouge, varie plus “sagement”, tout en reproduisant les “tendances” du signal bruité, et surtout, elle existe pour chaque point de la courbe! Elle ne croise que 2 fois le seuil à $50.

A quoi cela va-t-il ressembler dans le cas de notre ADC? Eh bien c’est strictement la même chose. Chaque fois qu’on fait l’acquisition d’un échantillon, on en calcule la moyenne mobile, et on utilise cette valeur pour prendre la décision sur le signal de sortie.

Moyenne mobile: méthode naïve.

De manière simplifiée, on va avoir besoin d’un tableau pour stocker l’historique des N dernières valeurs. On aura besoin d’un index pour se souvenir de notre position, et d’une routine qui calcule la moyenne.

#define Nbits 4
#define N (1<<Nbits) //Du coup, N est une puissance de 2, ici 16
uchar hist[N];
uchar hist_index;

void init(void) {
    for(hist_index=0;hist_index<N;hist_index++) {
        hist[hist_index]=0;
    }
    hist_index=0;
}

void adc_interrupt(void) {
    uchar sample,i,filtered;
    ushort mm;

    //pseudo fonction pour lire l'ADC
    sample = adc_read();

    //stockage
    hist[hist_index] = sample;
    hist_index++;
    if(hist_index==N) { //gestion buffer circulaire
        hist_index = 0;
    }

    //Calcul brut de la MM
    mm=0;
    for(i=0;i<N;i++) {
        mm += hist[i];
    }

    //maintenant on peut jouer avec la moyenne mobile mm
    //si on veut qu'elle soit du même ordre de grandeur que les samples:
    filtered = mm >> Nbits
}

Cette méthode marche. Mais si N devient grand, le nombre d’additions à réaliser pour le calcul de mm devient énorme, et la routine d’interruption prendra énormément de temps. Il ne faut pas utiliser ce code, il est nul.

Moyenne mobile: Méthode pratique

On va repartir sur quelques bouts de maths : pas bien difficile, des suites. Imaginez des tableaux et des indices, c’est la même chose.

Définissons le signal en sortie de l’ADC comme une suite numérique.

s(n) == valeur lue de l’ADC à l’étape N

Appelons K le nombre de valeurs qu’on veut utiliser dans notre moyenne.

La moyenne mobile est alors définie comme

mm(n) = somme des s(n) pour les K dernières valeurs / K

Oublions tout de suite la division par K (elle est triviale et optionnelle, surtout si K est une puissance de 2)

Donc, on fait la somme des K dernières valeurs :

mm(n) = s(n) + s(n-1) + s(n-2) … + s(n-K+1)

Si cette expression vous perturbe, vous n’avez qu’à imaginer que s(n) est aussi définie pour des valeurs négatives de n, ou alors, qu’elle vaut zéro en ces points. En tout cas, ne vous prenez pas le chou avec ça. Dans le programme ça disparaitra.

C’est en gros le calcul qu’on a fait dans l’algo naïf.

 

Prenons un peu de recul et écrivons le calcul que nous allons faire à l’étape suivante pour calculer mm(n+1) :

mm(n+1) = s(n+1) + s(n-1+1) + s(n-2+1) … + s(n-K+1+1)

soit:

mm(n+1) = s(n+1) + s(n) + … + s(n-K+2)

Allez, ajoutons, pour le fun des maths, s(n-K+1) de chaque coté. Ca ne change rien n’est ce pas? Si A=B, alors A+truc = B+même_truc.

mm(n+1) + s(n-K+1) = s(n+1) + s(n) + … + s(n-K+2)+ s(n-K+1)

Quelle surprise… Mais c’est presque cette vieille mm(n) ! Mais si:

mm(n+1) + s(n-K+1) = s(n+1) + s(n) + … + s(n-K+2)+ s(n-K+1)

Au boulot:

mm(n+1) + s(n-K+1) = s(n+1) + mm(n)

Maintenant, on enlève s(n-K+1) de chaque coté :

mm(n+1) = s(n+1) + mm(n) – s(n-K+1)

C’est gagné. Observons ce qu’on a obtenu:

 

La nouvelle moyenne mobile, c’est  l’ancienne moyenne mobile, plus la valeur la plus récente, moins la plus ancienne.

 

Et ça, c’est vachement cool, parce que ça veut dire que dans mon programme, je n’ai plus que DEUX additions à faire à chaque IRQ, au lieu de calculer toute la somme à chaque fois!

Concrètement on a besoin de quoi? La valeur plus récente, on vient de la lire depuis l’ADC. Il faut se souvenir des anciennes valeurs. Combien d’anciennes valeurs? Exactement K. Après K étapes, ma plus ancienne valeur sera la valeur récente actuelle.

  1. A l’init, ma moyenne vaut zéro, et mes anciennes valeurs valent zéro (ça règle le problème des indices négatifs).
  2. A chaque irq, je lis un sample de l’ADC.
  3. Je mets à jour ma moyenne mobile
  4. Je stocke ce sample dans le tableau des anciennes valeurs
  5. Profit!

Allez, je termine le suspense avec le petit programme d’exemple qui implémente ça:

#define Nbits 4
#define N (1<<Nbits) //Du coup, N est une puissance de 2, ici 16
uchar hist[N];
uchar plus_ancien; //(@skywodd: pas volatile, touchée que dans l'irq après l'init)
volatile ushort mm; //stockage pour la moyenne, doit être une var globale

void init(void) {
    //effacer les valeurs anciennes
    for(plus_ancien=0;plus_ancien<N;plus_ancien++) {
        hist[plus_ancien] = 0;
    }
    //on commencera à stocker à cet offset
    plus_ancien = 0;
    mm = 0;
}

void irq_adc(void) {
    //acquisition
    sample = acd_read();

    //calcul MM
    mm = mm + sample - hist[plus_ancien];

    //cette plus ancienne valeur n'est plus utile, on y stocke la plus récente
    hist[plus_ancien] = sample;
    plus_ancien ++;
    if(plus_ancien==N) { //gestion du buffer circulaire
        plus_ancien=0;
    }
    //La MM est déjà prête à utiliser
}

Au début, avant d’avoir effectivement lu N (ou K) valeurs, la mm va augmenter lentement, il y a donc un petit régime transitoire avant que les choses se stabilisent. C’est inévitable, c’est une caractéristique des filtres.

Certains vont me dire, mais ça marche comme un filtre ! Eh bah oui, on a un genre de filtre passe bas à réponse impulsionnelle infinie (un mot compliqué pour dire que si tous vos samples ADC sont nuls sauf le premier, (et que la précision est infinie), alors tous les samples qui suivent auront une valeur non nulle…).

Les pros du traitement de signal me contrediront sans doute, mais c’est  pas grave. “Avec les mains”, ça marche, le signal est lissé.

Et voila, maintenant vous saurez dire que vous savez utiliser un ATTiny (ou autre :)) comme un DSP !

La combinaison

Bien sûr , si votre source de données a un débit élevé, rien ne vous empêche de calculer une moyenne mobile de valeurs préalablement décimées… Je sais pas trop si il y a quelque chose à y gagner…

 

A+ pour de nouveaux algos!