Filtre numerique RII, RIF….digital filter...FFT...atmega328, ESP32

un résumé de ce poste qui est parti dans tous les sens est sur ce lien

un wiki suivant a été suppimé
https://digital-filter-atmega328-esp32.fandom.com/fr/wiki/Wiki_Filtre_numerique_digital_filter,_RII,_RIF,_atmega328,_ESP32

S’il y a de nombreux cours, livre sur internet sur les filtres numériques, il n’y a pas beaucoup d’exemples sur le net, ni de vulgarisation.

Voici quelques lien sur les filtres :

http://www.f-legrand.fr/scidoc/docimg/numerique/filtre/filtrenum/filtrenum.html
il y a sur ce lien, sur le filtrage avec Arduino
http://www.f-legrand.fr/scidoc/docimg/sciphys/arduino/filtrage2/filtrage2.html

Mais la vulgarisation n’est pas simple, car il y a de nombreuses méthodes pour déterminer les coefficients d’un filtre numérique à partir d’un gabarit désiré ou à partir d’un filtre analogique.
Evidemment, il y a des softs qui déterminent les coefficients assez facilement tel que :
Matlab filterBuilder pour Mathcad Plus avec des fonctions spécifiques
des soft en C Iowa Hills Free Software Download Page
Je n’ai pas réussi à faire fonctionner ce lien
jayduino - Arduino stuff written by me

Pas mieux avec ce lien ou il manque un help
http://t-filter.engineerjs.com/

Or, le traitement numérique est crucial pour isoler un signal du bruit, avoir des valeurs moyennes mais en gardant une certaine dynamique d’un capteur…..

Mais, les résultats du filtre vont dépendre de la précision du calcul, de la période d’échantillonnage, donc du temps de calculs du processeur et du compilateur.
Souvent, les étudiants apprennent l’étude des filtres avec Matlab qui peut travailler à des fréquences d’échantillonnages très grandes et avec des précisions de calculs à plus de 18 chiffres après la virgule.

Mais avec un Arduino nano, on peut facilement observer les limites du traitement numérique d’un petit processeur ?

Avec l’IDE d’Arduino, la déclaration en float, n’a que 7 chiffres significatifs par excès :
dont voici un exemple : 0.1234567 ou 1.123456 ou 12.345

La base d’un programme sur les filtres numériques est une lecture d’un signal d’entrée analogique fournit par un GBF et une sortie qui peut être le terminal, le traceur série, ou une PWM. L’afficheur LCD n’est pas utile mais permet de debugger le programme parfois.
Un Arduino nano sera utilisé avec une PWM avec une fréquence de 32kHz bien supérieure à la période d’échantillonnage fixé arbitrairement à 1ms, et un RC de (10kohms et une capacité de 33nF) donc une fréquence de coupure de 482Hz atténuera l'ondulation du signal carré de la PWM comme on peut l’observer sur la figure suivante :

Le programme de base avec un gain unitaire est le suivant avec Te=1ms dans une routine interruption :
Ce programme de base dure 150microseconde, donc cela laisse du temps à faire du traitement numérique dans cette routine.

#include <LiquidCrystal.h>
#include <SoftwareSerial.h>
#include <TimerOne.h>


#define PWM3   3      //   timer2    
#define LED13    13       
LiquidCrystal lcd(9, 8, 4, 5, 6, 7);   // LiquidCrystal lcd(rs, en, d4, d5, d6, d7);
// Configuration des variables


unsigned    int    temps=0; 

            float  entree=0;
            float  entree1=0;
            float  entree2=0;
         
            float  sortie=0;
            float  sortie1=0;
            float  sortie2=0;
            
            float  out=0;
            const float b1 =2;
            const float b2 =1;
            
            const float a1 =-1.911;
            const float a2 =0.915;
            
            const float gain =1000;   
                    

void setup() {
  pinMode(LED13, OUTPUT);
  pinMode(PWM3,OUTPUT);

  Timer1.initialize(1000);           // initialize timer1, and set a 0,1 second period =>  100 000  pour 0.01s  10 000
  Timer1.attachInterrupt(callback);   // attaches callback() as a timer overflow interrupt
  lcd.begin(16, 2);                   //modifier pour un afficheur 20x4
  Serial.begin(9600); 

  TCCR2B = (TCCR2B & 0b11111000) | 0x01;         //pin 3  32khz    http://playground.arduino.cc/Main/TimerPWMCheatsheet
  
}



// Interruptions  tous les 1ms fait par le timer1***********************************
void callback()  {

digitalWrite(LED13,HIGH);  //permet de mesurer à l'oscillo, le temps du calcul du filtre et le temps de la routine d'interruption
entree=analogRead(A0); //convertisseur 10 bits sous 5V
   
      
/* ---- exemple filtre passe pas Butterwoth   pour fc=10hz,  fechantillon=1000Hz------         

entree2=entree1;      //entree(n-2)
entree1=entree;       //entree(n-1)

sortie2=sortie1;      //sortie(n-2)
sortie1=sortie;        //sortie(n-1)

//sortie=entree+entree1*b1+entree2*b2+entree3*b3-sortie1*a1-sortie2*a2-sortie3*a3 ;   //filtre passe pas recursif ordre 3
//sortie=(entree+entree1*b1+entree2*b2-sortie1*a1-sortie2*a2) ;                         //filtre passe pas recursif ordre 2
*/

sortie=entree/4  ;                        //mise à l'echelle 10 bits en entrée et 8 bits en sortie
//out = map(out, 0, 1023, 0, 255);   //mise à l'echelle 10 bits en entrée et 8 bits en sortie
if ( sortie>254)  {sortie=254;}
analogWrite(PWM3,sortie);          //


digitalWrite(LED13,LOW);

 
/*   creation du fichier CSV  traçage dans Excel
  Serial.print(entree);  Serial.print(";");  Serial.print("\t");   //entrée
  Serial.print(sortie);  Serial.print("\t");  Serial.println(";");  //mise à la ligne dans le terminal
*/    

}//fin routine



///////////////////////////////////////////// Boucle correspondant à la fonction main 
void loop() {  
  lcd.setCursor(0,0); 
  lcd.print("filtre numerique   ");

  lcd.setCursor(0,1); 
  lcd.print("IUT GEIIsoissons   "); 
} // fin loop

Le fichier de simulation ISIS peut être téléchargé sur ce lien

Attention la PWM demande beaucoup de ressource de calcul pour ISIS par le PC, donc souvent il vaut mieux travailler en reel qu’en simulation.

Avec une période d’échantillonnage de 1kHz, un signal sinusoïdale de 200Hz a été mis, on peut observer les différences de sorties à cause qu’il n’y a que 5 échantillons de mesure de ce signal.

Photo du banc en reel à venir et

Sommaire
Voici la définition des filtres numeriques

il existe 2 types de filtres
• RIF : Filtre à réponse impulsionnelle finie dépendant que de l’entrée et peu seulement converger.

RII : (Filtre à réponse impulsionnelle infinie) dépendant de l’entrée et d’ancienne et sortie et peut diverger

Dans un premier temps, Nous allons présenter comment on troupe les coefficients de filtres récursif RII.
Puis, la courbe d’atténuation théorique sera tracé, et 5 valeurs en réelles seront prix pour vérifier la théorie.

  • Etude filtre passe bas de premier ordre RC,
    En mode sinusoïdale et en mode temporel avec un signal carré.
  • Filtre butterworth passe bas, d’ordre 1
    Filtre de Butterworth — Wikipédia
    En mode sinusoïdale et en mode temporel avec un signal carré.
  • Filtre butterworth passe bas, d’ordre 2
  • Filtre butterworth passe bas, d’ordre 3

Dans un deuxième temps, des Filtres numériques moyenneur RIF seront presentés
https://www.f-legrand.fr/scidoc/docimg/numerique/filtre/rif/rif.html
http://herve.boeglen.free.fr/Tsignal/chapitre3/chapitre3.htm

De nombreux essais vont etre realisés sur ce sujet pour expliquer le filtrage numerique par des etudiants.

Ce qui sera peut-être plus explicite qu’un spécialiste du domaine ou que des bibliothèques ou l’on ne connait pas les tenants et les aboutissants
https://playground.arduino.cc/Code/Filters/

Etant donné que le temps de la lecture de l’entrée analogique et de son écriture dure 150us. La routine d’interruption est imposée à 1ms par conséquent le filtre analogique RC
est à 482Hz qui permet d’avoir la valeur moyenne de la PWM de sortie imposée à 32Khz

Explication : On va commencer par injecter une tension quelconque inférieure à 5V en entrée analogique A0 de l’arduino
On visualisera en sortie de l’arduino nano D3 c’est-à-dire en entrée du filtre analogique RC le signal de la tension désirée. Puis on récupéra le signal de sortie pas vraiment nécessaire dans notre cas sur l’oscilloscope.

Voici les résultats observés :

La fréquence de coupure a 482Hz car fc=1/(2πRC)=1/(2π10^3〖33*10〗^(-9) )

Le programme de base avec un gain unitaire est le suivant avec Te=1ms dans une routine interruption :

Nous avons tout d’abord procédé à une série de simulation sur ISIS pour vérifier certains points du programme avant de téléverser le programme. Voici les différentes composantes de la simulation Isis. Vous trouverez le lien de téléchargement au dessus.

Il aurait été possible que l’Arduino lise un signal carré extérieur et de faire une sortie analogique des valeurs sur une PWM pour démontrer l’utilisation d’un filtre numérique sur un oscilloscope avec une routine d’interruption de 1ms pour faire les calculs du filtrage.
Mais, nous avons utilisé juste le terminal pour observer les calculs du filtre et vérifier que le float avec ces 7 chiffres significatif fonctionne bien comme on peut l’observer sur le code suivant qui ont été placé dans Excel pour observer les dynamiques du filtre.

pour faire un filtre RC du premier ordre, il faut un seul terme a1 correspondant à l’equation
sortie(z)=(entree*(1-a1))/(1+a1*z-1) avec a1=e-Te/taux
cela se concretisite par le programme suivant

#include <LiquidCrystal.h>
#include <SoftwareSerial.h>
#include <TimerOne.h>
#include <avr/wdt.h>   //chien de garde

LiquidCrystal lcd(9, 8, 4, 5, 6, 7);   // LiquidCrystal lcd(rs, en, d4, d5, d6, d7);
// Configuration des variables

unsigned    int     temps=0;
            float  entree=0;
            float  sortie=0;

// Filtre numerique premier ordre
void setup() {
    
  Timer1.initialize(1000);           // initialize timer1, and set a 0,1 second period =>  100 000  pour 0.01s  10 000
  Timer1.attachInterrupt(callback);   // attaches callback() as a timer overflow interrupt
  lcd.begin(16, 2);                   //modifier pour un afficheur 20x4
  Serial.begin(9600);

  TCCR2B = (TCCR2B & 0b11111000) | 0x01;         //pin 10  32khz    http://playground.arduino.cc/Main/TimerPWMCheatsheet
                                            
}
// Interruptions  tous les 0.001s fait par le timer1***********************************
void callback()  {
temps++;
sortie=entree*0.008+sortie*0.992;   //filtre passe pas    a1=0.992 pour une constante de temps de 
  Serial.print(temps);
  Serial.print(";");
  Serial.print("\t");
    
  Serial.print(entree);
  Serial.print(";");
  Serial.print("\t");
  
  Serial.print(sortie,2);
  Serial.print("\t");
  Serial.println(";");  //mise à la ligne
}//fin routine

///////////////////////////////////////////// Boucle correspondant à la fonction main
void loop() {  

if (temps>=50)         //1=singnal carré de 100ms
{
temps=0;  //reset : attention cette boucle doit etre infereir à son temps d'execution
if (entree==16384) {entree=0;}  else {entree=16384;}  
  
}//temps>100

  
} // fin loop

pour une constante de temps deisirée de 10ms et un échantillonnage à 1ms, le terme a1=0.91.
Sur la figure suivante, on peut observer que la constante de temps est environ celle théorique de 10ms, pour une amplitude de 1023 de l’entrée (10bits).

Nous allons mettre ce programme en pratique et vérifier que nous avons réellement 10ms pour la constante de temps

0.63 x 5V = 3.15 V
Sur la courbe on voit bien que le temps correspondant à 3.5V est ≈ 10ms

Avec le terme 0.98 du filtre échantillon à 1ms, On peut observer que la constante de temps est légèrement inférieure à la valeur théorique de 50ms, pour une amplitude de 16384 (14bits) de l’entrée.
Lorsque la constante de temps devient très grande devant la période du signal, la sortie du filtre donnera une oscillation avec une amplitude correspondant à l’équation suivant autour de la valeur moyenne du signal d’entrée

En pratique pour une fréquence d’entrée de 10 Hz nous avons bien la constante de temps qui dure 50ms


Nous avons le temps d’une montée qui correspond à la constante de temps et qui dure 1divison c’est à dire 50ms

Avec le terme 0.998 du filtre échantillon à 1ms, la constante de temps theorique

on peut observer que l’amplitude

Donc, il y a des petites différences entre la théorie et la simulation mais globalement, cela fonctionne correctement en pratique.

Remarque : ne jamais depasser la valeur de 1 pour a1, sinon c’est instable (divergence), d’où l’obligation de declarer en float a1
faut il des bibliotheques filtre pour cela ? ? ? ? ? ?

Suite étude de filtre passe bas de premier ordre RC
(En mode temporel avec un signal carré)

Photo du banc en réel

Visualisation du signal d’entrée (carrée) sur le Channel 1 fournis par le GBF à 100Hz

Matériels utilisés :

GBF, arduino nano (attention au bootloader), résistance de 10kΩ, condensateur de 33nF, un oscilloscope.
Afficheur LCD I2C (facile à câbler dans notre cas)

On observe un signal en sortie ici

• Atténuer du au coefficient de division, ici nous avons diviser l’entrée par 5, contrairement à notre premier programme
• Déphasé du à la lecture de la pin de sortie de l’arduino
• Avec une légère courbure au début du signal ascendant dû à la charge du condensateur
Le mode temporel nous permet de voir l’importance du filtre avec l’atténuation constatée et l’élimination des bruits ou impureté du signal.

Nous avions choisi cependant utilisé l’afficheur LCD I2C (avec A0 , A1,A2 non shunté ) pour faciliter le câblage à inclure dans sa bibliothèque pour la compilation et extraction en binaires compilés sur ISIS pour simulation.

Remarque très importante : N’oublier pas de prendre en considération le LCD I2C avec cependant les modifications nécessaires

Voici notre programme fonctionnel avec la modification et prise en compte de l’afficheur lcd I2C

#include <LiquidCrystal_I2C.h>
#include <SoftwareSerial.h>
#include <TimerOne.h>
#include <avr/wdt.h>   //chien de garde
#include "LowPower.h"     //https://github.com/rocketscream/Low-Power
#include <Wire.h>

LiquidCrystal_I2C lcd(0x27, 16, 2);     //A0, A1, A2 non shunter

// Configuration des variables


unsigned    int     temps=0;
            float  entree=0;
            float  sortie=0;
unsigned    int inc=0;



// Filtre numerique premier ordre
void setup() {
 pinMode(3, OUTPUT);   //PWM
 pinMode(13, OUTPUT);   //led
 lcd.init();  //et pas lcd begin comme cetraine biblio  // activer l'affichage
 lcd.display();   // allumer retroeclairage
 lcd.backlight();

  Timer1.initialize(1000);           // initialize timer1, and set a 0,1 second period =>  100 000  pour 0.01s  10 000
  Timer1.attachInterrupt(callback);   // attaches callback() as a timer overflow interrupt
  lcd.begin(16, 2);                   //modifier pour un afficheur 20x4
  Serial.begin(9600);

  TCCR2B = (TCCR2B & 0b11111000) | 0x01;         //pin 10  32khz    http://playground.arduino.cc/Main/TimerPWMCheatsheet
                                            

//wdt_enable(WDTO_15MS);
}



// Interruptions  tous les 0.001s fait par le timer1***********************************
void callback()  {
//temps++;

//if ( digitalRead(13)== 1 ) {digitalWrite(13,LOW);}  else {digitalWrite(13,HIGH);}

digitalWrite(13,HIGH);
entree=analogRead(A0);      //10 bits
//sortie=entree*0.008+sortie*0.992;   //filtre passe bas
sortie=inc++;
//sortie=entree/4;
analogWrite(3, sortie);     // 8bits
digitalWrite(13,LOW);

  /*
  Serial.print(temps);
  Serial.print(";");
  Serial.print("\t");
    
  Serial.print(entree);
  Serial.print(";");
  Serial.print("\t");
  
  Serial.print(sortie,2);
  Serial.print("\t");
  Serial.println(";");  //mise à la ligne
*/

//wdt_reset();
}//fin routine



///////////////////////////////////////////// Boucle correspondant à la fonction main
void loop() {  

lcd.setCursor(0,0);   //colonne, ligne
lcd.print("FILTRE NUMERIQUE");
lcd.setCursor(0,1);   //colonne, ligne
lcd.print("IUT GEII SOISSONS");
//lcd.print(entree,1);

 

  
} // fin loop

Remarque : Lors de la compilation de ce programme, il peut s’avérer qu’il y ait des erreurs de compilation. Pas de panique !!
Sans oublié de télécharger (inclure) toutes les librairies nécessaires
Cela n’a pas de rapport avec le code. Arduino souvent chose incompréhensible préfère que l’on compile avec le old bootloader surtout avec l’arduino Nano.
Pour se faire allez sur : outils>processeur>old bootloader.
Et recompiler.

TEST : La routine d’interruption d’une durée de 1ms

Nous avons alors créé une routine d’interruption avec le callback qui doit « normalement » durer 1ms
Est-ce vraiment le cas ???
Vérifions cela alors en simulation

Résultats observés :

La courbe du channel D correspond au front montant et descendant de la routine d’interruption
On voit bien que chaque division dure 1ms donc :
High ===> pendant 1ms
Low ==> pendant 1ms
Le temps d’une montée dure alors 1ms de même que le temps d’une descente : la routine d’interruption est alors correcte. YES !!!

TEST : temps d’incrémentation à chaque fois que l’on rentre dans la routine d’interruption

Nous avons comme expliqué précédemment une routine d’interruption qui dure 1ms. Dans cette routine nous avons écrit un bout de code qui doit incrémenter la sortie 256 fois (0 à 255). Nous allons alors vérifier le temps que met cette incrémentation à s’exécuter et nous assurer que cela ne dépasse pas la durée de la routine d’interruption

On peut alors observer (figure suivante) sur le channel B (courbe bleue) un signal triangulaire qui représente la courbe d’incrémentation de 0 à 255 qui s’effectue sue 13 divisions de 20 ms chacune
donc :
13*20ms=0.26ms < 1ms Yes !!!

Nous avons créé un filtre passe bas RC, nous allons alors faire un programme suiveur c’est-à-dire que la sortie recopie l’entrée donc un gain de 1
Mais comme nous avons une entrée sur 10 bits (1024) et une sortie à 8 bits (256) nous allons diviser l’entrée par 4 ( 1024/256)

Résultats observés :

Courbe en jaune = entrée A0
Courbe en bleu = sortie
La courbe en bleu représentant la sortie est en retard par rapport à la courbe en jaune dû au laps de temps du traitement de l’Analog read

Lorsqu’on fait un zoom sur la sortie on s’aperçoit qu’il y’a un hachage dû a la PWM.

Cette mesure est dû au fait que nous voulons en sortie 8 bits.
Il faut donc diviser 1024 par 4 pour avoir 8 bits en sorties.

Filtre butterworth passe bas, d’ordre 1
(En mode sinusoïdale et en mode temporel avec un signal carré)

Le calcul mathématique d’un butterworth numérique quelque soit son ordre va être présenté.
Dans le lien suivant, un article présente la méthode bilinéaire mais avec un filtre de coupure de 1 rad/s ?? Ce qui n’aide pas à la compréhension.
https://www.dsprelated.com/showarticle/1119.php

Quelques soient l’ordre du filtre butterWorth, il n’y a pas de gain positif, le lien dans Wikipédia présente bien ce filtre analogique et il est très facile de déterminer la courbe d’atténuation avec Mathcad.
Nous allons prendre une fréquence de coupure de 10hz et une fréquence d’échantillonnage de 100Hz

Afin de vérifier, l’atténuation du filtre, la courbe de fréquence est tracée. On retrouve bien une atténuation de 0.707 ou de -3db à la fréquence de coupure.

Plus, l’écart entre la fréquence d’échantillonnage et la fréquence de coupure est grande, plus l’ordre N du filtre pourra être grand. mais plus le gain devra être re petit et la résolution des zeros devront être très précise.
pour bien tester le filtre numérique, il faut générer un signal sinusoïdal en entrée et voir en sortie l’atténuation pour chaque fréquence testée

sinon, en utilisant, la décomposition de série de fourriers d’un signal carré, L’amplitude de l’ondulation sera égale à l’équation suivante exemple avec une fréquence de coupure de 1hz.

Si on utilise, un ordre 1, alors l’ondulation théorique sera de 14%

Si on utilise un ordre 3, ondulation sera de 0.1%

Maintenant, il reste à implanter dans l’Arduino.
Evidement les coefficients a(z) sont en float car il faut une certaine précision.

Pour un filtre butteworth d’ordre 3, avec une fréquence d’échantillonnage de 100Hz et de fréquence de coupure de 10hz,
Si l’on choisit les coefficients au centième près par défaut, le gain passe 1/55 à 1/53 et la courbe d’atténuation fréquentielle sera légèrement modifiée avec une fréquence de coupure légèrement décalée

Le temps de calcul du filtre avec l’ATMEL 328 est de 120us pour un troisième ordre
J’aurai pensé que le temps de calcul aurait été bien plus long. le compilateur doit faire des simplifications

Voici les résultats simuler avec Matlab dans un premier temps avec une fréquence d’échantillon de 100Hz, un signal carré de 10hz, et une fréquence de coupure de 10hz donc avec 10 échantillons par période du signal, ce n’est pas terrible

Si l’on essaye de faire avec une fréquence d’échantillon de 1000hz, il y a une divergence du signal de sortie même avec 3 chiffres ou 4 chiffres après la virgule des coefficients zéros du filtre.
C’est là que l’on voit qu’entre la théorie et la pratique, il y a un monde
Plus l’ordre est élevé et plus s’il y a un écart entre fe et fc et plus il faudra une précision sur les coefficients du filtre.
Remarque : Lorsqu’il y a un trop grand écart entre le gain statique avec les coefficients pris par défaut ou avec des coefficients avec 7 chiffres significatifs alors le filtrage devient divergent.

1 Like

Il est possible de faire des exemples en mode fréquentielle mais avec une amplitude du signal d’entre unipolaire (donc que positif)
Pour 1 ordre passe bas avec pour Fe =1000 Hz et le terme a1=0.81 donc pour une frequence de coupure de 33Hz

Pour une fréquence d’entrée de 10Hz, sur l’oscilloscope on peut observer aucune atténuation du signal d’entrée car l’amplitude est de 10 carreaux en sortie divisé par l’amplitude de 10 carreaux entrée

Pour une fréquence d’entrée de 50Hz, l’atténuation du signal d’entrée est de 0.6

Pour une fréquence d’entrée de 100Hz, l’atténuation est de 0.3 du signal d’entrée.

Pour une fréquence d’entrée de 500Hz, l’atténuation est de 0.1 du signal d’entrée.
on peut observer qu'il n'y a que 2 echantillons en sortie à cettte frequence qui est celle de la frequence d'echantillonage divisé par 2.

donc de verifier la courbe d'attenuation du filtre et de verifier les calculs des coefficients du filtre....

Pour un Butterworth PASSE BAS de n=2éme ordre avec une période d’échantillonnage toujours de 1ms, une fréquence de coupure de 50Hz, les coefficients du filtres sont les suivants
bo=1, b1=2, b2=1, a0=1, a1=-1.56, a2=0,64, gain=1/50

Les courbes dans Matlab sont les suivantes pour 10Hz, l’atténuation est de 1

Pour 50Hz, l’atténuation est de 0.707 correspondants bien à la fréquence de coupure

Pour 100Hz, l’atténuation est de 0.25 correspondant bien à l’équation suivante
Attenuation=[(Frequence coupure)/(frequence entree)]^n=(50/100)^2=0.25

Evidement avec l’arduino et les memes coeficient que precedent, les données sont identiques
pour 10Hz, l’atténuation est aussi de 1

Pour 50Hz, l’atténuation est de 0.707 correspondants bien à la fréquence de coupure

Pour 100Hz, l’atténuation est aussi de 0.25

#include <Wire.h>
//#include <LiquidCrystal_I2C.h>
#include <LiquidCrystal.h>
#include <SoftwareSerial.h>
#include <TimerOne.h>


#define PWM3   3      //   timer2    
#define LED13    13       
LiquidCrystal lcd(9, 8, 4, 5, 6, 7);   // LiquidCrystal lcd(rs, en, d4, d5, d6, d7);
// Configuration des variables
//LiquidCrystal_I2C lcd(0x27, 20, 4);     //A0, A1, A2 non shunter 20x4ligne


unsigned    int    temps=0; 

            float  entree=0;
            float  entree1=0;
            float  entree2=0;
         
            float  sortie=0;
            float  sortie1=0;
            float  sortie2=0;
            
            float  out=0;
            const float b1 =2;
            const float b2 =1;
            
            const float a1 =-1.56;
            const float a2 =0.64;
            

                    

void setup() {
  pinMode(LED13, OUTPUT);
  pinMode(PWM3,OUTPUT);

  Timer1.initialize(1000);           // initialize timer1, and set a 0,1 second period =>  100 000  pour 0.01s  10 000
  Timer1.attachInterrupt(callback);   // attaches callback() as a timer overflow interrupt
  lcd.begin(20, 4);                   //modifier pour un afficheur 20x4
/*
lcd.init();               //et pas lcd begin comme cetraine biblio
lcd.display();            // activer l'affichage
lcd.backlight();          // allumer retroeclairage  
*/  
  
  Serial.begin(9600); 

  TCCR2B = (TCCR2B & 0b11111000) | 0x01;         //pin 3  32khz    http://playground.arduino.cc/Main/TimerPWMCheatsheet
  
}



// Interruptions  tous les 1ms fait par le timer1***********************************
void callback()  {

digitalWrite(LED13,HIGH);  //permet de mesurer à l'oscillo, le temps du calcul du filtre et le temps de la routine d'interruption
entree=analogRead(A0); //convertisseur 10 bits sous 5V
  
      
// ---- exemple filtre passe pas Butterwoth   pour fc=10hz,  fechantillon=1000Hz------         

entree2=entree1;      //entree(n-2)
entree1=entree;       //entree(n-1)

sortie2=sortie1;      //sortie(n-2)
sortie1=sortie;       //sortie(n-1)


//sortie=entree+entree1*b1+entree2*b2+entree3*b3-sortie1*a1-sortie2*a2-sortie3*a3 ;   //filtre passe pas recursif ordre 3
sortie=(entree+entree1*b1+entree2*b2-sortie1*a1-sortie2*a2) ;                        //filtre passe pas recursif ordre 2
//sortie=(entree*0.19+sortie*0.81);                                 //filtre passe bas    a1=0.992 pour une constante de temps de 10ms  Te=1ms durre algo 150us                                                                                   //filtre passe pas recursif ordre 1
//sortie= entree  ;                                                                                     //gain unitaire
out=(sortie/1);
out=(out/4)  ;       //mise à l'echelle 10 bits en entrée et 8 bits en sortie                     
analogWrite(PWM3,out);          //    d


 
//   creation du fichier CSV  traçage dans Excel
//  Serial.print(entree);  Serial.print(";");  Serial.print("\t");   //entrée
//   Serial.print(sortie);  Serial.print("\t");  Serial.println(";");  //mise à la ligne dans le terminal
// 
digitalWrite(LED13,LOW);  

}//fin routine



///////////////////////////////////////////// Boucle correspondant à la fonction main 
void loop() {  
  lcd.setCursor(0,0); 
  lcd.print("Te=1ms");
  lcd.setCursor(0,1); 
  lcd.print("    ");
      lcd.setCursor(0,2);    
  lcd.print("   ");
    lcd.setCursor(0,3);  
//   lcd.print("Gain=1/50");
} // fin loop

La fréquence d’entrée ne peut être supérieure à fréquence d’échantillonnage divisé par 2, sinon, il y a des duplications spectrales du signal de sortie.
D’ailleurs l’amplitude de la sortie ne va plus diminuer mais augmenter au-delà d’une fréquence d’entrée supérieure à 500Hz dans le cas précèdent

Quelle serait la fréquence d’échantillon maximale ?
Etant donné que la conversion et des temps de calculs du filtre dure 0.2ms, la fréquence d’échantillon pourrait être de 4000hz

Pour un Butterworth passe bas de n=2éme ordre avec une période d’échantillonnage toujours de 0.25ms,
Donc timer1 à 250 dans le programme et plus 1000
Pour une fréquence de coupure de 50Hz, les coefficients du filtres sont les suivants
bo=1, b1=2, b2=1, a0=1, a1=-1.89, b2=0,895, gain=1/800

Pour 50Hz, avec cette nouvelle période et les coefficients, l’atténuation est d 0.707 correspondants bien à la fréquence de coupure

Pour 100Hz, l’atténuation est un peu inférieure à 0.25 correspondant à 0.2

Pour 1000Hz, l’atténuation est de 0.02=0.05V/2.5V et correspond bien la théorie
Atténuation=[(Frequence coupure)/(frequence entree)]^n=(50/1000)^2=0.025

A vous de jouer pour donner les coefficients avec toujours une période d’échantillonnage de 1ms.

  • pour une fréquence de coupure désirée de 10Hz pour un deuxième ordre puis faire les essais.
  • pour une fréquence de coupure désirée de 1Hz pour un deuxième ordre puis faire les essais.
  • pour une fréquence de coupure de 50Hz pour un troisième ordre puis faire les essais.

Pourrait-on demander à l’Arduino de faire le calcul des coefficients en rentrant juste la période d’échantillonnage, l’ordre du filtre, la fréquence de coupure désirée ???????…..

Étant donné qu’include math.h ne fait pas les calculs en complexe……
Voici le calcul en direct des coefficients pour un butterworth Deuxième ordre à partir de la fréquence d’échantillonnage et de la fréquence de coupure.
Il a fallu juste développer et simplifier les équations pour avoir a1 et A2

Voici les valeurs pour fréquence d’échantillon de 1000Hz et fréquence de coupure de 50hz comme précédent

Si l’on désire 1Hz, voilà les coefficients déterminé directement par le microcontrôleur.
malgre, les petits valeurs, le micro fait bien le filtre

Voici le code

#include <Wire.h>
//#include <LiquidCrystal_I2C.h>
#include <LiquidCrystal.h>
#include <SoftwareSerial.h>
#include <TimerOne.h>
#include <math.h>


#define PWM3   3      //   timer2    
#define LED13    13       
LiquidCrystal lcd(9, 8, 4, 5, 6, 7);   // LiquidCrystal lcd(rs, en, d4, d5, d6, d7);
// Configuration des variables
//LiquidCrystal_I2C lcd(0x27, 20, 4);     //A0, A1, A2 non shunter 20x4ligne


unsigned    int    temps=0; 

            float  entree=0;
            float  entree1=0;
            float  entree2=0;
         
            float  sortie=0;
            float  sortie1=0;
            float  sortie2=0;
            
            float  out=0;
            const float b1 =2;
            const float b2 =1;
            
            float a1=0;
            float a2=0;
            float denominateur=0;
            float Gain=0;

            float  fe=1000;  //frequence d'echantilonnage   
            float  fc=1;  //frequence de coupure desiréé     
            

                    

void setup() {
  pinMode(LED13, OUTPUT);
  pinMode(PWM3,OUTPUT);

  Timer1.initialize(1000);           // initialize timer1, and set a 0,1 second period =>  100 000  pour 0.01s  10 000
  Timer1.attachInterrupt(callback);   // attaches callback() as a timer overflow interrupt
  lcd.begin(20, 4);                   //modifier pour un afficheur 20x4
/*
lcd.init();               //et pas lcd begin comme cetraine biblio
lcd.display();            // activer l'affichage
lcd.backlight();          // allumer retroeclairage  
*/  
  
  Serial.begin(9600); 

  TCCR2B = (TCCR2B & 0b11111000) | 0x01;         //pin 3  32khz    http://playground.arduino.cc/Main/TimerPWMCheatsheet

denominateur=(fe*fe+fe*fc*PI*sqrt(2)+fc*fc*PI*PI);
a1=(fc*fc*PI*PI-fe*fe)*2/denominateur;
a2=(fe*fe*fe*fe+(fc*fc*fc*fc*PI*PI*PI*PI))/(denominateur*denominateur)         ;                             //pow(base, exponent)
Gain=(1+2+1)/(1+a1+a2);
  
}



// Interruptions  tous les 1ms fait par le timer1    fe=1000Hz***********************************
void callback()  {
if ( digitalRead(13)== 1 ) {digitalWrite(13,LOW);}  else {digitalWrite(13,HIGH);}
//digitalWrite(LED13,HIGH);  //permet de mesurer à l'oscillo, le temps du calcul du filtre et le temps de la routine d'interruption
entree=analogRead(A0); //convertisseur 10 bits sous 5V
  
      
// ---- exemple filtre passe pas Butterwoth   pour fc=10hz,  fechantillon=1000Hz------         

entree2=entree1;      //entree(n-2)
entree1=entree;       //entree(n-1)

sortie2=sortie1;      //sortie(n-2)
sortie1=sortie;       //sortie(n-1)


//sortie=entree+entree1*b1+entree2*b2+entree3*b3-sortie1*a1-sortie2*a2-sortie3*a3 ;   //filtre passe pas recursif ordre 3
sortie=(entree+entree1*b1+entree2*b2-sortie1*a1-sortie2*a2) ;                        //filtre passe pas recursif ordre 2
//sortie=(entree*0.19+sortie*0.81);                                 //filtre passe bas    a1=0.992 pour une constante de temps de 10ms  Te=1ms durre algo 150us                                                                                   //filtre passe pas recursif ordre 1
//sortie= entree  ;                                                                                     //gain unitaire
out=(sortie/Gain);
out=(out/4)  ;       //mise à l'echelle 10 bits en entrée et 8 bits en sortie                     
analogWrite(PWM3,out);          //    d

//digitalWrite(LED13,LOW);  

}//fin routine



///////////////////////////////////////////// Boucle correspondant à la fonction main 
void loop() {  
  lcd.setCursor(0,0); 
  lcd.print("fe=");
  lcd.print(fe,0);
    lcd.print(" fc=");
  lcd.print(fc,0);
  lcd.setCursor(0,1); 
  lcd.print("b0=1 b1=2 b2=1 a0=1");
  lcd.setCursor(0,2);    
  lcd.print("a1=");
  lcd.print(a1,3);
   lcd.print(" a2=");
   lcd.print(a2,3); 
   lcd.setCursor(0,3);  
   lcd.print("Gain=1/");
   lcd.print(Gain,1);



} // fin loop

c’est simple les math…il reste à verifier pour 10Hz, le deuxieme ordre et faire le troisieme ordre

Pour un Butterworth PASSE BAS de n=3éme ordre avec une période d’échantillonnage toujours de 1ms, une fréquence de coupure de 50Hz, les coefficients du filtres sont les suivants
bo=1, b1=3, b2=3, b3=1, a0=1, a2=-2.374, a2=1.93, a3=-0.53, gain=1/340

Les courbes en simulation dans ISIS sont les suivantes pour 10Hz, l’atténuation est de 1

Pour 50Hz, l’atténuation est de 0.707 correspondants bien à la fréquence de coupure

Pour 100Hz, l’atténuation est de bien de 0.125 correspondant bien à l’équation suivante
Attenuation=[(Frequence entree)/(frequence coupure)]^n=(50/100)^3=0.125

Le temps de calcul pour faire un 1er ordre est 0.18ms
Le temps de calcul pour faire un 2er ordre est de 0.2ms
Le temps de calcul pour faire un 3er ordre est de 0.25ms

Filtre butterworth passe bas, d’ordre 3

Effectivement include math.h ne fait pas les calculs en complexe……

C’est bien dommage du fait que la variable float s’arrête qu’à sept chiffres après la virgule. Du coup cela devient difficile pour le compilateur d’effectuer de gros calcul.

Il existe donc plusieurs façons de trouver les coefficients du filtre du troisième ordre soit en utilisant les logiciels de calcul prédéfinis , soit en développant les calculs nous-mêmes et l’intégrant sur le compilateur arduino…

On peut aussi utiliser une fonction sur matlab dans le but de trouver les différents coefficients a et b

Voici la fonction qui permet de calculer des coefficients sur matlab pour un butterworth troisième ordre à partir de la fréquence d’échantillonnage et de la fréquence de coupure.

% N= filter order
% fc= -3 dB frequency in Hz
% fs= sample frequency in Hz
% b = numerator coefficients of digital filter
% a = denominator coefficients of digital filter


function [b,a]= butter_synth(N,fc,fs);
%
if fc>=fs/2;
   error('fc must be less than fs/2')
end
% I.  Find poles of analog filter
k= 1:N;
theta= (2*k -1)*pi/(2*N);
pa= -sin(theta) + j*cos(theta);     % poles of filter with cutoff = 1 rad/s
%
% II.  scale poles in frequency
Fc= fs/pi * tan(pi*fc/fs);          % continuous pre-warped frequency
pa= pa*2*pi*Fc;                     % scale poles by 2*pi*Fc
%
% III.  Find coeffs of digital filter
% poles and zeros in the z plane
p= (1 + pa/(2*fs))./(1 - pa/(2*fs));      % poles by bilinear transform
q= -ones(1,N);                   % zeros
%
% convert poles and zeros to polynomial coeffs
a= poly(p);                   % convert poles to polynomial coeffs a
a= real(a);
b= poly(q);                   % convert zeros to polynomial coeffs b
K= sum(a)/sum(b);       % amplitude scale factor
b= K*b;

et voici le programme exécutable

N= 3           % ordre  du filtre
fc= 50;        % Hz freq de coupure
fe= 1000;       % Hz sample freq d'echantillonage te=1ms
[b,a]= butter_synth(N,fc,fe)

%calcul du gain
z=1;
a= 1+3*(z^-1)+3*(z^-2)+(z^-3) 
b=1-2.3741+1.9294-0.5321
gain=a/b

Etant donné que c’est un ordre 3, les coefficients de b sont pris par rapport au coef du triangle de pascal
D’où b0=1, b1=3, b2=3, b3=1.

Voici une capture des résultats trouvés après compilation SUR Matlab

a0=1, a1=-2.3741, a2=1.9294 , a3=-0.5321 et Gain=344

Il ne reste plus qu’à compiler sur arduino

voici le code

#include <Wire.h>
#include <LiquidCrystal_I2C.h>
//#include <LiquidCrystal.h>
#include <SoftwareSerial.h>
#include <TimerOne.h>
#include <math.h>


#define PWM3   3      //   timer2    
#define LED13    13      
//LiquidCrystal lcd(9, 8, 4, 5, 6, 7);   // LiquidCrystal lcd(rs, en, d4, d5, d6, d7);

//LiquidCrystal_I2C lcd(0x27, 16, 2);     //A0, A1, A2 non shunter
// Configuration des variables
LiquidCrystal_I2C lcd(0x27, 20, 4);     //A0, A1, A2 non shunter 20x4ligne


unsigned    int    temps=0;

            float  entree=0;
            float  entree1=0;
            float  entree2=0;
            float  entree3=0;
        
            float  sortie=0;
            float  sortie1=0;
            float  sortie2=0;
            float  sortie3=0;
            
            float  out=0;
            const float b1 =3;
            const float b2 =3;
            const float b3 =1;
            
            // valeur de a trouver sur matlab
            const float a1 =-2.3741;
            const float a2 =1.9294;
            const float a3 =-0.5321;
            float Gain=0;
            
            float  fe=1000;  //frequence d'echantilonnage  
            float  fc=100;  //frequence de coupure desiréé    
            

                    

void setup() {
  pinMode(LED13, OUTPUT);
  pinMode(PWM3,OUTPUT);

  Timer1.initialize(1000);           // initialize timer1, and set a 0,1 second period =>  100 000  pour 0.01s  10 000
  Timer1.attachInterrupt(callback);   // attaches callback() as a timer overflow interrupt
  lcd.begin(20, 4);                   //modifier pour un afficheur 20x4
/*
lcd.init();               //et pas lcd begin comme cetraine biblio
lcd.display();            // activer l'affichage
lcd.backlight();          // allumer retroeclairage  
*/  
  
  Serial.begin(9600);

  TCCR2B = (TCCR2B & 0b11111000) | 0x01;         //pin 3  32khz    http://playground.arduino.cc/Main/TimerPWMCheatsheet

Gain=(1+3+3+1)/(1+a1+a2+a3);
  
}



// Interruptions  tous les 1ms fait par le timer1    fe=1000Hz***********************************
void callback()  {
if ( digitalRead(13)== 1 ) {digitalWrite(13,LOW);}  else {digitalWrite(13,HIGH);}
//digitalWrite(LED13,HIGH);  //permet de mesurer à l'oscillo, le temps du calcul du filtre et le temps de la routine d'interruption
entree=analogRead(A0); //convertisseur 10 bits sous 5V
  
      
// ---- exemple filtre passe pas Butterwoth   pour fc=10hz,  fechantillon=1000Hz------        

entree3=entree2;      //entree(n-3)
entree2=entree1;      //entree(n-2)
entree1=entree;       //entree(n-1)

sortie3=sortie2;      //sortie(n-3)
sortie2=sortie1;      //sortie(n-2)
sortie1=sortie;       //sortie(n-1)


sortie=entree+entree1*b1+entree2*b2+entree3*b3-sortie1*a1-sortie2*a2-sortie3*a3 ;   //filtre passe pas recursif ordre 3

out=(sortie/Gain);
out=(out/4)  ;       //mise à l'echelle 10 bits en entrée et 8 bits en sortie                    
analogWrite(PWM3,out);          //    d

//digitalWrite(LED13,LOW);  

}//fin routine



///////////////////////////////////////////// Boucle correspondant à la fonction main
void loop() {  
  lcd.setCursor(0,0);
  lcd.print("fe=");
  lcd.print(fe,0);
    lcd.print(" fc=");
  lcd.print(fc,0);
  lcd.print(" b0=1");
  lcd.setCursor(0,1);
  lcd.print("b1=3 b2=3 b3=1 a0=1");
  lcd.setCursor(0,2);    
  lcd.print("a1=");
  lcd.print(a1,3);
   lcd.print(" a2=");
   lcd.print(a2,3);
   lcd.setCursor(0,3); 
   lcd.print("a3="); 
   lcd.print(a3,3);
   lcd.print(" Gain=1/");
   lcd.print(Gain,1);



} // fin loop

Pour une fréquence de coupure de Fc=10Hz et une période d’échantillonnage de 1ms
On a les résultats de simulations suivantes:
5V en sortie sur 5V en entrée … une atténuation de 1

Pour une fréquence de coupure de Fc=50Hz
3,5V/5V=0,7 en atténuation

Pour une fréquence de coupure de Fc=100Hz
0,5V/5V=0,1 en atténuation

Une deuxième méthode aussi est envisageable

Deuxième méthode

On obtient les mêmes résultats

voici le code avec les différents paramètres de a(1,2,3,4)

Cette méthode est fastidieuse mais elle nous permet après de trouver les coefficients de a sans passer par les calculs intermédiaires. Je dirai qu’elle est généralisée pour l’ordre 3 du filtre.

#include <Wire.h>
#include <LiquidCrystal_I2C.h>
//#include <LiquidCrystal.h>
#include <SoftwareSerial.h>
#include <TimerOne.h>
#include <math.h>


#define PWM3   3      //   timer2    
#define LED13    13      
//LiquidCrystal lcd(9, 8, 4, 5, 6, 7);   // LiquidCrystal lcd(rs, en, d4, d5, d6, d7);

//LiquidCrystal_I2C lcd(0x27, 16, 2);     //A0, A1, A2 non shunter
// Configuration des variables
LiquidCrystal_I2C lcd(0x27, 20, 4);     //A0, A1, A2 non shunter 20x4ligne


unsigned    int    temps=0;

            float  entree=0;
            float  entree1=0;
            float  entree2=0;
            float  entree3=0;
        
            float  sortie=0;
            float  sortie1=0;
            float  sortie2=0;
            float  sortie3=0;
            
            float  out=0;
            float denominateur=0;
            const float b1 =3;
            const float b2 =3;
            const float b3 =1;
            
            float Gain=0;
             float a1 =0;
             float a2 =0;
             float a3 =0;
            
            float  fe=1000;  //frequence d'echantilonnage  
            float  fc=50;  //frequence de coupure desiréé    
       

void setup() {
  pinMode(LED13, OUTPUT);
  pinMode(PWM3,OUTPUT);

  Timer1.initialize(1000);           // initialize timer1, and set a 0,1 second period =>  100 000  pour 0.01s  10 000
  Timer1.attachInterrupt(callback);   // attaches callback() as a timer overflow interrupt
  lcd.begin(20, 4);                   //modifier pour un afficheur 20x4
/*
lcd.init();               //et pas lcd begin comme cetraine biblio
lcd.display();            // activer l'affichage
lcd.backlight();          // allumer retroeclairage  
*/  
  
  Serial.begin(9600);

  TCCR2B = (TCCR2B & 0b11111000) | 0x01;         //pin 3  32khz    http://playground.arduino.cc/Main/TimerPWMCheatsheet

denominateur=(fe*fe+fe*fc*PI+fc*fc*PI*PI);
a1=-((fe*fe-fc*fc*PI*PI)*2/denominateur)-(fe-fc*PI)/(fe+fc*PI);
a2=(((fe*fe-fc*fc*PI*PI)*2)/denominateur)*((fe-fc*PI)/(fe+fc*PI))+((fe*fe*fe*fe+fc*fc*fe*fe*PI*PI)+(fc*fc*fc*fc*PI*PI*PI*PI))/(denominateur*denominateur);
a3=-(((fe*fe*fe*fe+fc*fc*fe*fe*PI*PI)+(fc*fc*fc*fc*PI*PI*PI*PI))/(denominateur*denominateur))*((fe-fc*PI)/(fe+fc*PI));
Gain=(1+3+3+1)/(1+a1+a2+a3);
  
}



// Interruptions  tous les 1ms fait par le timer1    fe=1000Hz***********************************
void callback()  {
if ( digitalRead(13)== 1 ) {digitalWrite(13,LOW);}  else {digitalWrite(13,HIGH);}
//digitalWrite(LED13,HIGH);  //permet de mesurer à l'oscillo, le temps du calcul du filtre et le temps de la routine d'interruption
entree=analogRead(A0); //convertisseur 10 bits sous 5V
  
      
// ---- exemple filtre passe pas Butterwoth   pour fc=10hz,  fechantillon=1000Hz------        

entree3=entree2;      //entree(n-3)
entree2=entree1;      //entree(n-2)
entree1=entree;       //entree(n-1)

sortie3=sortie2;      //sortie(n-3)
sortie2=sortie1;      //sortie(n-2)
sortie1=sortie;       //sortie(n-1)


sortie=entree+entree1*b1+entree2*b2+entree3*b3-sortie1*a1-sortie2*a2-sortie3*a3 ;   //filtre passe pas recursif ordre 3

out=(sortie/Gain);
out=(out/4)  ;       //mise à l'echelle 10 bits en entrée et 8 bits en sortie                    
analogWrite(PWM3,out);          //    d



}//fin routine



///////////////////////////////////////////// Boucle correspondant à la fonction main
void loop() {  
  lcd.setCursor(0,0);
  lcd.print("fe=");
  lcd.print(fe,0);
  lcd.print(" fc=");
  lcd.print(fc,0);
  lcd.print(" b0=1");
  lcd.setCursor(0,1);
  lcd.print("b1=3 b2=3 b3=1 a0=1");
  lcd.setCursor(0,2);    
  lcd.print("a1=");
  lcd.print(a1,3);
  lcd.print(" a2=");
  lcd.print(a2,3);
  lcd.setCursor(0,3); 
  lcd.print("a3="); 
  lcd.print(a3,3);
  lcd.print(" G=1/");
  lcd.print(Gain,1);

} // fin loop

Pour une fréquence de coupure de Fc=100Hz
0,5V/5V=0,1 en atténuation

Les résultats sont identiques.

Nous allons faire dans cette partie le calcul des coefficients du filtre Butterworth avec des fréquences de coupures différentes : 10Hz et 1Hz
Dans la partie présentation du filtre Butterworth c’est Mathcad qui a été utilisé pour déterminer les valeurs des coefficients du filtre selon l’ordre

Cependant à défaut d’avoir ce logiciel le programme pour le calcul des coefficients est donné avec Arduino précédemment donc nous allons juste changer les fréquences de coupure désirées puis faire les essais avec différentes fréquences d’entrée (qui ne doivent pas être supérieures à Fe/2 (Shannon) )

Rappel du code :

#include <Wire.h>
#include <LiquidCrystal_I2C.h>
//#include <LiquidCrystal.h>
#include <SoftwareSerial.h>
#include <TimerOne.h>
#include <math.h>


#define PWM3   3      //   timer2    
#define LED13    13     
LiquidCrystal_I2C lcd(0x27, 16, 2); 
//LiquidCrystal_I2C lcd(9, 8, 4, 5, 6, 7);   // LiquidCrystal lcd(rs, en, d4, d5, d6, d7);
// Configuration des variables
//LiquidCrystal_I2C lcd(0x27, 20, 4);     //A0, A1, A2 non shunter 20x4ligne


unsigned    int    temps=0;

            float  entree=0;
            float  entree1=0;
            float  entree2=0;
        
            float  sortie=0;
            float  sortie1=0;
            float  sortie2=0;
            
            float  out=0;
            const float b1 =2;
            const float b2 =1;
            
            float a1=0;
            float a2=0;
            float denominateur=0;
            float Gain=0;

            float  fe=1000;  //frequence d'echantilonnage  
            float  fc=10;  //frequence de coupure desiréé    
            

                    

void setup() {
  pinMode(LED13, OUTPUT);
  pinMode(PWM3,OUTPUT);

  Timer1.initialize(1000);           // initialize timer1, and set a 0,1 second period =>  100 000  pour 0.01s  10 000
  Timer1.attachInterrupt(callback);   // attaches callback() as a timer overflow interrupt
  lcd.begin(20, 4);                   //modifier pour un afficheur 20x4
/*
lcd.init();               //et pas lcd begin comme cetraine biblio
lcd.display();            // activer l'affichage
lcd.backlight();          // allumer retroeclairage  
*/  
  
  Serial.begin(9600);

  TCCR2B = (TCCR2B & 0b11111000) | 0x01;         //pin 3  32khz    http://playground.arduino.cc/Main/TimerPWMCheatsheet

denominateur=(fe*fe+fe*fc*PI*sqrt(2)+fc*fc*PI*PI);
a1=((fc*fc*PI*PI)-(fe*fe))*2/denominateur;
a2=(fe*fe*fe*fe+(fc*fc*fc*fc*PI*PI*PI*PI))/(denominateur*denominateur)         ;                             //pow(base, exponent)
Gain=(1+2+1)/(1+a1+a2);
  
}



// Interruptions  tous les 1ms fait par le timer1    fe=1000Hz***********************************
void callback()  {
if ( digitalRead(13)== 1 ) {digitalWrite(13,LOW);}  else {digitalWrite(13,HIGH);}
//digitalWrite(LED13,HIGH);  //permet de mesurer à l'oscillo, le temps du calcul du filtre et le temps de la routine d'interruption
entree=analogRead(A0); //convertisseur 10 bits sous 5V
  
      
// ---- exemple filtre passe pas Butterwoth   pour fc=10hz,  fechantillon=1000Hz------        

entree2=entree1;      //entree(n-2)
entree1=entree;       //entree(n-1)

sortie2=sortie1;      //sortie(n-2)
sortie1=sortie;       //sortie(n-1)


//sortie=entree+entree1*b1+entree2*b2+entree3*b3-sortie1*a1-sortie2*a2-sortie3*a3 ;   //filtre passe pas recursif ordre 3
sortie=(entree+entree1*b1+entree2*b2-sortie1*a1-sortie2*a2) ;                        //filtre passe pas recursif ordre 2
//sortie=(entree*0.19+sortie*0.81);                                 //filtre passe bas    a1=0.992 pour une constante de temps de 10ms  Te=1ms durre algo 150us                                                                                   //filtre passe pas recursif ordre 1
//sortie= entree  ;                                                                                     //gain unitaire
out=(sortie/Gain);
out=(out/4)  ;       //mise à l'echelle 10 bits en entrée et 8 bits en sortie                    
analogWrite(PWM3,out);          //    d

//digitalWrite(LED13,LOW);  

}//fin routine



///////////////////////////////////////////// Boucle correspondant à la fonction main
void loop() {  
  lcd.setCursor(0,0);
  lcd.print("fe=");
  lcd.print(fe,0);
    lcd.print(" fc=");
  lcd.print(fc,0);
  lcd.setCursor(0,1);
  lcd.print("b0=1 b1=2 b2=1 a0=1");
  lcd.setCursor(0,2);    
  lcd.print("a1=");
  lcd.print(a1,3);
   lcd.print(" a2=");
   lcd.print(a2,3);
   lcd.setCursor(0,3);  
   lcd.print("Gain=1/");
   lcd.print(Gain,2);



} // fin loop

Tous les essais se feront avec une fréquence d’échantillonnage de 1000Hz c’est à dire avec Te=1ms

Pour une fréquence de coupure de 10 Hz

Les coefficients obtenus sont :

b0=1 b1=2 b2=1 a0=1 a1 -1.911 et a2 =0.915 et le gain G=1/1000

Pour une fréquence de 10Hz l’atténuation est de 0.44

Pour une fréquence de 50Hz nous avons un atténuation de 0.08 = (0.2/(5x0.5))

Avec une fréquence de 100Hz la sortie est vraiment atténuée

Nous avons alors une atténuation correspondant à 0.02= (0.5x0.2)/5

Les valeurs de l’atténuation diffèrent en fonction de la fréquence choisie

Excellent travaille kabore pour Matlab, et pour avoir rentré le code
Pas facile, de mettre le code sous forme mathématique, c’est là que l’on voit les possibilités manquantes du compilateur IDE.
Il faudrait passer au 4eme ordre, fe=1000 Hz, Fcoupure=50Hz, juste pour savoir, si le micro peut toujours le faire ?

belle essai kaneD,
Normalemment, avec butterworth à la frequence de coupure l’atténuation est de 0.707
mais il est vrai 4carreau/7carreau=0.57 ????
l'entree est alimenté avec 7.1V donc l'entrée est saturé. donc voila l'erreur. D'ailleurs, le signal de sortie est legerement deformé et pas sinusoidale.
Attention, de mettre le zero des courbes toujours au meme endroit, car on croit qu’il y une composante continue à la sortie. etant donné, que la courbe sinosoidale d’entrée est unipolaire, le milieu est à 2.5V.

C’est bien beau, de faire que du filtre passe bas.
Mais faire comment on fait un filtre coupe bande pour laisser passer ?????????

En C++, il y a une library complexe.
Un soft peut faire tous des coefficients en ligne……celui-ci fonctionne tres bien
https://www-users.cs.york.ac.uk/~fisher/cgi-bin/mkfscript
on retrouve les valeurs du dexieme ordre


Il y a meme la courbe d’atténuation

Voici les coefficients pour Band-pass_filter entre 40Hz et 60Hz qu’il faudrait vérifier

il y a toujours des questions, c’est pénible.......mais on avance peu à peu

Bref, Pour mesurer l’atténuation et le déphasage du filtre numérique Arduino, il est possible de le faire directement avec le programme pour des signaux d’entrée alternatif
Pour le dephasage, on c’est inspirer des methodes de la mesure de puissance suivante
file:///C:/Users/geii/Downloads/Rapport_P6-3_2011_28.pdf

voici le code pour un passe bas avec frequence d’echantillonnage de 1kHz

#include <Wire.h>
//#include <LiquidCrystal_I2C.h>
#include <LiquidCrystal.h>
#include <SoftwareSerial.h>
#include <TimerOne.h>
//#include <math.h>


#define PWM3   3      //   timer2    
#define LED13    13       
LiquidCrystal lcd(9, 8, 4, 5, 6, 7);   // LiquidCrystal lcd(rs, en, d4, d5, d6, d7);
// Configuration des variables
//LiquidCrystal_I2C lcd(0x27, 20, 4);     //A0, A1, A2 non shunter 20x4ligne


unsigned    int    temps=0; 

            float  entree=0;
            float  entree1=0;
            float  entree2=0;
            float  entree3=0;
            float  entree4=0;

            float entreeMax;
            float entreeMin;
            float entreeMoy;
            float phase;
            float attenuation;
            float P;        //pour mesurer le dephasage
            float PMax;
            float PMin;
            float PMoy;
            
            float outMax;
            float Div;
            float factorPower;
//           float outMin;
//           float outMoy;
         
            float  sortie=0;
            float  sortie1=0;
            float  sortie2=0;
            float  sortie3=0;
            float  sortie4=0;
            
            float  out=0;

            const float b1 =4;            
            const float b2 =6;
            const float b3 =4;
            const float b4 =1;
            
            float a1=-3.181;
            float a2=3.86;
            float a3=-2.112;
            float a4=0.438;  
                       
            float denominateur=0;
            float Gain=3200;

            float  fe=1000;  //frequence d'echantilonnage   
            float  fc=50;  //frequence de coupure desiréé     
            

                    

void setup() {
  pinMode(LED13, OUTPUT);
  pinMode(PWM3,OUTPUT);

  Timer1.initialize(1000);           // initialize timer1, and set a 0,1 second period =>  100 000  pour 0.01s  10 000
  Timer1.attachInterrupt(callback);   // attaches callback() as a timer overflow interrupt
  lcd.begin(20, 4);                   //modifier pour un afficheur 20x4
/*
lcd.init();               //et pas lcd begin comme cetraine biblio
lcd.display();            // activer l'affichage
lcd.backlight();          // allumer retroeclairage  
*/  
  
  Serial.begin(9600); 

  TCCR2B = (TCCR2B & 0b11111000) | 0x01;         //pin 3  32khz    http://playground.arduino.cc/Main/TimerPWMCheatsheet
 
}



// Interruptions  tous les 1ms fait par le timer1    fe=1000Hz***********************************
void callback()  {
//if ( digitalRead(13)== 1 ) {digitalWrite(13,LOW);}  else {digitalWrite(13,HIGH);}
digitalWrite(LED13,HIGH);  //permet de mesurer à l'oscillo, le temps du calcul du filtre et le temps de la routine d'interruption
temps++;

entree=analogRead(A0); //convertisseur 10 bits sous 5V      https://www.gammon.com.au/adc
  
      
// ---- exemple filtre passe pas Butterwoth   pour fc=50hz,  fechantillon=1000Hz------    
entree4=entree3;      //entree(n-3)     
entree3=entree2;      //entree(n-3)
entree2=entree1;      //entree(n-2)
entree1=entree;       //entree(n-1)

sortie4=sortie3;      //sortie(n-3)
sortie3=sortie2;      //sortie(n-3)
sortie2=sortie1;      //sortie(n-2)
sortie1=sortie;       //sortie(n-1)

sortie=(entree+entree1*b1+entree2*b2+entree3*b3+entree4*b4-sortie1*a1-sortie2*a2-sortie3*a3-sortie4*a4)    ;
//sortie=(entree+entree1*b1+entree2*b2+entree3*b3-sortie1*a1-sortie2*a2-sortie3*a3) ;   //filtre passe pas recursif ordre 3
//sortie=(entree+entree1*b1+entree2*b2-sortie1*a1-sortie2*a2) ;                        //filtre passe pas recursif ordre 2
//sortie=(entree*0.19+sortie*0.81);                                 //filtre passe bas    a1=0.992 pour une constante de temps de 10ms  Te=1ms durre algo 150us                                                                                   //filtre passe pas recursif ordre 1
//sortie= entree  ;                                                                                     //gain unitaire
out=(sortie/Gain);
        
//rafraichissement toutes les 1s
if (temps>=1000)   {entreeMax=0;outMax=0;entreeMin=1023;PMin=1023;PMax=0;temps=0;}
if (entree> entreeMax) {entreeMax = entree; }  // Enregistrer la valeur maximale du capteur 
if (out> outMax) {outMax = out; }  // Enregistrer la valeur maximale du capteur 
//if (out< outMin) {outMin = out; }  // Enregistrer la valeur mini  
if (entree< entreeMin) {entreeMin = entree; }  // Enregistrer la valeur mini  
entreeMoy=(entreeMax+entreeMin)/2;
P=(entree-entreeMoy)*(out-entreeMoy);               //power instantanée  outmoy=entreemoy
if (P< PMin) {PMin = P; }
if (P> PMax) {PMax = P; }
//outMoy=(outMax+outMin)/2;

out=(out/4)  ;       //mise à l'echelle 10 bits en entrée et 8 bits en sortie 
if (out<0) {out=0;}
if (out>255) {out=255;}                                     
analogWrite(PWM3,out);
digitalWrite(LED13,LOW);  

}//fin routine



///////////////////////////////////////////// Boucle correspondant à la fonction main 
void loop() {  
/*  lcd.setCursor(0,0); 
  lcd.print("fe=");
  lcd.print(fe,0);
    lcd.print(" fc=");
  lcd.print(fc,0);
  lcd.setCursor(0,1); 
  lcd.print("b0=1 b1=3 b2=3 b3=1");
  lcd.setCursor(0,2);    
  lcd.print("a1=");
  lcd.print(a1,3);
   lcd.print(" a2=");
   lcd.print(a2,3);
   lcd.setCursor(0,3);  
  lcd.print("a3=");
     lcd.print(a3,3);  
   lcd.print(" G=1/");
   lcd.print(Gain,0);
*/
PMoy=((PMax+PMin)/2);   //valable seulement si 
Div=(outMax-entreeMoy)*(entreeMax-entreeMoy);
factorPower=(((PMoy*2)/Div));

lcd.setCursor(0,0); 
lcd.print("Emax=");     //amplitude max
lcd.print(entreeMax-entreeMoy,0);
lcd.print("   ");
lcd.setCursor(0,1); 
lcd.print("FactorPower");
lcd.print(factorPower,2);
lcd.print("   ");
  
lcd.setCursor(0,2);
lcd.print("attenuation=");
lcd.print((outMax-entreeMoy)/(entreeMax-entreeMoy));
lcd.print("  ");

phase=(acos(factorPower))*180/3.1416;   //dephasage par la methode des puissance

//phase=(acos(0.78)*180/3.1416);
  lcd.setCursor(0,3);
  lcd.print("phase=");
  lcd.print(phase,2);
  lcd.print((char)223);
  lcd.print("  ");


} // fin loop

Voici les simulations, on peut observer avec l’oscilloscope pour la fréquence de 4Hz, que l’atténuation est bien de 0.99 et le déphasage de 14° avec un écart de (0.5carreau/6carreau)*180°.
le signal vert est l’entrée et le jaune est la sortie.

On peut observer avec l’oscilloscope pour la fréquence de 25Hz, que l’atténuation est bien de 0.77 et le déphasage de 85° avec un écart de (1.9carreau/4carreau)*180°

On peut observer avec l’oscilloscope pour la fréquence de 50Hz, que l’atténuation est bien de 0.45 et le déphasage de 180° mais que la mesure donne 144° car la precision n’est pas terrible

pour la fréquence de 69Hz, que l’atténuation est bien de 0.19 et le déphasage de le déphasage de -52° avec un écart de (5carreau/7carreau)*180° , le signe moins n’est pas indiqué.

Pour la mesure déphasage, on aurait pu multiplier l’entrée et la sortie est utilisé un filtre passe bas qui est bien mieux que la mesure de la puissance max et mini pour en déduire la puissance moyenne.
De plus, la mesure de la puissance max et min est valable seulement si les signaux sont alternatifs. (pas de valeur moyenne)

Exemple de filtre numerique passe Bas sous excel pour mesurer la valeur moyenne d’une puissance monophasée

ce fichier peut etre telechargé ici

Pour améliorer la mesure de la phase donc de la mesure de la puissance moyenne précédente.
Voici, un nouveau programme qui met en exergue l’utilisation de filtre numérique
A1 est la tension qui doit être mesuré
A2 est la mesure de courant mesuré par un ACS712 (5A)
Pour faire la simulation, les 2 signaux sont identiques

La mesure de la fréquence d’entrée sera aussi mesurée.
Pour avoir la puissance active moyenne, il faut juste multiplier tension et courant puis filtrer l’harmonique 2 créé par cette multiplication. Il faut que cette multiplication (puissance instantanée) ait des valeurs négatives pour les déphasages supérieurs à 90°.
D’où l’utilisation de AOP pour avoir un offset de 2.5V.
La fréquence de coupure du filtre est de 5Hz, le programme peut donc fonctionner jusqu’à 500Hz thoriquement avec la frequence d’echantillonnafe de 1kHz

Pour une forme sinus, sans déphasage, avec les amplitudes suivantes la puissance active est bien de (1V*1A)/2=0.5WAtt

Pour une forme sinus, toujours sans déphasage, avec les amplitudes (0.5V*0.5A)/2=0.125Watt=0.125VA

Pour un retard de 2.5ms, donc déphasage de 45°, le facteur de puissance devrait être de 0.707 à la place de 0.73 de la mesure

Pour un retard de 5ms, donc un déphasage de 90°, le facteur de puissance est bien de 0.

Pour un retard de 6ms, donc un déphasage de 108°, le facteur de puissance est bien de -0.3

Pour un retard de 10ms, donc un déphasage de 180°, le facteur de puissance est bien de -1

Pour une fréquence de 250Hz, donc un déphasage de 90°, le facteur de puissance est de 0.17, il y a une erreur un peu plus importante.

Pour une fréquence de 450Hz, cela marche très bien aussi et même jusqu’à 1kHz correspondant à la fréquence d’échantillonnage grâce à la duplication spectrale

Par conséquent, la précision de la mesure du déphasage du facteur de puissance est bien meilleure qu’avec la méthode précédente

Voici le programme

#include <Wire.h>
//#include <LiquidCrystal_I2C.h>
#include <LiquidCrystal.h>
#include <SoftwareSerial.h>
#include <TimerOne.h>
//#include <math.h>


#define PWM3   3      //   timer2    
#define LED13    13       
LiquidCrystal lcd(9, 8, 4, 5, 6, 7);   // LiquidCrystal lcd(rs, en, d4, d5, d6, d7);
// Configuration des variables
//LiquidCrystal_I2C lcd(0x27, 20, 4);     //A0, A1, A2 non shunter 20x4ligne


unsigned    int    temps=0; 

            float  entree=0;
            float  entree1=0;
            float  entree2=0;

            float Tension;
            float Courant;
            float TensionMax;
            float CourantMax;
            float GAinCourant=2;
            float factorPower;
         
            float  sortie=0;
            float  sortie1=0;
            float  sortie2=0;
           
            float  powermoy=0;   //average power
            float ApparePower;
            
           
            float a1=-1.956;     //fc5hz
            float a2=0.957;
                       
            float denominateur=1;
            float Gain=1;

            float  fe=1000;  //frequence d'echantilonnage   
            float  fc=5;  //frequence de coupure desiréé     
            
      unsigned int TO;
      unsigned int Time=1;
        

void setup() {
  pinMode(LED13, OUTPUT);
  pinMode(PWM3,OUTPUT);

  Timer1.initialize(1000);           // initialize timer1, and set a 0,1 second period =>  100 000  pour 0.01s  10 000
  Timer1.attachInterrupt(callback);   // attaches callback() as a timer overflow interrupt
  lcd.begin(20, 4);                   //modifier pour un afficheur 20x4
/*
lcd.init();               //et pas lcd begin comme cetraine biblio
lcd.display();            // activer l'affichage
lcd.backlight();          // allumer retroeclairage  
*/  
  
 Serial.begin(9600); 

TCCR2B = (TCCR2B & 0b11111000) | 0x01;         //pin 3  32khz    http://playground.arduino.cc/Main/TimerPWMCheatsheet

attachInterrupt(0, interrup2, FALLING);   // broche2
TO = micros();


    denominateur=(fe*fe+fe*fc*PI*sqrt(2)+fc*fc*PI*PI);
    a1=(fc*fc*PI*PI-fe*fe)*2/denominateur;
    a2=(fe*fe*fe*fe+(fc*fc*fc*fc*PI*PI*PI*PI))/(denominateur*denominateur)         ;                             //pow(base, exponent)
    Gain=(1+2+1)/(1+a1+a2);

}



void interrup2() // la fonction appelée par l'interruption externe n°0
{
Time=(micros()-TO);  //mesure du temps
TO = micros();
}





// Interruptions  tous les 1ms fait par le timer1    fe=1000Hz***********************************
void callback()  {
//if ( digitalRead(13)== 1 ) {digitalWrite(13,LOW);}  else {digitalWrite(13,HIGH);}
//digitalWrite(LED13,HIGH);  //permet de mesurer à l'oscillo, le temps du calcul du filtre et le temps de la routine d'interruption
temps++;

Tension=analogRead(A0);         //convertisseur 10 bits sous 5V      https://www.gammon.com.au/adc
Courant=analogRead(A1);         //courant     
//rafraichissement mesure toutes les 1s
entree=(((Courant-512)/GAinCourant)*(Tension-500));               //power instantanée  outmoy=entreemoy

//rafraichissement toutes les 2s
if (temps>=2000)   {CourantMax=0;TensionMax=0;temps=0;  }
if (Tension> TensionMax) {TensionMax = Tension; }  // Enregistrer la valeur maximale  amplitude
if (Courant> CourantMax) {CourantMax = Courant; }  // Enregistrer la valeur maximale


      
// -- filtre passe pas Butterwoth  2ordre ,  fechantillon=1000Hz------    
entree2=entree1;      //entree(n-2)
entree1=entree;       //entree(n-1)

sortie2=sortie1;      //sortie(n-2)
sortie1=sortie;       //sortie(n-1)

sortie=(entree+entree1*2+entree2-sortie1*a1-sortie2*a2) ;        //filtre passe pas recursif ordre 2  correspondant la puissance moyenne     
powermoy=(sortie/Gain);                                 // powermoy= ((CourantMax*TensionMax)/1024*2)*facteurpower                                              

//visuliation de la puissance intantannée inutile
//powermoy=(entree/(4))  ;      
//powermoy=(powermoy/4);          //mise à l'echelle 10 bits en entrée et 8 bits en sortie
//if (powermoy<0) {powermoy=0;}
//if (powermoy>255) {powermoy=255;}                                     
//analogWrite(PWM3,powermoy);

//digitalWrite(LED13,LOW);  

}//fin routine



///////////////////////////////////////////// Boucle correspondant à la fonction main 
void loop() {  
 
lcd.setCursor(0,0); 
lcd.print("F(Hz)");     //
fe=1000000/Time;        //redressement double alternance time en microseconde
lcd.print(fe);
lcd.print("   ");

ApparePower=(((CourantMax-512)/GAinCourant)*(TensionMax-500))/(2);
factorPower=((powermoy)/(ApparePower));
lcd.setCursor(0,1); 
lcd.print("FactorPower");
lcd.print(factorPower,2);
lcd.print("   ");
  
lcd.setCursor(0,2);
lcd.print("ActivPower=");
lcd.print((powermoy/41943),2);   //il y a deja une division de 1024 pour le filtre et avoir 7 chiffres significatif   powermoy
lcd.print("Watt");

lcd.setCursor(0,3);
lcd.print("ApparPower=");
lcd.print(ApparePower/41943,2);   //5*5/1024*1024   ApparePower
lcd.print("VA");


} // fin loop

On aurait pu faire le calcul de la fréquence de coupure en en fonction de la mesure de la fréquence du signal. Sachant que le temps des instructions de la routine d’interruption est de 0.4ms.

Sinon, il y a PZEM-004T mais qui fonctionne seulement de 45 à 65Hz

sinon, le calcul de la valeur moyenne glissante avec des données dans un tableau pour mesurer en TRMS serait à tester…

Pour bien comprendre une mesure moyenne par tableau, voici un petit programme pour mesurer la valeur moyenne ou efficace d’un signal.
La longueur de la table va dépendre de la fréquence de signal.
Mais pourquoi ne pas faire une table de 200 valeurs « int », sur 0.2 s avec la période d’échantillon de 1ms.
car l’ATmega328P a 2koctet de SRAM.

faut-il utiliser une moyenne glissante ?
Un exemple de la moyenne glissante, ici, mais il aurait du utiliser un filtre passe bas, cela prend moins de mémoire.

A la place de faire une moyenne glissante autant faire par paquet de mesure de 0.2s et faire un affichage avec cette période.
De plus, additionner 200 valeurs pour faire une somme, cela prend 2ms correspondant à l’incrémentation dans la table et l’addition 6µs.

Pour une fréquence d’échantillonnage de 1000Hz, alors que dans la mesure analogique et la valeur dans la table dure seulement 0.1ms donc on pourrait augmenter cette valeur.
Avec un signal carré de 50Hz et un rapport cyclique de 0.75, il y aura 20 échantillons utiles de mesure correspondant à 1000Hz/50Hz tout le reste est périodique et redondant.
Donc une précision 5% sur la mesure moyenne.
Avec une fréquence de 250Hz, la précision passe à 25%, ce n’est pas pas tolérable

En changeant, le prescaleur à 16 pour analogread, on passe à 20µs en lecture.
Donc, la fréquence d’échantillonnage peu passer à 10Khz pour avoir une meilleure précision.
Mais cela augmente énormément, le nombre de valeur dans la table.
Par conséquent, la table sera fixe et la fréquence d’échantillonnage variera sur la fréquence du signal a mesuré pour avoir 0.5% sur la mesure mais avec une limitation de la fréquence d’échantillonnage à 10Khz

Pour 50hz, précision à 0.5%, la mesure 1023*0.75 donne bien 768.

Pour 500hz, précisons de 5%, la valeur moyenne mesurée se dégrade légèrement mais sans osciller

Pour 1KHz, précision de 10%, la valeur moyenne oscille entre 692 et 844

Pour 10Hz, la précision est bonne aussi, théoriquement de 0.5%

Voici le code

#include <Wire.h>
//#include <LiquidCrystal_I2C.h>
#include <LiquidCrystal.h>
#include <SoftwareSerial.h>
#include <TimerOne.h>
//#include <math.h>


#define PWM3   3      //   timer2    
#define LED13    13       
LiquidCrystal lcd(9, 8, 4, 5, 6, 7);   // LiquidCrystal lcd(rs, en, d4, d5, d6, d7);
// Configuration des variables
//LiquidCrystal_I2C lcd(0x27, 20, 4);     //A0, A1, A2 non shunter 20x4ligne


unsigned    int    temps=0; 


            float Tension;
            float Courant;
            float GAinCourant=2;

            float somme;
            float Vaverage;

  
            float  fmesure=1;  //frequence de la mesure
            
            float TO;
            float Time=1;

      unsigned int table[200];
               byte tablemax=200; 
                       
            unsigned int Te=100;       

void setup() {
  pinMode(LED13, OUTPUT);
  pinMode(PWM3,OUTPUT);

  Timer1.initialize(100);           // initialize timer1, and set a 0,1 second period =>  100 000  pour 0.01s  10 000
  Timer1.attachInterrupt(callback);   // https://www.pjrc.com/teensy/td_libs_TimerOne.html
  lcd.begin(20, 4);                   //modifier pour un afficheur 20x4
  
 Serial.begin(9600); 

TCCR2B = (TCCR2B & 0b11111000) | 0x01;         //pin 3  32khz    http://playground.arduino.cc/Main/TimerPWMCheatsheet

attachInterrupt(0, interrup2, FALLING);   // broche2
TO = micros();

//ADCSRA |= bit (ADPS0) | bit (ADPS1) | bit (ADPS2);   // 128
ADCSRA &= ~(bit (ADPS0) | bit (ADPS1) | bit (ADPS2)); // clear prescaler bits
  ADCSRA |= bit (ADPS2);                           //  16 
}



void interrup2() // la fonction appelée par l'interruption externe n°0
{
Time=(micros()-TO);  //mesure du temps
TO = micros();
}





// Interruption fait par le timer1    fe=10KHz***********************************
void callback()  {
//if ( digitalRead(13)== 1 ) {digitalWrite(13,LOW);}  else {digitalWrite(13,HIGH);}
//digitalWrite(LED13,HIGH);  //permet de mesurer à l'oscillo, le temps du calcul du filtre et le temps de la routine d'interruption
temps++;
Tension=analogRead(A0);         //convertisseur 10 bits sous 5V      https://www.gammon.com.au/adc
table[temps]=Tension;
//digitalWrite(13,LOW);
}//fin routine



///////////////////////////////////////////// Boucle correspondant à la fonction main 
void loop() {  
 
lcd.setCursor(0,0); 
lcd.print("F(Hz)");     //
fmesure=1000000/Time;        //redressement double alternance time en microseconde
lcd.print(fmesure,0);
lcd.print("   ");

 //changement de la frequence d'echnatillonnnage en fonction de la frequence de mesure 
if (fmesure>=50)    {Timer1.initialize(100);}     // Fechantillon=10KHz
if (fmesure<50)    {Te=-20*fmesure+1100;Te=Te/2;Timer1.initialize(Te);}     //

if (temps>=tablemax)  {
Timer1.stop();
for(int i=0;i<tablemax; i++)  {somme+=table[i];
Serial.print(table[i]); Serial.print(";");   
}
Serial.println(";");
 Vaverage=somme/tablemax;
 temps=0;
Timer1.restart();                                               
}

lcd.setCursor(0,1); 
lcd.print("Vmoyen");     //
lcd.print(Vaverage*5/1024,2);
lcd.print("V   ");


} // fin loop

Donc, il est possible d’avoir la mesure de la puissance moyenne avec une table de 200 valeurs, avec une assez bonne précision jusqu’à 1kHz.
La table de 200 valeurs utilise 40% de la mémoire SRAM atmega328,

Avec cette méthode, il est possible de mesurer la valeur moyenne de n’importe quel type de signal
Exemple pour un signal mono alternance, puisque le micro ne peut lire les valeurs négatives. La tension moyenne correspond bien de Amplitude/PI=4V/3.14=1.25V

Pour un signal triangulaire, la tension moyenne correspondant à l’amplitude/4=4V/4=1V.

Donc, on peut faire la valeur efficace vraie TRMS avec cette methode, mais on sort du sujet

Une autre façon de mesurer le facteur de puissance ou le déphasage.
De façon tres simple avec 2 AOPs……

Filtre butterworth passe bas, d’ordre 4

Utilisation du programme exécutable matlab pour retrouver les valeurs des différentes coefficients a et b du filtre, comme pour l’ordre 3 expliqué à la première page.

Cette méthode parce que en développant les calculs nous-mêmes et l’intégrant sur le compilateur arduino, cela devient de plus en plus complexe. Sauf si dans l’avenir les developpeurs arduino arrivent a pallier ce problème de la variable float.

Voici le programme qui permet de calculer des coefficients sur matlab pour un butterworth 4 ordre à partir de la fréquence d’échantillonnage (Fe=1000) et de la fréquence de coupure (Fc=50hz)

N= 4           % ordre  du filtre
fc= 50;        % Hz freq de coupure
fe= 1000;       % Hz sample freq d'echantillonage te=1ms
[b,a]= butter_synth(N,fc,fe)

Etant donné que c’est un ordre 4, les coefficients de b sont pris par rapport au coef du triangle de pascal
D’où b0=1, b1=4, b2=6, b3=4, b4=1

Il ne reste plus qu’à compiler sur arduino

voici le code

#include <Wire.h>
#include <LiquidCrystal_I2C.h>
#include <SoftwareSerial.h>
#include <TimerOne.h>
//#include <math.h>


#define PWM3   3      //   timer2    
#define LED13    13      
// Configuration des variables
LiquidCrystal_I2C lcd(0x27, 20, 4);     //A0, A1, A2 non shunter 20x4ligne


unsigned    int    temps=0;

            float  entree=0;
            float  entree1=0;
            float  entree2=0;
            float  entree3=0;
            float  entree4=0;

        
            float  sortie=0;
            float  sortie1=0;
            float  sortie2=0;
            float  sortie3=0;
            float  sortie4=0;
            
            float  out=0;

            const float b1 =4;            
            const float b2 =6;
            const float b3 =4;
            const float b4 =1;
            
            float a1=-3.1806;
            float a2=3.8612;
            float a3=-2.1122;
            float a4=0.4383;  
                      
            float denominateur=0;
            float Gain=0;

            float  fe=1000;  //frequence d'echantilonnage  
            float  fc=50;  //frequence de coupure desiréé    
            

                    

void setup() {
  pinMode(LED13, OUTPUT);
  pinMode(PWM3,OUTPUT);

  Timer1.initialize(1000);           // initialize timer1, and set a 0,1 second period =>  100 000  pour 0.01s  10 000
  Timer1.attachInterrupt(callback);   // attaches callback() as a timer overflow interrupt
  lcd.begin(20, 4);                   //modifier pour un afficheur 20x4
/*
lcd.init();               //et pas lcd begin comme cetraine biblio
lcd.display();            // activer l'affichage
lcd.backlight();          // allumer retroeclairage  
*/  
  
  Serial.begin(9600);

  TCCR2B = (TCCR2B & 0b11111000) | 0x01;         //pin 3  32khz    http://playground.arduino.cc/Main/TimerPWMCheatsheet

 Gain=(1+4+6+4+1)/(1+a1+a2+a3+a4);
}



// Interruptions  tous les 1ms fait par le timer1    fe=1000Hz***********************************
void callback()  {
//if ( digitalRead(13)== 1 ) {digitalWrite(13,LOW);}  else {digitalWrite(13,HIGH);}
digitalWrite(LED13,HIGH);  //permet de mesurer à l'oscillo, le temps du calcul du filtre et le temps de la routine d'interruption
temps++;

entree=analogRead(A0); //convertisseur 10 bits sous 5V      https://www.gammon.com.au/adc
  
      
// ---- exemple filtre passe pas Butterwoth   pour fc=50hz,  fechantillon=1000Hz------    
entree4=entree3;      //entree(n-3)    
entree3=entree2;      //entree(n-3)
entree2=entree1;      //entree(n-2)
entree1=entree;       //entree(n-1)

sortie4=sortie3;      //sortie(n-3)
sortie3=sortie2;      //sortie(n-3)
sortie2=sortie1;      //sortie(n-2)
sortie1=sortie;       //sortie(n-1)

sortie=(entree+entree1*b1+entree2*b2+entree3*b3+entree4*b4-sortie1*a1-sortie2*a2-sortie3*a3-sortie4*a4)    ;
//sortie=(entree+entree1*b1+entree2*b2+entree3*b3-sortie1*a1-sortie2*a2-sortie3*a3) ;   //filtre passe pas recursif ordre 3
//sortie=(entree+entree1*b1+entree2*b2-sortie1*a1-sortie2*a2) ;                        //filtre passe pas recursif ordre 2
//sortie=(entree*0.19+sortie*0.81);                                 //filtre passe bas    a1=0.992 pour une constante de temps de 10ms  Te=1ms durre algo 150us                                                                                   //filtre passe pas recursif ordre 1
//sortie= entree  ;                                                                                     //gain unitaire
out=(sortie/Gain);
      

out=(out/4)  ;       //mise à l'echelle 10 bits en entrée et 8 bits en sortie
if (out<0) {out=0;}
if (out>255) {out=255;}                                    
analogWrite(PWM3,out);
digitalWrite(LED13,LOW);  

}//fin routine



///////////////////////////////////////////// Boucle correspondant à la fonction main
void loop() {  
   lcd.setCursor(0,0);
   lcd.print("G=1/");
   lcd.print(Gain,0);

  //lcd.print(fe,0);
    lcd.print(" fc=");
  lcd.print(fc,0);
  
  lcd.setCursor(0,1);
  lcd.print("b0=1 b1=3 b2=3 b3=1");
  lcd.setCursor(0,2);    
  lcd.print("a1=");
  lcd.print(a1,3);
   lcd.print(" a2=");
   lcd.print(a2,3);
   lcd.setCursor(0,3);  
  lcd.print("a3=");
  lcd.print(a3,3);  
  lcd.print(" a4=");
  lcd.print(a4,3);  

} // fin loop

Pour une fréquence de coupure de Fc=10Hz et une période d’échantillonnage de 1ms
On a les résultats de simulations suivantes:
4,75V en sortie sur 5V en entrée … une atténuation de 0.9

Pour une fréquence de coupure de Fc=50Hz
3.5V/5V=0,7 en atténuation

Pour une fréquence de coupure de Fc=70Hz
1.4V/5V=0,28 en atténuation

L’ordre 4 tient ses promesses, il filtre mieux que les précédents. Et on constate que l’arduino arrive aussi a gérer le calcul d’ordre 4.

Nous avons eu à faire pas mal d’essai avec Butterworth mais seulement avec le filtre passe bas sous différents ordres ( de 1 à 4), nous avons alors trouvé intéressant de faire des essais le filtre Butterwrth passe bande.

On faisait le calcul de coefficients avec Mathcad et qui donnait le résultat dépendant de la fréquence de coupure et d’échantillonnage.

Il existe En C++, il y a une librarie complexe.
Un soft qui peut faire tous des coefficients en ligne et donc nous n’auront qu’a rentrer le type de filtre désiré, les fréquences de coupures ,la fréquence d’échantillonnage et l’ordre du filtre.

Dans notre cas nous avons choisi un filtre type Butterworth passe bande second ordre avec une fréquence de coupure basse de 40Hz une fréquence de coupure haute de 60Hz et une Fréquence d’échantillonnage de 1000.

Ainsi il nous fournit les coefficients du filtre ainsi que le gain

https://www-users.cs.york.ac.uk/~fisher/cgi-bin/mkfscript

Voici le code sur Arduino

#include <Wire.h>
#include <LiquidCrystal_I2C.h>
//#include <LiquidCrystal.h>
#include <SoftwareSerial.h>
#include <TimerOne.h>
#include <math.h>


#define PWM3   3      //   timer2   
#define LED13    13     
LiquidCrystal_I2C lcd(0x27, 16, 2);
//LiquidCrystal_I2C lcd(9, 8, 4, 5, 6, 7);   // LiquidCrystal lcd(rs, en, d4, d5, d6, d7);
// Configuration des variables
//LiquidCrystal_I2C lcd(0x27, 20, 4);     //A0, A1, A2 non shunter 20x4ligne


unsigned    int    temps=0;

            float  entree=0;
            float  entree1=0;
            float  entree2=0;
            float  entree3=0;
            float  entree4=0;
       
            float  sortie;
            float  sortie1=0;
            float  sortie2=0;
            float  sortie3=0;
            float  sortie4=0;
           
            float  out=0;
            const float b1 =0;
            const float b2 =-2;
            const float b3=0;
            const float b4=1;

            float a1=3.642;
            float a2=-5.1461;
            float a3=3.3324;
            float a4=-0.8371;
            float denominateur=0;
            float Gain=0;

            float  fe=1000;  //frequence d'echantilonnage 
            float  fc=50;  //frequence de coupure desiréé   
           

                   

void setup() {
  pinMode(LED13, OUTPUT);
  pinMode(PWM3,OUTPUT);

  Timer1.initialize(1000);           // initialize timer1, and set a 0,1 second period =>  100 000  pour 0.01s  10 000
  Timer1.attachInterrupt(callback);   // attaches callback() as a timer overflow interrupt
  lcd.begin(20, 4);                   //modifier pour un afficheur 20x4


  Serial.begin(9600);

  TCCR2B = (TCCR2B & 0b11111000) | 0x01;         //pin 3  32khz    http://playground.arduino.cc/Main/TimerPWMCheatsheet

Gain=(2.761027558e+02);  //gain donné par le soft pour le calcul des coefficients
 
}


// Interruptions  tous les 1ms fait par le timer1    fe=1000Hz***********************************
void callback()  {
if ( digitalRead(13)== 1 ) {digitalWrite(13,LOW);}  else {digitalWrite(13,HIGH);}
digitalWrite(LED13,HIGH);  //permet de mesurer à l'oscillo, le temps du calcul du filtre et le temps de la routine d'interruption
entree=analogRead(A0); //convertisseur 10 bits sous 5V
 
     
// ---- filtre passe bande Butterwoth , fechantillon=1000Hz------       
entree4=entree3;      //entree(n-4)
entree3=entree2;       //entree(n-3)
entree2=entree1;      //entree(n-2)
entree1=entree;       //entree(n-1)

sortie4=sortie3;      //sortie(n-4)
sortie3=sortie2;       //sortie(n-3)
sortie2=sortie1;      //sortie(n-2)
sortie1=sortie;       //sortie(n-1)


  

sortie=(b4*entree4)+(b2*entree2)+entree +(a4*sortie4)+(a3*sortie3)+(a2*sortie2)+(a1*sortie1);  //filtre passe bande second ordre
//gain unitaire

out=(sortie/Gain);
out=(out/4)  ;       //mise à l'echelle 10 bits en entrée et 8 bits en sortie         
if (out<0) {out=0;}
if (out>255) {out=255;}           
analogWrite(PWM3,out);          //    

digitalWrite(LED13,LOW); 

}//fin routine



///////////////////////////////////////////// Boucle correspondant à la fonction main
void loop() { 
   lcd.setCursor(0,0);
  
   lcd.print("b0=1 b1=0");
 
  lcd.setCursor(0,1);
  lcd.print(" b2=-2 b3=0 b4=1");
  lcd.setCursor(0,2);   
  lcd.print("a1=");
  lcd.print(a1,3);
   lcd.print(" a2=");
   lcd.print(a2,3);
   lcd.setCursor(0,3); 
  lcd.print("a3=");
  lcd.print(a3,3); 
  lcd.print(" a4=");
  lcd.print(a4,3); 

   



} // fin loop

Nous allons faire des simulations sur ISIS pour vérifier les atténuations de la sortie par rapport à l’entrée en choisissant différentes valeurs de fréquences d’entrée

Pour une fréquence d’entré de 10Hz nous avons une très forte atténuation

sortie/ entrée = (0.1/2) * 0.5 = 0.05
L’atténuation est de 0.05 donc très grande

Vérifions la valeur de l’atténuation pour une fréquence d’entrée de 50Hz

Pour une fréquence de 50Hz l’atténuation est quasi-nulle
4.85*5 = 0.96 donc sensiblement égale à 1

Rentrons à présent les fréquences de coupures hautes et basses et vérifions l’atténuation de la sortie

Pour F=40Hz


Nous avons une atténuation de 0.36 soit 1.8/5

Pour F=60Hz


L’atténuation ici correspond également à 0.36

Donc avec le filtre passe bande nous avons la même atténuation aux deux fréquences de coupures choisies

Une atténuation nulle à la fréquence centrale de 50Hz et une atténuation trop forte pour 10Hz

Cependant nous devrions obtenir pour ces fréquences de coupures une atténuation de 0.707 correspondant au gain à -3dB mais on obtient cette valeur qu’avec les fréquences de 46 et 54 Hz soient plus 6Hz de la fréquence basse et -6Hz de la fréquence haute ce qui est un peu bizarre

Pour 46Hz nous avons ce résultat

3.5/5 = 0.7

Pour F=54 Hz

Nous obtenons une attenuation de 3.7/5 = 0.74

Le filtre passe bande ne fonctionne pas alors car nous n’avons pas les atténuations attendues

Sur les figures precedentes, on peut observer qu’il n’y a pas la partie negative sur la sortie.
C’est normal car c’est un passe bande, donc la sortie du filtre n’aura pas de valeur de moyenne autour de 2.5V ne sera pas present
donc il faut rajourer à la PWM la constante 127.(voir le code)

Mais, effectivement, le passe bande avec une frequence de 1000Hz demande des coefficients avec 4 à 5 chiffres apres la virgules sinon cela diverge.
Voici le code, attention, les coefficients sont opposés par rapport à votre programme

#include <Wire.h>
//#include <LiquidCrystal_I2C.h>
#include <LiquidCrystal.h>
#include <SoftwareSerial.h>
#include <TimerOne.h>
//#include <math.h>


#define PWM3   3      //   timer2    
#define LED13    13       
LiquidCrystal lcd(9, 8, 4, 5, 6, 7);   // LiquidCrystal lcd(rs, en, d4, d5, d6, d7);
// Configuration des variables
//LiquidCrystal_I2C lcd(0x27, 20, 4);     //A0, A1, A2 non shunter 20x4ligne


unsigned    int    temps=0; 

            float  entree=0;
            float  entree1=0;
            float  entree2=0;
            float  entree3=0;
            float  entree4=0;
         
            float  sortie=0;
            float  sortie1=0;
            float  sortie2=0;
            float  sortie3=0;
            float  sortie4=0;
            
            float  out=0;
float entreeMax;
float outMax;
float attenuation;

            
 //        Fc 40hz et 60 Hz  fe=1000Hz  ok  
             float a1=-3.6427;
            float a2=5.14618;
            float a3=-3.33247;
            float a4=0.83718;
            float Gain=300;     //
            
 //        Fc 10hz et 90 Hz  fe=1000Hz  divergence         
 /*           float a1=-3.24;   // -3.24;
            float a2=4;
            float a3=-2.25;    //-2.25
            float a4=0.5;
            float Gain=21;           
*/               
     
 //        Fc 10hz et 90 Hz  fe=500Hz  ok mais attenue aussi haute vitesse
 /*            float a1=-2.46;       //-2.46
            float a2=2.31;
            float a3=-1.08;
            float a4=0.25;
            float Gain=8;
 */            
 //      Fc 25Hz et 75 Hz  fe=500Hz  divergence
 /*           float a1=-3.43;       //-3.43
            float a2=4.54;
            float a3=-2.7;
            float a4=0.64;
            float Gain=50;
 */         
 //      Fc 10Hz et 50 Hz  fe=500Hz ok mais attenue aussi haute vitesse
  /*           float a1=-3.17;       
            float a2=3.88;
            float a3=-2.2;
            float a4=0.5;
            float Gain=30;   */
                    

void setup() {
  pinMode(LED13, OUTPUT);
  pinMode(PWM3,OUTPUT);

  Timer1.initialize(1000);           //   pour 0.001s  1000   0.002s 2000
  Timer1.attachInterrupt(callback);   // attaches callback() as a timer overflow interrupt
  lcd.begin(20, 4);                   //modifier pour un afficheur 20x4
/*
lcd.init();               //et pas lcd begin comme cetraine biblio
lcd.display();            // activer l'affichage
lcd.backlight();          // allumer retroeclairage  
*/  
  
  Serial.begin(9600); 

  TCCR2B = (TCCR2B & 0b11111000) | 0x01;         //pin 3  32khz    http://playground.arduino.cc/Main/TimerPWMCheatsheet

}



// Interruptions  tous les 1ms fait par le timer1    fe=1000Hz***********************************
void callback()  {
//if ( digitalRead(13)== 1 ) {digitalWrite(13,LOW);}  else {digitalWrite(13,HIGH);}
//digitalWrite(LED13,HIGH);  //permet de mesurer à l'oscillo, le temps du calcul du filtre et le temps de la routine d'interruption
entree=analogRead(A0); //convertisseur 10 bits sous 5V
  
      
// ---- exemple filtre passe pas Butterwoth   ,  fechantillon=1000Hz------    
entree4=entree3;      //entree(n-4)
entree3=entree2;      //entree(n-3)
entree2=entree1;      //entree(n-2)
entree1=entree;       //entree(n-1)

sortie4=sortie3;
sortie3=sortie2;      //sortie(n-3)
sortie2=sortie1;      //sortie(n-2)
sortie1=sortie;       //sortie(n-1)

sortie=(entree4-2*entree2+entree-a4*sortie4-a3*sortie3-a2*sortie2-a1*sortie1);  //filtre passe bande second ordre
//sortie=(entree+entree1*b1+entree2*b2+entree3*b3+entree4*b4-sortie1*a1-sortie2*a2-sortie3*a3-sortie4*a4)    ;
//sortie=(entree+entree1*b1+entree2*b2+entree3*b3-sortie1*a1-sortie2*a2-sortie3*a3) ;   //filtre passe pas recursif ordre 3
//sortie=(entree+entree1*b1+entree2*b2-sortie1*a1-sortie2*a2) ;                        //filtre passe pas recursif ordre 2
//sortie=(entree*0.19+sortie*0.81);                                 //filtre passe bas    a1=0.992 pour une constante de temps de 10ms  Te=1ms durre algo 150us                                                                                   //filtre passe pas recursif ordre 1
//out= entree  ;                                                                                     //gain unitaire
out=(sortie/Gain);
out=((out/4)+127);       //mise à l'echelle 10 bits en entrée et 8 bits en sortie   
if (out<0) {out=0;}
if (out>255) {out=255;}                  
analogWrite(PWM3,out);          //    
//Serial.println(out);
temps++;
//rafraichissement toutes les 0.2s
if (temps>=200)   {entreeMax=512;outMax=127;temps=0;}
if (entree> entreeMax) {entreeMax = entree; }  // Enregistrer la valeur maximale du capteur
if (out> outMax) {outMax = out; } 
//digitalWrite(LED13,LOW);  
}//fin routine



///////////////////////////////////////////// Boucle correspondant à la fonction main 
void loop() {  
  lcd.setCursor(0,0);    
  lcd.print("passe bande 2eme");
  lcd.setCursor(0,1);
  attenuation=((outMax-127)*4/(entreeMax-512));    
  lcd.print("A=");
  lcd.print(attenuation,2);
   lcd.print("  ");



} // fin loop

un petit rappel du filtre passe bande theorique du deuxieme ordre analogique, avec son coefficient de qualité surnomée Q

On peut observer sur la figure suivante avec une fe=1000hz et les frequences de coupures de 40Hz et 60Hz
Pour la fréquence d’entrée de 25Hz, atténuation correspond bien à la théorie
Atténuation=(Frequence entree)/(frequence coupureQ)=(25/(502.5))=0.2

A 40Hz, on a bien une atténuation de 0.707
A 50Hz, il y avait une petite saturation, de la sortie donc à la place de prendre le gain théorique de 276, nous avons choisi un gain de 300

A 60Hz, on a bien une attenuation de 0.707

A 80Hz, on a bien une attenuation tres importante

Conclusions, nous avions bien un filtre passe bande autour de 50Hz avec une bande passante de 20Hz
Si l’on ne prend que 2 chiffres après les virgules, il y a une divergence de la sortie.

Le temps de calcul de ce deuxième ordre dure 0.3ms donc inferieur à notre routine d’interruption de 1ms