Calcul de moyenne mobile rapide et efficace en mémoire

Je cherche une solution économe en temps et en mémoire pour calculer une moyenne mobile en C. Je dois éviter la division, car je suis sur un PIC 16 qui ne dispose pas d'une unité de division dédiée.

Pour le moment, je stocke toutes les valeurs dans un tampon circulaire et stocke et met simplement à jour la somme à chaque nouvelle valeur. C'est vraiment efficace, mais utilise malheureusement la majeure partie de ma mémoire disponible ...

Comme d'autres l'ont mentionné, vous devriez envisager un filtre IIR (réponse impulsionnelle infinie) plutôt que le filtre FIR (réponse impulsionnelle finie) que vous utilisez actuellement. Il y a plus que cela, mais à première vue, les filtres FIR sont implémentés sous forme de convolutions explicites et de filtres IIR avec équations.

Le filtre IIR que j'utilise beaucoup dans les microcontrôleurs est un filtre passe-bas unipolaire. C'est l'équivalent numérique d'un simple filtre analogique RC. Pour la plupart des applications, celles-ci auront de meilleures caractéristiques que le filtre de boîte que vous utilisez. La plupart des utilisations d’un filtre de boîte que j’ai rencontrées résultent de l’absence de vigilance de la part de la classe de traitement du signal numérique et non du fait qu’elles ont besoin de leurs caractéristiques particulières. Si vous souhaitez simplement atténuer les hautes fréquences que vous savez être du bruit, un filtre passe-bas unipolaire est préférable. La meilleure façon d’implémenter numériquement l’un dans un microcontrôleur est généralement:

FILT <- FILT + FF (NOUVEAU - FILT)

FILT est un morceau d'état persistant. C'est la seule variable persistante dont vous avez besoin pour calculer ce filtre. NEW est la nouvelle valeur que le filtre est mis à jour avec cette itération. FF est la fraction de filtre qui ajuste la "lourdeur" du filtre. Regardez cet algorithme et voyez que pour FF = 0 le filtre est infiniment lourd puisque la sortie ne change jamais. Pour FF = 1, il n’ya vraiment aucun filtre car la sortie ne fait que suivre l’entrée. Les valeurs utiles sont entre les deux. Sur les petits systèmes, vous choisissez FF égal à 1/2 ^N sorte que la multiplication par FF puisse être accomplie comme un décalage à droite de N bits. Par exemple, FF pourrait être 1/16 et multiplié par FF, donc un décalage à droite de 4 bits. Sinon, ce filtre ne nécessite qu'une soustraction et une addition, bien que les nombres doivent généralement être plus larges que la valeur d'entrée (plus de précision numérique dans une section séparée ci-dessous).

Je prends habituellement les lectures A / D bien plus rapidement que nécessaire et applique deux de ces filtres en cascade. Il s’agit de l’équivalent numérique de deux filtres RC en série et s’atténue de 12 dB / octave au-dessus de la fréquence de décélération. Cependant, pour les lectures A / D, il est généralement plus pertinent d'examiner le filtre dans le domaine temporel en tenant compte de sa réponse progressive. Cela vous indique à quelle vitesse votre système verra un changement lorsque ce que vous mesurez change.

Pour faciliter la conception de ces filtres (c’est-à-dire choisir le FF et choisir le nombre de filtres à utiliser en cascade), j’utilise mon programme FILTBITS. Vous spécifiez le nombre de bits de décalage pour chaque FF dans la série de filtres en cascade, qui calcule la réponse à l'étape et les autres valeurs. En fait, je le fais habituellement via mon script wrapper PLOTFILT. Cela lance FILTBITS, qui crée un fichier CSV, puis trace le fichier CSV. Par exemple, voici le résultat de "PLOTFILT 4 4":

Les deux paramètres à PLOTFILT signifient qu'il y aura deux filtres en cascade du type décrit ci-dessus. Les valeurs de 4 indiquent le nombre de bits de décalage pour réaliser la multiplication par FF. Les deux valeurs FF sont donc 1/16 dans ce cas.

La trace rouge correspond à la réponse pas à pas de l'unité et constitue l'élément principal à examiner. Par exemple, cela vous indique que si l'entrée change instantanément, la sortie du filtre combiné s'installera à 90% de la nouvelle valeur en 60 itérations. Si vous vous souciez du temps de stabilisation de 95%, vous devez attendre environ 73 itérations, et pour le temps de stabilisation de 50%, seulement 26 itérations.

Le tracé vert vous montre le résultat d’un pic unique d’amplitude complète. Cela vous donne une idée de la suppression du bruit aléatoire. Il semble qu'aucun échantillon ne provoque plus de 2,5% de changement dans la sortie.

La trace bleue donne une impression subjective de ce que fait ce filtre avec le bruit blanc. Ce test n’est pas rigoureux, car rien ne garantit exactement le contenu exact des nombres aléatoires choisis comme entrée de bruit blanc pour cette série de PLOTFILT. C'est seulement pour vous donner une idée approximative de la quantité de poudre que vous allez écraser et de sa douceur.

PLOTFILT, peut-être FILTBITS, et de nombreux autres éléments utiles, en particulier pour le développement de microprogrammes PIC, sont disponibles dans la version logicielle des outils de développement PIC sur la page Téléchargements de logiciels .

Ajouté à propos de la précision numérique

Je vois dans les commentaires et maintenant une nouvelle réponse qu'il y a un intérêt à discuter du nombre de bits nécessaires pour implémenter ce filtre. Notez que la multiplication par FF créera de nouveaux bits de Log ₂ (FF) en dessous du point binaire. Sur les petits systèmes, FF est généralement choisi égal à 1/2 ^{N, de} sorte que cette multiplication est effectivement réalisée par un décalage à droite de N bits.

FILT est donc généralement un nombre entier à point fixe. Notez que cela ne change en rien le calcul du point de vue du processeur. Par exemple, si vous filtrez les lectures A / D 10 bits et que N = 4 (FF = 1/16), vous avez besoin de 4 bits de fraction inférieurs aux lectures A / D de 10 bits d'entiers. Dans la plupart des processeurs, vous feriez des opérations sur les entiers 16 bits en raison des lectures A / N 10 bits. Dans ce cas, vous pouvez toujours effectuer exactement les mêmes opérations sur les entiers 16 bits, mais commencez par les lectures A / D décalées de 4 bits. Le processeur ne connaît pas la différence et n'en a pas besoin. Faire le calcul sur des entiers 16 bits entiers fonctionne que vous considériez qu’il s’agit d’entiers de 12,4 points fixes ou de vrais entiers 16 bits (16,0 points fixes).

En général, vous devez ajouter N bits à chaque pôle de filtrage si vous ne souhaitez pas ajouter de bruit en raison de la représentation numérique. Dans l'exemple ci-dessus, le deuxième filtre de deux aurait 10 + 4 + 4 = 18 bits pour ne pas perdre d'informations. En pratique, sur une machine 8 bits, cela signifie que vous utiliseriez des valeurs 24 bits. Techniquement, seul le deuxième pôle sur deux aurait besoin de la valeur la plus large, mais pour simplifier les microprogrammes, j'utilise généralement la même représentation, et donc le même code, pour tous les pôles d'un filtre.

En général, j’écris un sous-programme ou une macro pour effectuer une opération de filtrage, puis l’applique à chaque pôle. Que ce soit un sous-programme ou une macro dépend de l'importance des cycles ou de la mémoire programme dans ce projet particulier. Quoi qu’il en soit, j’utilise un état de travail pour passer NEW dans le sous-programme / macro, qui met à jour FILT, mais le charge également dans le même état de travail de NEW. Cela facilite l’application de plusieurs pôles puisque le FILT mis à jour d’un pôle est le NOUVEAU du suivant. Lorsqu’un sous-programme, il est utile d’avoir un pointeur sur FILT à l’entrée, qui est mis à jour juste après FILT à la sortie. De cette façon, le sous-programme fonctionne automatiquement sur des filtres consécutifs en mémoire s’il est appelé plusieurs fois. Avec une macro, vous n'avez pas besoin d'un pointeur puisque vous transmettez l'adresse pour pouvoir l'utiliser à chaque itération.

Exemples de code

Voici un exemple de macro décrite ci-dessus pour un PIC 18:

////////////////////////////////////////////////////////// ////////////////////////////////////
//
// Macro FILTER filt
//
// Met à jour un pôle de filtrage avec la nouvelle valeur dans NEWVAL. NEWVAL est mis à jour pour
// contient la nouvelle valeur filtrée.
//
// FILT est le nom de la variable d'état du filtre. Il est supposé être de 24 bits
// large et dans la banque locale.
//
// La formule pour mettre à jour le filtre est la suivante:
//
// FILT <- FILT + FF (NEWVAL - FILT)
//
// La multiplication par FF est réalisée par un décalage à droite des bits FILTBITS.
//
/ filtre macro
  /écrire
         dbankif lbankadr
         movf [arg 1] +0, w; NEWVAL <- NEWVAL - FILT
         subwf newval + 0
         movf [arg 1] +1, w
         souswfb newval + 1
         movf [arg 1] +2, w
         souswfb newval + 2

  /écrire
  / loop n filtbits; une fois pour chaque bit pour déplacer NEWVAL à droite
         rlcf newval + 2, w; décale NEWVAL d'un bit à droite
         rrcf newval + 2
         rrcf newval + 1
         rrcf newval + 0
    / endloop

  /écrire
         movf newval + 0, w; ajoute une valeur décalée dans le filtre et enregistre dans NEWVAL
         addwf [arg 1] +0, w
         movwf [arg 1] +0
         movwf newval + 0

         movf newval + 1, w
         addwfc [arg 1] +1, w
         movwf [arg 1] +1
         movwf newval + 1

         movf newval + 2, w
         addwfc [arg 1] +2, w
         movwf [arg 1] +2
         movwf newval + 2
  / endmac

Et voici une macro similaire pour un PIC 24 ou dsPIC 30 ou 33:

////////////////////////////////////////////////////////// ////////////////////////////////////
//
// Macro FILTER ffbits
//
// Met à jour l'état d'un filtre passe-bas. La nouvelle valeur d'entrée est dans W1: W0
// et l'état du filtre à mettre à jour est pointé par W2.
//
// La valeur de filtre mise à jour sera également renvoyée dans W1: W0 et W2 indiqueront
// vers la première mémoire après l'état du filtre. Cette macro peut donc être
// invoqué successivement pour mettre à jour une série de filtres passe-bas en cascade.
//
// La formule de filtre est:
//
// FILT <- FILT + FF (NOUVEAU - FILT)
//
// où la multiplication par FF est effectuée par un décalage arithmétique à droite de
// FFBITS.
//
// ATTENTION: W3 est mis à la corbeille.
//
/ filtre macro
  / var new ffbits entier = [arg 1]; récupère le nombre de bits à décaler

  /écrire
  / write "; Effectue un filtrage passe-bas à un pôle, bits de décalage =" ffbits
  /écrire " ;"

         sous w0, [w2 ++], w0; NOUVEAU - FILT -> W1: W0
         subb w1, [w2--], w1

         lsr w0, # [v ffbits], w0; décale le résultat dans W1: W0 à droite
         sl w1, # [- 16 ffbits], w3
         ior w0, w3, w0
         asr w1, # [vffbits], w1

         ajouter w0, [w2 ++], w0; ajouter FILT pour obtenir le résultat final dans W1: W0
         addc w1, [w2--], w1

         mov w0, [w2 ++]; écrit le résultat dans l'état du filtre, avance le pointeur
         mov w1, [w2 ++]

  /écrire
  / endmac

Ces deux exemples sont implémentés sous forme de macros à l'aide de mon préprocesseur d'assembleur PIC , qui est plus capable que l'une ou l'autre des installations de macro intégrées.

Si vous pouvez vivre avec la limitation d'une puissance de deux éléments à la moyenne (par exemple 2,4,8,16,32 etc.), la division peut être facilement et efficacement effectuée sur un micro à faible performance sans division dédiée car il peut être fait comme un peu de décalage. Chaque changement de droite est une puissance de deux, par exemple:

avg = sum >> 2; //divide by 2^2 (4)

ou

avg = sum >> 3; //divide by 2^3 (8)

etc.

Il existe une réponse pour un filtre à moyenne mobile véritable (ou "filtre de boîte à lettres") avec moins de besoins en mémoire, si cela ne vous dérange pas de sous-échantillonner. C'est ce qu'on appelle un filtre en peigne intégrateur en cascade (CIC). L'idée est que vous avez un intégrateur dont vous prenez les différences sur une période de temps donnée, et le principal moyen d'économiser de la mémoire est qu'en réduisant l'échantillonnage, il n'est pas nécessaire de stocker toutes les valeurs de l'intégrateur. Il peut être implémenté en utilisant le pseudocode suivant:

function out = filterInput(in)
{
   const int decimationFactor = /* 2 or 4 or 8 or whatever */;
   const int statesize = /* whatever */
   static int integrator = 0;
   static int downsample_count = 0;
   static int ringbuffer[statesize];
   // don't forget to initialize the ringbuffer somehow
   static int ringbuffer_ptr = 0;
   static int outstate = 0;

   integrator += in;
   if (++downsample_count >= decimationFactor)
   {
     int oldintegrator = ringbuffer[ringbuffer_ptr];
     ringbuffer[ringbuffer_ptr] = integrator;
     ringbuffer_ptr = (ringbuffer_ptr + 1) % statesize;
     outstate = (integrator - oldintegrator) / (statesize * decimationFactor);
   }
   return outstate;
}

Votre longueur moyenne mobile effective est, decimationFactor*statesizemais il vous suffit de conserver des statesizeéchantillons. Il est évident que vous pouvez obtenir de meilleures performances si votre statesizeet decimationFactorsont des puissances de 2, de sorte que les opérateurs de la division et le reste se remplacés par des quarts de travail et masque ands.

Post-scriptum: Je suis d'accord avec Olin sur le fait que vous devriez toujours considérer les filtres IIR simples avant un filtre à moyenne mobile. Si vous n'avez pas besoin des valeurs de fréquence zéro d'un filtre de wagon couvert, un filtre passe-bas à 1 ou 2 pôles fonctionnera probablement très bien.

D'autre part, si vous filtrez à des fins de décimation (en prenant une entrée à fréquence d'échantillonnage élevée et en la moyenant pour une utilisation par un processus à fréquence réduite), un filtre CIC peut être exactement ce que vous recherchez. (surtout si vous pouvez utiliser etatsize = 1 et éviter le tampon circulaire avec une seule valeur d'intégrateur précédente)

Il existe une analyse approfondie des calculs à l'aide du filtre IIR de premier ordre décrit par Olin Lathrop à propos de l' échange de pile de traitement des signaux numériques (inclut de nombreuses images jolies.) L'équation de ce filtre IIR est la suivante:

y [n] = αx [n] + (1-α) y [n-1]

Cela peut être implémenté en utilisant uniquement des entiers et aucune division en utilisant le code suivant (il faudra peut-être un certain débogage car je tapais à partir de la mémoire.)

/**
*  @details    Implement a first order IIR filter to approximate a K sample 
*              moving average.  This function implements the equation:
*
*                  y[n] = alpha * x[n] + (1 - alpha) * y[n-1]
*
*  @param      *filter - a Signed 15.16 fixed-point value.
*  @param      sample - the 16-bit value of the current sample.
*/

#define BITS 2      ///< This is roughly = log2( 1 / alpha )

short IIR_Filter(long *filter, short sample)
{
    long local_sample = sample << 16;

    *filter += (local_sample - *filter) >> BITS;

    return (short)((*filter+0x8000) >> 16);     ///< Round by adding .5 and truncating.
}

Ce filtre se rapproche de la moyenne mobile des K derniers échantillons en définissant la valeur de alpha sur 1 / K. Faites cela dans le code précédent #defineen BITSallant à LOG2 (K), c’est-à-dire que K = 16 est réglé BITSsur 4, que K = 4 est réglé BITSsur 2, etc.

(Je vérifierai le code répertorié ici dès que j'aurai un changement et modifierai cette réponse si nécessaire.)

Voici un filtre passe-bas unipolaire (moyenne mobile, avec fréquence de coupure = CutoffFrequency). Très simple, très rapide, fonctionne très bien et presque pas de surcharge de mémoire.

Remarque: Toutes les variables ont une portée au-delà de la fonction de filtre, à l'exception de celles passées dans newInput.

// One-time calculations (can be pre-calculated at compile-time and loaded with constants)
DecayFactor = exp(-2.0 * PI * CutoffFrequency / SampleRate);
AmplitudeFactor = (1.0 - DecayFactor);

// Filter Loop Function ----- THIS IS IT -----
double Filter(double newInput)
{
   MovingAverage *= DecayFactor;
   MovingAverage += AmplitudeFactor * newInput;

   return (MovingAverage);
}

Remarque: il s'agit d'un filtre à une étape. Plusieurs étapes peuvent être mises en cascade pour augmenter la netteté du filtre. Si vous utilisez plus d'une étape, vous devrez ajuster DecayFactor (en fonction de la fréquence de coupure) pour compenser.

Et évidemment, tout ce dont vous avez besoin, c'est de placer ces deux lignes n'importe où, elles n'ont pas besoin de leur propre fonction. Ce filtre a un temps de montée avant que la moyenne mobile ne représente celle du signal d’entrée. Si vous devez éviter ce temps de montée, vous pouvez simplement initialiser MovingAverage à la première valeur de newInput au lieu de 0 et espérer que le premier newInput ne sera pas une valeur aberrante.

(CutoffFrequency / SampleRate) a une plage comprise entre 0 et 0,5. DecayFactor est une valeur comprise entre 0 et 1, généralement proche de 1.

Les flotteurs à simple précision conviennent à la plupart des choses, je préfère simplement les doubles. Si vous devez vous en tenir à des nombres entiers, vous pouvez convertir DecayFactor et Facteur d’amplitude en nombres entiers fractionnaires, dans lesquels le numérateur est stocké en tant qu’entier et le dénominateur est une puissance entière égale à 2 (vous pouvez donc décaler le bit dénominateur plutôt que de devoir se diviser pendant la boucle de filtrage). Par exemple, si DecayFactor = 0,99 et que vous souhaitez utiliser des entiers, vous pouvez définir DecayFactor = 0,99 * 65536 = 64881. Et à chaque fois que vous multipliez par DecayFactor dans votre boucle de filtre, il suffit de décaler le résultat >> 16.

Pour plus d'informations à ce sujet, un excellent livre en ligne, chapitre 19 sur les filtres récursifs: http://www.dspguide.com/ch19.htm

PS Pour le paradigme de la moyenne mobile, une approche différente de la configuration de DecayFactor et AmplitudeFactor peut être plus pertinente pour vos besoins. Supposons que vous souhaitiez la moyenne des 6 éléments précédents, en faisant la moyenne, en ajoutant discrètement 6 éléments et en les divisant par 6, vous pouvez donc régler AmplitudeFactor sur 1/6 et DecayFactor sur (1.0 - AmplitudeFactor).

Vous pouvez approximer une moyenne mobile pour certaines applications avec un simple filtre IIR.

le poids est 0..255 valeur, valeurs élevées = délai plus court pour la moyenne

Value = (newvalue * weight + value * (256-weight)) / 256

Pour éviter les erreurs d'arrondi, la valeur doit normalement être longue et vous devez uniquement utiliser des octets d'ordre supérieur en tant que valeur "réelle".

Tous les autres ont formulé des commentaires approfondis sur l'utilité de l'IIF par rapport au FIR et à la division par deux. Je voudrais juste donner quelques détails de mise en œuvre. Ce qui suit fonctionne bien sur de petits microcontrôleurs sans FPU. Il n'y a pas de multiplication, et si vous gardez N une puissance de deux, toute la division est un transfert de bits à cycle unique.

Tampon en anneau FIR de base: conserve un tampon en cours d’exécution comportant les N dernières valeurs et un état récapitulatif en cours d'exécution de toutes les valeurs du tampon. Chaque fois qu'un nouvel échantillon entre, soustrayez la valeur la plus ancienne dans la mémoire tampon de SUM, remplacez-le par le nouvel échantillon, ajoutez le nouvel échantillon à SUM et générez SUM / N.

unsigned int Filter(unsigned int sample){
    static unsigned int buffer[N];
    static unsigned char oldest = 0;
    static unsigned long sum;

    sum -= buffer[oldest];
    sum += sample;
    buffer[oldest] = sample;
    oldest += 1;
    if (oldest >= N) oldest = 0;

    return sum/N;
}

Mémoire en anneau IIR modifiée: conservez la somme cumulée des N dernières valeurs. Chaque fois qu'un nouvel échantillon entre, SUM - = SUM / N, ajoute le nouvel échantillon et génère SUM / N.

unsigned int Filter(unsigned int sample){
    static unsigned long sum;

    sum -= sum/N;
    sum += sample;

    return sum/N;
}

Comme l'a dit mikeselectricstuff , si vous devez réellement réduire vos besoins en mémoire et que votre réponse impulsionnelle est exponentielle (au lieu d'une impulsion rectangulaire), je choisirais un filtre à moyenne mobile exponentielle . Je les utilise beaucoup. Avec ce type de filtre, vous n'avez besoin d'aucun tampon. Vous n'êtes pas obligé de stocker N échantillons précédents. Juste un. Ainsi, vos besoins en mémoire sont réduits d'un facteur N.

En outre, vous n'avez besoin d'aucune division pour cela. Seulement des multiplications. Si vous avez accès à l'arithmétique à virgule flottante, utilisez des multiplications à virgule flottante. Sinon, multiplier et décaler les nombres entiers vers la droite. Cependant, nous sommes en 2012 et je vous recommanderais d'utiliser des compilateurs (et des MCU) vous permettant de travailler avec des nombres à virgule flottante.

En plus d'être plus efficace en termes de mémoire et plus rapide (vous n'avez pas à mettre à jour les éléments dans un tampon circulaire), je dirais que c'est aussi plus naturel , car une réponse impulsionnelle exponentielle correspond mieux au comportement de la nature, dans la plupart des cas.

Un problème avec le filtre IIR, presque touché par @olin et @sercercat mais apparemment ignoré par d’autres, est que l’arrondi vers le bas introduit une certaine imprécision (et éventuellement un biais / troncature): en supposant que N est une puissance de deux et que seul l’arithmétique entière est utilisé, le décalage à droite élimine systématiquement les LSB du nouvel échantillon. Cela signifie que, quelle que soit la durée de la série, la moyenne ne les prendra jamais en compte.

Par exemple, supposons une série qui décroît lentement (8,8,8, ..., 8,7,7,7, ... 7,6,6,) et supposons que la moyenne est bien 8 au début. Le premier échantillon "7" portera la moyenne à 7, quelle que soit la force du filtre. Juste pour un échantillon. Même histoire pour 6 ans, etc. Imaginons maintenant le contraire: la série monte. La moyenne restera sur 7 pour toujours, jusqu'à ce que l'échantillon soit assez grand pour le faire changer.

Bien sûr, vous pouvez corriger le "biais" en ajoutant 1/2 ^ N / 2, mais cela ne résoudra pas vraiment le problème de précision: dans ce cas, la série décroissante restera indéfiniment à 8 jusqu'à ce que l'échantillon soit 8-1. / 2 ^ (N / 2). Pour N = 4 par exemple, tout échantillon supérieur à zéro conservera la moyenne inchangée.

Je pense qu'une solution à cela impliquerait de conserver un accumulateur des LSB perdus. Mais je ne me suis pas rendu assez loin pour avoir le code prêt, et je ne suis pas sûr que cela ne nuirait pas au pouvoir de l'IIF dans d'autres cas de séries (par exemple, si 7,9,7,9 serait en moyenne à 8 alors) .

@Olin, votre cascade en deux étapes aurait également besoin d'explications. Voulez-vous dire conserver deux valeurs moyennes avec le résultat de la première entrée dans la seconde à chaque itération? Quel est l'avantage de cela?