ECG AD8232 nano programme du monitoring, BPM heart rate

International Français Réalisations et Projets Finis
1

il existe des shields qui permettant de faire l’Electrocardiogramme portatif pas cher

1. la bibliographie
Il y a de nombreux blogs qui explique ce circuit intégré, mais sans aller bien loin dans le programme
Ce blog permet d’avoir juste un moniteur grâce à la liaison série

De même ici, qui explique un peu mieux le capteur
https://learn.sparkfun.com/tutorials/ad8232-heart-rate-monitor-hookup-guide?_ga=2.54916523.1050236860.1613662700-154441868.1613662700

un article « Prototype Low-cost Portable Electrocardiogramme (ECG) Based on Arduino-Uno with Bluetooth Feature »
https://aip.scitation.org/doi/pdf/10.1063/1.5139392
mais sans partage du programme, un autre article
« Towards IoT and ML Driven Cardiac Status Prediction System »
https://www.researchgate.net/publication/336988184_Towards_IoT_and_ML_Driven_Cardiac_Status_Prediction_System
mesure BPM (battement per minute heart rate)

un ECG IOT avec un ESP32

Malgré une grosse recherche pas trouver de programme open source complet intéressant pour ce capteur ????

2. Des questions de démarrages
Quelle est la sensibilité fait par le capteur AD832 ? Donc quelle est le signal récupéré ?
Quelle fréquence d’échantillon doit-on utiliser ?
Est-ce que ce signal est fort perturbé par la CEM ? Peut-on mesurer les temps QST ?
Peut-on utiliser la fonction pulseIn (Signal_GBF, LOW); pour avoir le BPM ?
Sur le forum, il y aurait quelques soucis d’utilisation de ce capteur « Problem with AD8232 ECG WaveForm »
https://forum.arduino.cc/?topic=706502#msg4749936
https://forum.arduino.cc/index.php?topic=405315.0

Il ne faut pas trop bouger car toute activité musculaire provoque des bio potentiels supérieurs à celle du cœur et vont saturer le capteur.

Le signal du battement du cœur est entre 0.66Hz (40BPM periode=151 centi seconde=1510milliseconde) à 3.33Hz (200BPM periode=30 centi second) donc si le signal a une fréquence bien plus haut, il y a une erreur. D’ailleurs l’ampli a une bande passante de 0.4Hz à 30Hz
Le gain est de 60dB donc une Amplification de 1000


En enregistrant les données dans un fichier CSV pouvant être lu par Excel, on peut s’apercevoir que les points RST sont très rapides et avec une fréquence d’échantillonnage de 0.01s un seul échantillon est mesuré
En faisant la dérivée d’Euler on peut détecter très facilement RST.
Pour lire le battement cardiaque, la valeur absolue d’Euler donne encore une meilleure précision comme on peut l’observer sur la figure suivante


Suite à differents essais, le signal peut être très variable en fonction du placement des électrodes et du type de personne.

Mais grâce à la dérivée d’Euler, il est facile de détecter le BPM

Mais, si la valeur moyenne n’est pas aux alentours de 1.65V (337 dec), le BPM ne sera pas correcte
La valeur moyenne peut être mesurée sur 3 secondes
Avec l’équation suivante Vmoy= SommeV/temps
Il est possible de mesurer la sensibilité du signal tous les 3 secondes pour validée le BPM ou pas

3. Programme de base V1

Donc voici notre premier programme qui donne le monitoring sur le traceur série ou dans un fichier CSV, le BPM en PJ
Le signal avec un peu de perturbation de 50Hz sur le traceur serie

4. Conclusion :
Grâce au math et au filtre numérique, il est possible de détecter le BPM assez facilement.
Mais si l’on désire mesurer les temps PR et ST, il faut des temps d’échantillons plus rapide donc utilisé une ESP32 OLED.
mais pour avoir des données pertinentes, il va falloir utiliser des filtres.

5. Perspectives
Souvent, il y a une mesure des impulsions cardiaques sur plusieurs heures pour savoir s’il y a des anomalies.
Par conséquent, le monitoring doit enregistrées les données.
Avec une liaison Bluetooth HC06, un smartphone pourrait faire le monitoring et la sauvegarde des données.
Un filtre numérique permettrait d’atténuer les perturbations électromagnétiques du 50Hz du signal bio potentiel.
Le BPM évolue sur plusieurs minutes, donc un filtre numérique passe bas devrait être utilisé.

Il est possible de simuler le programme sous ISIS en simulant le capteur avec un datafile ce qui permet de mieux vérifier le fonctionnement du programme.

ECG.ino (3.36 KB)

2

voila la moniteur cardio sur smarthphone
mais il faut configurer d’abord le HC06 avec la commande AT+UART=57600,0,0


En effet le Bluetooth à 9600 baud prend plus de temps que la période d’échantillonnage 0.01s que l’on a choisi imposé par la routine d’interruption.

De plus à la place d’utiliser les broches RX 0 et TX 1 pour la communication avec le Bluetooth, l’utilisation des broches 10 et 11 est plus pertinente avec la ligne SoftwareSerial hc06(10, 11);
Ce qui permet de ne pas avoir de conflit de communication lorsqu’une nouvelle version de programme est transférée…

voici le programme

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

SoftwareSerial hc06(10, 11);   //nano  
#define LED13    13       
LiquidCrystal lcd(9, 8, 4, 5, 6, 7);   // LiquidCrystal lcd(rs, en, d4, d5, d6, d7);

#define PI 3.141592
unsigned    int    temps=0;
unsigned    int    entree=0;
unsigned    int   entreen1=0;
unsigned    int  entreeMax=0;
unsigned    int  entreemini=0;
 int derivation;
 float derivationmoy;
 int detection;
 int detection1;
 int deriveImpulse;

unsigned int  centi;
unsigned int  periode;
byte tflag;
bool flag=0;
 
unsigned   int  sensibilite=0; 
unsigned   int  BPM=0;

bool LOm;
bool LOp;
float Inaverage=337;
float Inaverage1;

           
void setup() {
  pinMode(LED13, OUTPUT);
  pinMode(2, INPUT); // Setup for leads off detection LO +
  pinMode(3, INPUT); // Setup for leads off detection LO -

  
  Timer1.initialize(10000);           //   pour 0.001s  1000   0.002s 2000     0.01s 10000
  Timer1.attachInterrupt(callback);   // attaches callback() as a timer overflow interrupt
  lcd.begin(20, 4);                   //modifier pour un afficheur 20x4

  Serial.begin(57600);   //     broche RX, TX 0 et 1  moniteur serie
  hc06.begin(57600);    //AT+UART=57600,0,0   broche 10 et 11


lcd.setCursor(0,0);    
lcd.print("ECG  IUT soissons");
lcd.setCursor(0,1);    
lcd.print("BPM");

}



// Interruptions  tous les 0.01s fait par le timer1    fe=100Hz***********************************
void callback()  {
centi++;
temps++;

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

  entreen1=entree;
  entree=analogRead(A0);         //convertisseur 10 bits sous 5V 

  derivation=abs(entree-entreen1);

  derivationmoy=derivationmoy+derivation;  //la derivation moyenne est souvent du au signal 50Hz
  derivationmoy=derivationmoy/temps;        //
  detection=derivation-derivationmoy;       //retire la perturbation du 50Hz

 tflag++;
if (detection>170) {periode=centi;  centi=0;flag=1;tflag=0;
    if (periode<200 && periode>30)   {BPM=6000/periode;  }}    //60*100/periode
if (flag==1 && tflag>=5) {flag=0;}       //delay de 0.05

//LOm=digitalRead(2);
//LOp=digitalRead(3);
 
//traceur serie    
//Serial.print(entree); Serial.print(",");Serial.println(detection);   //Serial.print(",");Serial.println(BPM);  
//fichier bluetooth electronic    https://www.keuwl.com/apps/bluetoothelectronics/         http://www.keuwl.com/apps/bluetoothelectronics/userguide/remote_code.html
  if (hc06.available()) {
  Serial.write(hc06.read());  
 }


hc06.print("*G");hc06.print(String(entree)+","+String(detection));hc06.println("*");   //15ms à 9600 bauds  mais 300us à 115200 bauds
hc06.print("*b");hc06.print(String(BPM));hc06.println("*");
hc06.print("*p");hc06.print(String(periode));hc06.println("*");


//rafraichissement tous les 5s          
if (temps>=500)   {temps=0; }

 

}//fin routine




void loop() { 

} // fin loop

Avec 115200 bauds, la routine d’interruption dure 3ms
Avec 57600 bauds, la routine d'interruption dure 5ms, donc il reste du temps pour du filtre numérique du troisième ordre.

Maintenant, il va falloir essayer différents filtres pour améliorer la mesure de l’ECG

3
  1. Traitement numérique de l’ECG

Différents filtres numériques et leurs études, ont été postés sur ce lien
Filtre numerique RII, RIF….digital filter...FFT...atmega328, ESP32
mais il n’y a pas d’explication sur les filtres coupes bandes numeriques ????????

Mais quel est le spectre du signal d’un ECG ?
Quelle fréquence de coupure utilisée ? Quel ordre de filtre faut il utilser ?
Est ce l’ATMEGA 328 va pouvoir faire le job ou faudra-t-il prendre un processeur plus puissant ?

Voici quelques bibliographies sur le signal de l’ECG et de leur traitement :
Un article intéressant donnant les valeurs de filtres du second ordre numérique bande passante
http://dobigeon.perso.enseeiht.fr/teaching/signal/MODAP_TS_URSAFE.pdf

Une thèse « Traitement des signaux à échantillonnage irrégulier Application au suivi temporel de paramètres cardiaques »

Cette thèse montre la difficulté de faire le traitement du signal d’un ECG.

Pour mettre au point et vérifier le traitement numérique et son bon fonctionnement à partir d’un signal qui est enregistré dans un fichier CSV. Il est possible d’utiliser les softs suivants

  • utilisation de MATLAB
  • utilisation d’un tableur tel que EXCEL
  • utilisation d’un datafile dans ISIS qui permet de simuler arduino nano
  • utiliser directement le processeur en temps reels mais ce sera bein moins pertinent que les autres possiblilités car le signal d’entrée sera différent

comme tout le monde ne peuvent utiliser MATLAB et ISIS qui ne sont pas freeware
le tableau Excel ou google sheet sera utilisé.

La perturbation électrique la plus importante est bien sur le 50Hz, donc il faut réaliser un filtre coupe bande ou rejecteur ou à encoche

https://www.sonelec-musique.com/electronique_bases_filtre_audio.html
il y a une belle explication ici
http://w3.cran.univ-lorraine.fr/perso/hugues.garnier/Enseignement/TdS/J-TdS_Conception_filtres_num.pdf

avec une frequence de coupure de 2*fe, les coeficients z-1 s’en vont et demande peu de calcul.
On peut observer que le coeficient a permet d’avoir plus ou moins d’attenuation, mais provoque un gain statique diffèrent de 1.
Il faut au moins 4 échantillons pour que le filtre fonctionne, donc il faut au moins une fréquence d’échantillonnage de 200Hz avec une routine d’interruption de 5ms

Sur la courbe suivante avec fe=200Hz, on peut observer la formule du filtre dans la case C4.
La courbe bleu est le signal non traité, la courbe orange est celle filtré ou il n’y a plus le signal 50Hz, la courbe grise correspond à la dérivée d’Euler pour détecter les BPM grâce à l’onde RS
On peut remarquer que l’on distingue aussi l’onde T facilement

Il n’y a pas beaucoup de changement sur le signal de ECG avec un coefficient a de 0.7ou 0.9


Un filtre passe bas de 10hz avec une fréquence d’échantillonnage de 100Hz est idiot car il va atténuer la dynamique rapide de l’onde RS comme on peut l’observer sur la figure suivante et la detection du BPM sera plus difficile à faire.


Voici le nouveau programme avec le filtre réjecteur coupe bande et la frequence d’echantillon de 200Hz

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

SoftwareSerial hc06(10, 11);   //nano  
#define LED13    13       
LiquidCrystal lcd(9, 8, 4, 5, 6, 7);   // LiquidCrystal lcd(rs, en, d4, d5, d6, d7);

#define PI 3.141592
unsigned    int    temps=0;
unsigned    int    entree=0;
unsigned    int   entreen1=0;
unsigned    int   entreen2=0;

float filtre;
float filtren1;
float filtren2;


 int derivation;

unsigned int  centi;
unsigned int  periode;
byte tflag;
bool flag=0;
 
unsigned   int  sensibilite=0; 
unsigned   int  BPM=0;

bool LOm;
bool LOp;
float Inaverage=337;
float Inaverage1;

           
void setup() {
  pinMode(LED13, OUTPUT);
  pinMode(2, INPUT); // Setup for leads off detection LO +
  pinMode(3, INPUT); // Setup for leads off detection LO -

  
  Timer1.initialize(5000);           //   pour 0.001s  1000   0.002s 2000     0.01s 10000
  Timer1.attachInterrupt(callback);   // attaches callback() as a timer overflow interrupt
  lcd.begin(20, 4);                   //modifier pour un afficheur 20x4

  Serial.begin(57600);   //     broche RX, TX 0 et 1  moniteur serie
  hc06.begin(57600);    //AT+UART=57600,0,0   broche 10 et 11    routine dure 5ms
//  hc06.begin(115200);   //routine dure 3ms


lcd.setCursor(0,0);    
lcd.print("ECG  IUT soissons");
lcd.setCursor(0,1);    
lcd.print("BPM");
}



// Interruptions  tous les 0.005s fait par le timer1    fe=200Hz***********************************
void callback()  {
centi++;
temps++;

//if ( digitalRead(13)== 1 ) {digitalWrite(13,LOW);}  else {digitalWrite(13,HIGH);}
digitalWrite(13,HIGH);
  entreen2=entreen1;
  entreen1=entree;
  entree=analogRead(A0);         //convertisseur 10 bits sous 5V 

  filtren2=filtren1;
  filtren1=filtre;
  filtre=entree+entreen2-0.49*filtren2;  //0.7^2=0.49 

  derivation=4*abs(filtre-filtren1);   //6 coeficient d'amplification

 tflag++;
if (derivation>300) {periode=centi;  centi=0;flag=1;tflag=0;
    if (periode<400 && periode>60)   {BPM=12000/periode;  }}    //60*200/periode      <2s && >0.3s
if (flag==1 && tflag>=5) {flag=0;}       //delay de 0.05


//traceur serie    
Serial.print(entree); Serial.print(",");Serial.print(filtre,0);Serial.print(",");Serial.print(derivation);Serial.print(",");Serial.println(BPM);  
//fichier bluetooth electronic     http://www.keuwl.com/apps/bluetoothelectronics/userguide/remote_code.html
/*
  if (hc06.available()) {
  Serial.write(hc06.read());  }


//hc06.print("*G");hc06.print(String(entree)+","+String(detection));hc06.println("*");   //15ms à 9600 bauds  mais 300us à 115200 bauds   
//hc06.print("*b");hc06.print(String(BPM));hc06.println("*");                            //avec 57600 bauds, la routine d'interruption dure 5ms
//hc06.print("*p");hc06.print(String(periode));hc06.println("*");
*/

//rafraichissement tous les 5s          
if (temps>=500)   {temps=0; }
digitalWrite(13,LOW);

}//fin routine

void loop() { 

} // fin loop

Voici le signal reçu en temps réel, par le traceur série avec filtre et sans filtre……lorsqu’on ne bouge pas
On peut observer que l’atténuation du 50Hz est bien effective et que la détection du BPM est stable, ainsi que la détection de l’onde T.

Par contre, si l’on bouge et qu'il y a les biopotentiels des autres muscles (respiration, bras), la détection est plus aléatoire de l’onde T qui a la même fréquence que le mouvement n'est plus identifiable.

A cause des biopotentiels du mouvement, il faudrait aussi filtrer l’onde de détection par un filtre passe bas de fréquence de coupure de 10Hz.

il y a une etude de filtre RII coupe bande sur Github avec un autre AOP

4

Filtrage numérique du signal de détection
Sur les figures précédentes, on peut observer que le signal de détection du BMP correspondant à la variation de l’onde rapide RS est bruité à cause de sa valeur absolue.
De plus il y a l’onde T qui pourrait être mieux détecté.

Pour cela, un filtre passe bas du seconde ordre avec une fréquence de coupure de 15Hz a été choisi.
Le programme permet de changer la fréquence et les coefficients sont déterminés par le programme.
Ce poste démontre comment il est possible de tester le filtre.

Mar 25, 2020 cela a été démontré sur ce poste Re: Filtre numerique RII, RIF….

Etant donné que nous avons enregistré, le signal ECG brute dans un fichier ISIS, il est possible de simuler le programme dans proteus electronic ISIS pour vérifier son bon fonctionnement.
Il y a un affichage des coefficients du filtre passe bas en fonction de la fréquence de coupure sur la figure suivante
Par contre en simulation, il n’y a pas de traceur série
mais il aurait été possible du simuler un SSD oled

Sinon faire une sortie analogique d’une variable pour l’observer à l’oscilloscope

Sinon faire en réel, ce que fait le filtre sur la détection (valeur absolue de la dérivée), on peut observer sur la figure suivante le signal de détection (vert) et celui qui est filtré passe bas numérique (orange)
Donc, il est donc possible de mesurer les temps PR, ST, QRS….par contre l’onde U n’est pas facile à observer.
Une autre solution est de retirer la valeur moyenne sur l’ECG car en fonction de la respiration ou des mouvements, cette valeur moyenne varie.


Électrocardiogramme est bien expliqué dans Wikipédia et en fonction du placement des électrodes
[https://fr.wikipedia.org/wiki/Électrocardiographie#:~:text=L'électrocardiographie%20(ECG)%20est,et%20la%20conduction%20des%20influx](http://\\\\\"http://\\"https://fr.wikipedia.org/wiki/Électrocardiographie#:~:text=L'électrocardiographie (ECG) est,et la conduction des influx"").

Par conséquent, même avec un atmega328, il est possible de mesurer les temps PR, ST, QRS et de diagnostiquer des effets anormaux du cœur.

Tout peut être envisagé en programme et traitement du signal, il faut juste faire les bons choix entre les avantages et inconvénients des temps de calculs de fréquence d’échantillonnage, du nombre d’échantillonnage

Exemple sur la figure suivante l’ECG en bleu, sa dérivée de l’onde en rouge et filtré en vert avec le même filtre passe bas précèdent.


Des flags permettent de connaitre les temps des dynamiques de l’ECG à partir de chaque changement de pente


le nouveau programme de détection du BPM à partir de la valeur négative de la dérivée précédente et de la mesure du tableau de variation de l’ECG avec des flags en piece jointe

Maintenant, il faudrait détecter par le programme, les anomalies du signal ECG….
Mais, ce n’est pas si facile de faire un programme d’ECG lorsqu’on n’est pas spécialiste du cœur ?

Remarques :
Il y a plusieurs possibilités de dérivées numériques du signal qui peut être effectué pour déterminer les différentes phases de l’onde du signal.

Quel est le pas optimal de la dérivée donc est ce que la fréquence d’échantillonnage est bien choisi ?
est ce que notre choix de dérivée d’Euler est bien choisi ? quelle est l’erreur de mesure ?
https://fr.wikipedia.org/wiki/D%C3%A9riv%C3%A9e
thèse : Dérivation numérique : synthèse, application et intégration
file:///C:/Users/geii/Downloads/TH_T2240_mdridi.pdf
Détermination du pas optimal dans le calcul des dérivées
https://core.ac.uk/download/pdf/194545998.pdf

Aurait-on pu aussi déterminer le BPM et les autres caractéristiques avec la Fast Fourier Transfert ?

ECG_V3.ino (4.37 KB)

5

FFT de l’ECG
pas facile de comprendre les articles sur l’ECG quand on n’est pas spécialiste

https://regi.tankonyvtar.hu/hu/tartalom/tamop412A/2011_0079_jobbagy_biomedical/ch06s06.html

Comme on n’est jamais mieux servi que par soi-même,
Etant donné qu’Excel fait la fast fourrier transfert et que le signal est dans ce logiciel autant,
http://blog.aannaque.org/2020/04/analyse-de-fourrier-sur-excel.html#

Voici la FFT du signal ECG sans le filtre coupe bande de 50Hz, avec 256 échantillons*0.005=1,28s.
Il faut un nombre d’échantillon correspondant à la période du signal qui est de 1Hz (60BPM)
Le pas du domaine fréquentielle la fft correspond à la fréquence de coupure divisé par le nombre d’échantillon 200H/256=0.78Hz
On peut observer que l’amplitude du signal 100 de 50Hz est prépondérante par rapport au signal de l’ECG


Avec le filtre coupe bande la fréquence, la raie des 50Hz disparait. On peut observer que l’amplitude du signal est faible par rapport à la valeur moyenne de 456.

La fréquence de 1Hz du BPM n’est pas trop retrouvée à cause du « phénomène de fuite » car la fréquence d’échantillon n’est pas un multiple du nombre d’échantillon.
C’est pour cela qu’il faut que la fréquence d’échantillon soit un multiple du BPM ou avoir un nombre d’échantillon très grand.


En prenant 2048 données, avec un pas fréquentielle des raies tous les 0.1Hz ; on peut observer que l’amplitude du signal du BPM, correspondant à l’onde RS qui la raie la plus importante.
Donc la précision de la BPM par la mesure de la fréquence correspondra à l’équation suivante :
Précision BPM=0.1Hz*60=±6 BPM ce qui n’est pas terrible.


Avec 2048 valeurs à la fréquence d’échantillon de 200Hz, cela enregistre 5.12s de données.
Le processeur a besoin de mémoire mais au niveau du temps de calcul, ce sera possible avec un ATMega 328
La dérivée du signal en comptant, le temps entre chaque variation est donc bien plus rapide et precis que la méthode de la FFT

À mon avis, il serait intéressant de mesurer la variation du BPM sur 10 minutes, pour savoir s'il y a une perte du signal ou pas.
Etant donné qu'il faut 4 à 5 minutes pour qu'il y a une variation du BPM de 60 à 90 lors d'un effort, il faudrait certainement filtrer la valeur du BPM avec un fréquence de coupure 30s

6

Avec une fréquence d’échantillonnage=1s et une fréquence de coupure de 0.016Hz, d’un filtre passe bas sur le BPM
Voici les coefficients du filtre passe bas et la simulation sous ISIS

Voici l’étude fréquentielle théorique de ce filtre numerique qui confirme bien les valeurs des coefficients

Voici le nouveau programme

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

SoftwareSerial hc06(10, 11);   //nano  
#define LED13    13       
LiquidCrystal lcd(9, 8, 4, 5, 6, 7);   // LiquidCrystal lcd(rs, en, d4, d5, d6, d7);

#define PI 3.141592
unsigned    int    temps=0;
unsigned    int    entree=0;
unsigned    int   entreen1=0;
unsigned    int   entreen2=0;

float filtre;
float filtren1;
float filtren2;

float derivation;
float derivation1;
float derivation2;
float derivationABS;

unsigned int  centi;
unsigned int centiST;
float ST;      //seconde
unsigned int  periode; 

bool flag=0;
bool flag1=0;
bool flag2=0;

 
unsigned   int  sensibilite=0; 
float BPM=0;  //battement per minute
float BPMmo=60;
float BPMmoy;

float BPMn1=60;
float BPMn2=60;
float BPMFiltre;
float BPMFiltre1=60;
float BPMFiltre2=60;
float BPMFiltreC=60;

bool LOm;
bool LOp;
float Inaverage=337;
float Inaverage1;


            float a1=1;
            float a2=1;

            float detectionfiltre;
            float detectionfiltren1;
            float pentedetection;
            
            float sortie;
            float sortie1;
            float sortie2;
                       
            float denominateur=0;
            float Gain=1;

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

           
void setup() {
  pinMode(LED13, OUTPUT);
  pinMode(2, INPUT); // Setup for leads off detection LO +
  pinMode(3, INPUT); // Setup for leads off detection LO -

  
  Timer1.initialize(5000);           //   pour 0.001s  1000   0.002s 2000     0.01s 10000
  Timer1.attachInterrupt(callback);   // attaches callback() as a timer overflow interrupt
  lcd.begin(20, 4);                   //modifier pour un afficheur 20x4

  Serial.begin(57600);   //     broche RX, TX 0 et 1  moniteur serie
//  hc06.begin(57600);    //AT+UART=57600,0,0   broche 10 et 11    routine dure 5ms

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);
/*
  lcd.setCursor(0,0); lcd.print("fe=");  lcd.print(fe,1);  lcd.print(" fc=");  lcd.print(fc,3);
  lcd.setCursor(0,1); lcd.print("b0=1 b1=2 b2=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("G=1/"); lcd.print(Gain,0);
  */

}

// Interruptions  tous les 0.005s fait par le timer1    fe=200Hz***********************************
void callback()  {
centi++;
temps++;
if (flag==1) {centiST++;}

//if ( digitalRead(13)== 1 ) {digitalWrite(13,LOW);}  else {digitalWrite(13,HIGH);}
//digitalWrite(13,HIGH);
  entreen2=entreen1;
  entreen1=entree;
  entree=analogRead(A0);         //convertisseur 10 bits sous 5V 

  filtren2=filtren1;
  filtren1=filtre;                       //0.7^2=0.49
  filtre=entree+entreen2-0.49*filtren2;  //filtre rejecteur

  derivation=4*(filtre-filtren1);   //4 coeficient d'amplification
if (derivation>=900) {derivation=900;} 
if (derivation<-500) {derivation=-500;} 
 
//  derivationABS=abs(derivation);
//  if (derivationABS>=900) {derivationABS=900;} 
  
  sortie2=sortie1;
  sortie1=sortie;
  sortie=(derivation+derivation1*2+derivation2+sortie1*1.36-sortie2*0.52) ;    //filtre passe pas recurssif ordre 2  fe=200  fc=15
detectionfiltren1=detectionfiltre;
detectionfiltre=sortie/3;     //  detectionfiltre=sortie/Gain;
if (detectionfiltre>=900) {detectionfiltre=900;} 
if (detectionfiltre<-500) {detectionfiltre=-500;}   

 
if (detectionfiltre<-300) {periode=centi;  centi=0;flag=1;centiST=0;        //mesure le BPM en detectant l'onde QRS negatif
    if (periode<400 && periode>60)   {BPM=12000/periode;  }}    //60*200/periode      <2s && >0.3s
if (detectionfiltre>100 && flag==1) {flag1=1;  }
if (detectionfiltre<-30 && flag1==1) {if (centiST<40) {ST=centiST/200;}  centiST=0; flag=0;flag1=0;}  // mesure temps ST
 

//traceur serie    
//Serial.print(filtre,0); Serial.print(",");Serial.print(detectionfiltre,0);Serial.print(",");Serial.print(pentedetection,0);//Serial.print(",");Serial.print(flag*200);Serial.print(",");Serial.println(BPM);  
//Serial.print(filtre,0); Serial.print(",");Serial.print(detectionfiltre,0);Serial.print(",");Serial.println(flag*400);

//rafraichissement tous les 1s          
if (temps>=200)   {
  digitalWrite(13,HIGH);
  BPMn2=BPMn1;
  BPMn1=BPM;
  BPMFiltre2=BPMFiltre1;
  BPMFiltre1=BPMFiltreC;
  BPMFiltreC=(BPM+BPMn1*2+BPMn2-BPMFiltre1*a1-BPMFiltre2*a2); 
 // BPMFiltre=BPMFiltreC/Gain;
  BPMFiltre=BPMFiltreC/415;
  Serial.print(BPM);Serial.print(",");Serial.println(BPMFiltre,1);temps=0; flag=0;flag1=0;BPMmo=0;digitalWrite(13,LOW);    }


//digitalWrite(13,LOW);

}//fin routine



void loop() { 

} // fin loop

Le résultat fourni par le traceur série avec la mesure du BPM qui est très stable malgré que l’on bouge légèrement devant l’ordi.

Quelle serait la différence de résultat de la mesure du BPM entre le filtre passe bas et le calcul d’une valeur moyenne de BPM sur 60s ?

7

Avant de répondre à la question précédente :
Il faut faire une fréquence d’échantillonnage de plus de 1 Hz car la fréquence du BPM peut être d’inférieur à 0.5Hz avec un BPM de 30 par exemple.
Donc, une période d’échantillonnage de 6s (0.16Hz pour le filtre passe bas) a été choisie arbitrairement

Pour répondre à la question précédente, une valeur moyenne du BMP qui sera réinitialisé tous 60s, va être utilisé.
Ce filtre moyenneur est très facile à programmer par rapport au filtre passe pas du seconde ordre car c’est une intégration numérique sur un certain temps.
Voici les nouveaux coefficients du filtre passe bas du BPM


Compter les BPM pendant 60s est divisé par le nombre de mesure, cela correspond à un filtre RII (réponse impulsionnel infinie) avec tous les coefficients à 1.
Dans ce cas, la fréquence de coupure correspond à l’équation suivante :
Fc=(Fechantillon0.5)/nbr d’echantillon=(1Hz0.5)/61=0.008Hz (120s)

D’ailleurs, voici l’étude fréquentielle du filtre moyenner qui a bien une fréquence de coupure de 0.008Hz (120s) avec alors que celle du filtre passe bas précédent est de 60s.


Voici les résultats en réel. On peut apercevoir que la détection du BPM peut faire des erreurs car le testeur bouge mais qui sont passablement filtré avec une amplitude de ±15BPM pour le filtre passe bas et ±10BPM pour le filtre moyenneur.


Donc pour le filtre passe bas numérique a avec une fréquence de coupure de 1/120s a été testé pour par rapport au filtre moyenneur pour faire une comparaison
On peut remarquer que le filtre moyenneur en temps réell au début de la réinitialisation peut avoir des informations erronées par rapport à celle qui donne la valeur au moins 61 echantillons.


Conclusion le filtre moyenneur sur une minute est plus efficace que le filtre passe pas pour minimiser les erreurs de mesures
Mais il est possible aussi de melanger le filtre moyenneur et le filtre passe bas


voici le nouveau programme

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

SoftwareSerial hc06(10, 11);   //nano  
#define LED13    13       
LiquidCrystal lcd(9, 8, 4, 5, 6, 7);   // LiquidCrystal lcd(rs, en, d4, d5, d6, d7);

#define PI 3.141592
unsigned    int    temps=0;
unsigned    int    Minute=0;
unsigned    int    entree=0;
unsigned    int   entreen1=0;
unsigned    int   entreen2=0;

float filtre;
float filtren1;
float filtren2;

float derivation;
float derivation1;
float derivation2;
float derivationABS;

unsigned int  centi;
unsigned int centiST;
float ST;      //seconde
unsigned int  periode; 

bool flag=0;
bool flag1=0;
bool flag2=0;

 
unsigned   int  sensibilite=0; 
float BPM=60;  //battement per minute
float n=0;
float BPMmo;
float BPMmoy;
float BPMmoy1;

float BPMn1=60;
float BPMn2=60;
float BPMFiltre;
float BPMFiltre1=60;
float BPMFiltre2=60;
float BPMFiltreC=60;

bool LOm;
bool LOp;
float Inaverage=337;
float Inaverage1;


            float a1=1;
            float a2=1;

            float detectionfiltre;
            float detectionfiltren1;
            float pentedetection;
            
            float sortie;
            float sortie1;
            float sortie2;
                       
            float denominateur=0;
            float Gain=1;

            float  fe=0.166;  //frequence d'echantilonnage  periode 6S 
            float  fc=0.0083;  //frequence de coupure desiréé  60s

           
void setup() {
  pinMode(LED13, OUTPUT);
  pinMode(2, INPUT); // Setup for leads off detection LO +
  pinMode(3, INPUT); // Setup for leads off detection LO -

  
  Timer1.initialize(5000);           //   pour 0.001s  1000   0.002s 2000     0.01s 10000
  Timer1.attachInterrupt(callback);   // attaches callback() as a timer overflow interrupt
  lcd.begin(20, 4);                   //modifier pour un afficheur 20x4

  Serial.begin(57600);   //     broche RX, TX 0 et 1  moniteur serie
//  hc06.begin(57600);    //AT+UART=57600,0,0   broche 10 et 11    routine dure 5ms

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);
/*
  lcd.setCursor(0,0); lcd.print("fe=");  lcd.print(fe,2);  lcd.print(" fc=");  lcd.print(fc,4);
  lcd.setCursor(0,1); lcd.print("b0=1 b1=2 b2=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("G=1/"); lcd.print(Gain,0);
*/

}

// Interruptions  tous les 0.005s fait par le timer1    fe=200Hz***********************************
void callback()  {
centi++;   //detection 
temps++;   //detection chaque echantillon
Minute++;  //detection Minute
if (flag==1) {centiST++;}

//if ( digitalRead(13)== 1 ) {digitalWrite(13,LOW);}  else {digitalWrite(13,HIGH);}
//digitalWrite(13,HIGH);
  entreen2=entreen1;
  entreen1=entree;
  entree=analogRead(A0);         //convertisseur 10 bits sous 5V 

  filtren2=filtren1;
  filtren1=filtre;                       //0.7^2=0.49
  filtre=entree+entreen2-0.49*filtren2;  //filtre rejecteur

  derivation=4*(filtre-filtren1);   //4 coeficient d'amplification
if (derivation>=900) {derivation=900;} 
if (derivation<-500) {derivation=-500;} 
 
//  derivationABS=abs(derivation);
//  if (derivationABS>=900) {derivationABS=900;} 
  
  sortie2=sortie1;
  sortie1=sortie;
  sortie=(derivation+derivation1*2+derivation2+sortie1*1.36-sortie2*0.52) ;    //filtre passe pas recurssif ordre 2  fe=200  fc=15
detectionfiltren1=detectionfiltre;
detectionfiltre=sortie/3;     //  detectionfiltre=sortie/Gain;
if (detectionfiltre>=900) {detectionfiltre=900;} 
if (detectionfiltre<-500) {detectionfiltre=-500;}   

 
if (detectionfiltre<-300) {periode=centi;  centi=0;flag=1;centiST=0;        //mesure le BPM en detectant l'onde QRS negatif
    if (periode<400 && periode>60)   {BPM=12000/periode;n++;BPMmo=BPM+BPMmo;BPMmoy=BPMmo/n;}}    //60*200/periode      <2s && >0.3s
if (detectionfiltre>100 && flag==1) {flag1=1;  }
if (detectionfiltre<-30 && flag1==1) {if (centiST<40) {ST=centiST/200;}  centiST=0; flag=0;flag1=0;}  // mesure temps ST
 

//traceur serie    
//Serial.print(filtre,0); Serial.print(",");Serial.print(detectionfiltre,0);Serial.print(",");Serial.print(pentedetection,0);//Serial.print(",");Serial.print(flag*200);Serial.print(",");Serial.println(BPM);  
//Serial.print(filtre,0); Serial.print(",");Serial.print(detectionfiltre,0);Serial.print(",");Serial.print(flag1*100,0);Serial.println(",");

//rafraichissement tous les 6s          
if (temps>=1200)   {
  digitalWrite(13,HIGH);
  BPMn2=BPMn1;
  BPMn1=BPMmoy1;
  BPMFiltre2=BPMFiltre1;
  BPMFiltre1=BPMFiltreC;
  BPMFiltreC=(BPMmoy1+BPMn1*2+BPMn2-BPMFiltre1*a1-BPMFiltre2*a2); 
  BPMFiltre=BPMFiltreC/Gain;

//  Serial.print(BPM,0);Serial.print(",");Serial.print(BPMmoy,0);Serial.print(",");Serial.print(BPMFiltre,0);Serial.print(",");Serial.print(BPMmoy1,0);Serial.println(",");
  Serial.print(BPM,0);Serial.print(",");Serial.print(BPMFiltre,0);Serial.println(",");
temps=0; flag=0;flag1=0;digitalWrite(13,LOW);    
                    }

if (Minute>=12000)   {BPMmoy1=BPMmoy;BPMmo=60;n=0;Minute=0;}                    

//digitalWrite(13,LOW);

}//fin routine



void loop() { 
//lcd.setCursor(8,3); lcd.print("BPM"); lcd.print(BPMFiltre);

} // fin loop

il aurait aussi été possible de faire un tableau de 60 valeurs du BPM pour en faire une moyenne glissante
Moyenne mobile — Wikipédia.

8

programmation pour la mesure du BPM avec

valeur moyenne glissante (Auto-Regressive Mobile Average)

Le calcul de la valeur moyenne glissante utilise N de mémoire en utilisant l’équation récurrente suivante par contre n’utilise pas beaucoup de temps de calcul en utilisant la première équation par rapport à la troisième équation

La programmation peut donc être complément différent en fonction de l’équation mathématique utilisé.

Voici quelques liens comment coder une valeur moyenne glissante

Mais quelles sont les caractéristiques de ce filtre moyenneur mobile ?
Quelle sera les limites de la fréquence d’entrée en fonction du nombre d’échantillon ?
Quelle sera la fréquence de coupure ?

Peu d’informations dans la littérature
http://paristech.institutoptique.fr/site.php?id=753&fileid=9818
[http://www.isir.upmc.fr/UserFiles/File/zarader/Cours_Traitement%20_Signal_P2.pdf](http://www.isir.upmc.fr/UserFiles/File/zarader/Cours_Traitement _Signal_P2.pdf)
ARMA — Wikipédia

Voici les caractéristiques du filtre en fonction du nombre d’échantillon qui est ici de 10 echantillons


Pour que la valeur moyenne fonctionne, la relation entre la fréquence minimum, la fréquence d’échantillonnage et le nombre d’échantillon doit correspondre à l’inéquation suivante :

nbr d’echantillon*periode d’echantillonnage>>periode du signal traité

quelque soit le nombre d’echantillon, l’atténuation correspond à un filtre du seconde ordre avec une fréquence de coupure de
Fc=(Fechantillon*0.5)/nbr d’echantillon

Pour simplifier voici le programme d’un filtre moyenneur seule à partir de l’équation précédente pour un tableau de 10 entier, donc qui va utiliser un nombre d’octet 102+index1+average2+total4=28 octet

Mais le compilateur dit qu’il utilise 41 octets de variable sur les 2k, c’est quand même dérisoire

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

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

uint8_t  temps=0; 
uint16_t  average;
const byte numReadings  = 10;  //https://www.aranacorp.com/fr/implementation-de-la-moyenne-glissante-dans-arduino/
uint16_t  readings [numReadings];  //https://www.arduino.cc/en/Tutorial/BuiltInExamples/Smoothing
uint8_t  readIndex  = 0;
uint32_t  total  = 0;
                  
           
void setup() {
  pinMode(LED13, OUTPUT);
  
  Timer1.initialize(5000);           //   pour 0.001s  1000   0.002s 2000     0.01s 10000
  Timer1.attachInterrupt(callback);   // attaches callback() as a timer overflow interrupt
//  lcd.begin(20, 4);                   //modifier pour un afficheur 20x4

Serial.begin(57600);   //     broche RX, TX 0 et 1  moniteur serie
//lcd.setCursor(0,0); lcd.print("IUT GEII Soissons");  
}

// Interruptions  tous les 0.005s fait par le timer1    fe=200Hz***********************************
void callback()  {

//if ( digitalRead(13)== 1 ) {digitalWrite(13,LOW);}  else {digitalWrite(13,HIGH);}
//digitalWrite(13,HIGH);
temps++;
if (temps>=20)  {                             //tous les 0.1s
      total = total - readings[readIndex];         //filtre moyenneur glissant, lissage du BPM
      readings[readIndex] = analogRead(A0);
      total = total + readings[readIndex];
      readIndex++;
      if (readIndex >= numReadings) {readIndex = 0; }
      average = total / numReadings; 
//traceur serie        
Serial.print(readings[readIndex]); Serial.print(",");Serial.print(average);Serial.println(","); 
/*
for (byte i=0; (i<numReadings);i++) {Serial.print(readings [i]);Serial.print(";");}    
Serial.println(",");
*/
      temps=0;     }                             
//digitalWrite(13,LOW);

}//fin routine


void loop() { 
 
//lcd.setCursor(0,3); lcd.print("average smooth"); lcd.print(average);lcd.print("    ");

} // fin loop

Vérification du fonctionnement du programme,
Voici quelques simulations pour vérifier les caractéristiques du filtre à moyenne mobile

Avec une période d’échantillonnage 0.005s et 10 échantillons donc période du signal à traiter au minimale 50ms (20Hz)
La précision de la valeur moyenne pour un signal carré sera de 0.1 par défaut, le temps de l’algo est de 0.2ms
Fc=(Fechantillon0.5)/nbr d’echantillon=(2000.5)/10=10Hz

D’ailleurs, un signal carré de 20Hz, avec un rapport cyclique de 25%, donnera seulement une valeur moyenne de 102420%=204 et pas 102425%=256

Si le signal est de 40Hz la précision passe à 0.2 par défaut, car le micro ne prend que 5 échantillons.

Pour un signal sinusoïdal avec un offset de 2.5V et une amplitude de 2.5V de 20Hz, la valeur moyenne est bien stable à 512 comme on peut l’observer sur la figure suivante.
Alors que si on avait utilisé un filtre passe bas à 20Hz, il y aurait une petite ondulation de la valeur moyenne autour de de 512.

Pour 60 échantillons et une fréquence d’échantillon de 200hz, la fréquence minimale de mesure sera de 3.33Hz
la fréquence de coupure Fc=(Fechantillon0.5)/nbr d’echantillon=(2000.5)/60=1,66Hz
Pour signal carré avec une fréquence de 5Hz, le nombre d’échantillon pris en compte sera de 40 donc la précision sera de 2.5% mais il y a une fluctuation autour de la valeur moyenne pour la fréquence de 5Hz.


Evidement avec l’ECG, le nombre d’échantillon est une valeur moyenne d’un nombre de détection de BPM sur 1 minute et ne correspond pas exactement à la valeur moyenne précédente mais s’en rapproche.
voici les résultats
La courbe bleu est le BPM, en rouge la valeur moyenne du BPM avec 60 échantillons, en vert l’index du tableau pour faire le moyennage avec un envoie des données tous les 6 secondes. L’incrémentation de l’index permet de visualiser qu’il y a bien un certain nombre de mesure du BPM dans un temps correspondant au battement et qu’il n’y a pas de perte du signal de detection.
Avec 60 echantillons et un BMP autour de 1hz, on peut observer que le « moyenneur glissant » filtre moyennement correctement car l’amplitude des variations est de ±5 BMP autour de la mesure.

Pour 120 échantillons ce qui prend 40% de la mémoire d’un Atmega 328, La courbe est plus lissé que la courbe precedente

La littérature parle de « valeur moyenne glissante pondérée linéairement » ou le poids de chaque terme décroit linéairement mais que le gain statique reste à 1 ??????

Mais quel est l’avantage de cette pondération par rapport à la valeur moyenne glissante ?

ECG_filtre_moyenneur.ino (4.35 KB)

9

Les thermes b pour la pondération de la valeur moyenne glissante, sont déterminés par l’équation suivante, leur somme doit être de 1 pour avoir un gain du filtre unitaire. 2M correspond au nombre d’échantillon.

Pour 10 échantillons voici les coefficients, on peut observer sur la courbe rouge de l’atténuation avec la moyenne pondérée linaire est un peu plus faible que l’atténuation avec la moyenne classique.
Donc pas viable pour notre application.

Il y a pléthore d’article sur l’ECG….mais souvent cela manque d’explication et de vulgarisation....
En voici 1 qui n’est pas trop compliqué et qui explique assez bien en utilisant la dérivée première et seconde sur signal de l'ECG en utilisant la convolution de 2 moyennes glissantes pondérées de 5 points pour la détection de l’onde R, donc du BPM
« Evolution temporelle et fréquentielle de l’ECG : analyse battement par battement »…




La précision de mesure avec la dérivée première et seconde n’est pas présentée….les échelles non plus
Mais, ce serait intéressant de tester….

10

This topic was automatically closed 120 days after the last reply. New replies are no longer allowed.