Pendule : résonance paramétrique

2.a. Schéma général

La boussole est un barreau aimanté de longueur 60 mm, comportant en son milieu une encoche qui lui permet de reposer de manière stable sur une pointe verticale. Le barreau peut donc tourner autour d'un axe vertical, ce qui en fait une boussole pouvant s'aligner avec un champ magnétique horizontal. Le champ magnétique est produit par deux bobines plates (bobines de Helmholtz). Chacune de ces bobines comporte 100 spires de rayon 60 mm (réparties sur 5 couches) occupant une largeur 25 mm. Les centres des bobines sont espacées de 90 mm, ce qui laisse une zone de longueur environ 60 mm entre les deux bobines, dans laquelle la boussole est placée.

Les deux bobines, branchées en série, sont alimentées par un pont de transistors DMOS dont le fonctionnement est décrit dans Onduleur diphasé MLI pour bobines. Le signal à modulation de largeur d'impulsion (MLI) qui permet de piloter le pont DMOS est délivré par une carte Arduino MEGA. L'alimentation de puissance délivre une tension de 12 V et peut débiter jusqu'à 5 A. Le programme de l'arduino (décrit plus loin) génère le signal MLI, qui permet d'obtenir dans les bobines un courant d'intensité :

Le champ magnétique entre les deux bobines (approximativement uniforme) a donc la forme suivante :

Soit J le moment d'inertie du barreau par rapport à son axe de rotation et m son moment magnétique. L'équation différentielle du mouvement est :

qui est bien de la forme $\mathrm{(2)}$ si on pose : ${\omega }_0^2=\frac{m{B}_1}{J}$ et $h=\frac{B_2}{B_1}=\frac{I_2}{I_1}$. Le paramètre h est donc facilement contrôlable par le programme arduino mais la pulsation propre ${\omega }_0$ doit être déterminée expérimentalement, en analysant les petites oscillations lorsque B₂=0.

L'enregistrement du mouvement de rotation de la boussole se fait avec une caméra placée au dessus du dispositif, dont l'axe optique coïncide avec l'axe de rotation. Il s'agit d'une caméra Basler USB (modèle acA1300-200uc) capable de fournir 200 images par secondes en pleine résolution, et plus à une résolution plus basse. Dans cette expérience, la cadence est fixée à 200 images par secondes. La caméra est munie d'une prise de déclenchement externe, qui permet de piloter l'acquisition des images depuis la carte arduino, afin de synchroniser la vidéo avec la variation du champ magnétique, donc de connaître l'état exact du champ magnétique pour chaque image de la vidéo.

2.b. Programme arduino

Le programme (arduino MEGA ou UNO) génère le signal MLI et le signal de déclenchement de la caméra. Il permet en fait de générer deux signaux de commande (pour deux onduleurs) bien que cette application n'en utilise qu'un seul.

 
#include "Arduino.h"
#define NECHANT 128 // nombre d'échantillons de la table du signal
#define SHIFT_ACCUM 25
#define SET_FREQ 100
#define SET_AMP 101
#define START_SYNCHRO 102
#define STOP_SYNCHRO 103
#define PERIOD 100 // période du signal MLI en microsecondes (10 kHz)
#define NSYNCHRO 25 // période du signal de synchro = 2*NSYNCHRO*PERIOD

uint32_t icr;
uint32_t table_onde_A[NECHANT];
uint32_t table_onde_B[NECHANT];
uint32_t accumA,accumB,increm;
uint16_t diviseur[6] = {0,1,8,64,256,1024};
volatile uint8_t flag;
uint8_t strob = 0;
uint32_t indexA,indexB;
uint32_t compteur;
uint8_t start_synchro, synchro;

    			

Les deux signaux MLI (signaux A et B) sont générés au moyen du Timer 1. La fonction suivante configure et déclenche ce timer, avec une période donnée en microsecondes.

 
void init_pwm_timer1(uint32_t period) {
    TCCR1A = 0;
    TCCR1A |= (1 << COM1A1); //Clear OCnA/OCnB/OCnC on compare match, set OCnA/OCnB/OCnC at BOTTOM (non-inverting mode)
    TCCR1A |= (1 << COM1B1);
#if defined(__AVR_ATmega2560__) || defined(__AVR_ATmega32U4__)
    TCCR1A |= (1 << COM1C1);
#endif
    TCCR1B = 1 << WGM13; // phase and frequency correct pwm mode, top = ICR1
    int d = 1;
    icr = (F_CPU/1000000*period/2);
    while ((icr>0xFFFF)&&(d<6)) { // choix du diviseur d'horloge
        d++;
        icr = (F_CPU/1000000*period/2/diviseur[d]);
   } 
   ICR1 = icr; // valeur maximale du compteur
   TIMSK1 = 1 << TOIE1; // overflow interrupt enable
   sei(); // activation des interruptions
   TCNT1 = 0; // mise à zéro du compteur
   TCCR1B |= d; // déclenchement du compteur
}
    			

Les deux signaux sont générés sur les sorties suivantes :

• Arduino MEGA : D11 pour le signal A et D12 pour le signal B.
• Arduino UNO ou YUN : D9 pour le signal A et D10 pour le signal B.

Le timer génère une interruption à chaque dépassement de la valeur du compteur, c'est-à-dire à chaque période. Le gestionnaire d'interruption, donné ci-dessous effectue l'incrémentation des deux accumulateurs de phase et modifie le rapport cyclique des deux signaux en modifiant les registres OCR1A et OCR1B avec les valeurs stockées dans les tableaux table_onde_A et table_onde_B.

Lorsque l'indice de la première table s'annule alors que la variable start_synchro est égale à 1, la génération du signal de synchronisation (pour la caméra) est déclenchée par la mise à 1 de la variable synchro. La génération de ce signal consiste à incrémenter une variable compteur jusqu'à la valeur NSYNCHRO. À chaque fois que ce compteur revient à zéro, on fait basculer la sortie 7 (dont l'état est mémorisé par la variable flag). En conséquence, un signal de période 2*NSYNCHRO*PERIOD est généré sur la sortie D7. Ce signal sert à déclencher l'enregistrement des images par la caméra. Dans ce code, la période est 2*25*100=5000 us, soit une fréquence de 200 Hz.

 
ISR(TIMER1_OVF_vect) { // Timer 1 Overflow interrupt
  accumA += increm;
  accumB += increm;
  indexA = accumA >> SHIFT_ACCUM;
  OCR1A = table_onde_A[indexA];
  indexB = accumB >> SHIFT_ACCUM;
  OCR1B = table_onde_B[indexB];
  if ((indexA==0)&&(start_synchro==1)) {
    compteur = 0; // première image au début du cycle (champ A nul et croissant)
    flag = 0;
    start_synchro = 0;
    synchro = 1;
  }
  if (synchro) {
    if (compteur==0) {
      
      if (flag) {
    #if defined(__AVR_ATmega2560__)
        PORTH &= ~(1 << PORTH4); 
    #else
        PORTD &= ~(1 << PORTD7); 
    #endif   
        flag = 0;
      }
      else {
   #if defined(__AVR_ATmega2560__)
        PORTH |= (1 << PORTH4);  
   #else
        PORTD |= (1 << PORTD7);      
   #endif
        flag = 1;
      }
    }
    compteur += 1;
    if (compteur==NSYNCHRO) compteur=0;
  }
}
    			  

La fonction suivante remplit le tableau pour le signal A avec les échantillons de valeur de OCR1A (seuil de bascule de la sortie) qui permettent d'obtenir une modulation MLI sinusoïdale, avec une amplitude et un décalage donné. Ces deux valeurs sont comprises entre 0 et 1 et l'amplitude plus le décalage ne doit pas dépasser 1.

 
void set_sinus_table_A(float amp, float offset) {
  int i;
  float dt = 2*PI/NECHANT;
  for(i=0; i<NECHANT; i++) {
    table_onde_A[i] = icr*0.5*(1.0+offset+amp*sin(i*dt));
  }  
}    			  
    			  

La fonction suivante fait la même chose pour le tableau du signal B.

 
void set_sinus_table_B(float amp, float offset) {
  int i;
  float dt = 2*PI/NECHANT;
  for(i=0; i<NECHANT; i++) {
    table_onde_B[i] = icr*0.5*(1.0+offset+amp*sin(i*dt));
  }  
}    			  
    			  

La fonction setup configure les sorties, initialise les variables et configure le Timer. Une fréquence de modulation par défaut de 5 Hz est appliquée. La phase de chaque signal MLI est contrôlée par un accumulateur de phase (une phase donnée par un nombre 32 bits à virgule fixe). L'accumulateur du signal B a une valeur initiale qui permet d'obtenir un signal B en quadrature par rapport au signal A (mais le signal B n'est pas utilisé dans cette expérience).

void setup() {
    char c;
    Serial.begin(115200);
    Serial.setTimeout(0);
    c = 0;
    Serial.write(c);
    c = 255;
    Serial.write(c);
    c = 0;
    Serial.write(c);
    pinMode(8,OUTPUT); 
    digitalWrite(8,HIGH); // ENABLE A
    pinMode(6,OUTPUT);
    digitalWrite(6,HIGH); // ENABLE B
#if defined(__AVR_ATmega2560__) // arduino MEGA
    pinMode(11,OUTPUT); // INA
    pinMode(12,OUTPUT); // INB
#else // arduino UNO et YUN
    pinMode(9,OUTPUT); // INA
    pinMode(10,OUTPUT); // INB
#endif
	pinMode(7,OUTPUT); // sortie synchro pour la caméra (port PD7 sur UNO, PH4 sur MEGA)
    digitalWrite(7,LOW);
    
    float frequence = 5; // en Hz
    accumA = 0;
    accumB = ((uint32_t)(NECHANT * 0.25)) << SHIFT_ACCUM;
    increm = (uint32_t) (((float)(0xFFFFFFFF))*((float)(frequence)*1e-6*(float)(PERIOD))); // incrément de l'accumulateur de phase
    flag=0;
    compteur = 0;
    start_synchro = 0;
    synchro = 0;
    init_pwm_timer1(PERIOD);
    set_sinus_table_A(0.2,0.5);
}
    			    

La fonction loop gère la communication série avec le PC, afin de traiter les ordres de changement de fréquence, de changement d'amplitude et d'offset du signal A, de démarrage et d'arrêt du signal de synchronisation de la caméra.

  void loop() {
  char com;
  uint32_t c1,c2,c3,c4;
  uint32_t freq; // fréquence en mHz
  float frequence; // fréquence en Hz
  uint16_t A; // amplitude
  uint16_t Off; // offset
  float amp,offset;
    if (Serial.available()>0) {
       
          com = Serial.read();
          if (com==SET_FREQ) {
              while (Serial.available()<4) {};
              c1 = Serial.read();
              c2 = Serial.read();
              c3 = Serial.read();
              c4 = Serial.read();
              freq = ((c1 << 24) | (c2 << 16) | (c3 << 8) | c4);
              frequence = (float)(freq);
              frequence /= 1000.0;
              increm = (uint32_t) (((float)(0xFFFFFFFF))*((float)(frequence)*1e-6*(float)(PERIOD))); // incrément de l'accumulateur de phase
          }
          else if (com==SET_AMP) {
              while (Serial.available()<2) {};
              c1 = Serial.read(); 
              c2 = Serial.read();
              A = (c1 << 8) | c2;
              amp = (float) A;
              while (Serial.available()<2) {};
              c1 = Serial.read();
              c2 = Serial.read();
              Off = (c1 << 8) | c2;
              offset = (float) Off;
              set_sinus_table_A(amp/100,offset/100);
          }
          else if (com==START_SYNCHRO) {
             start_synchro = 1;
            
          }
          else if (com==STOP_SYNCHRO) {
             synchro=0;
             start_synchro = 0;
             digitalWrite(7,LOW);
          }
       
    }
}  			     
    			     

2.c. Programme de commande

La configuration du signal MLI et le déclenchement du signal de synchronisation sont effectués depuis un script python (par liaison série avec l'arduino). La classe ci-dessous comporte les fonctions permettant d'envoyer des informations ou des ordres au programme arduino :

import serial
import numpy
    			 
class Arduino():
    def __init__(self,port):
        self.ser = serial.Serial(port,baudrate=115200)
        c_recu = self.ser.read(1)
        while ord(c_recu)!=0:
            c_recu = self.ser.read(1)
        c_recu = self.ser.read(1)
        while ord(c_recu)!=255:
            c_recu = self.ser.read(1)
        c_recu = self.ser.read(1)
        while ord(c_recu)!=0:
            c_recu = self.ser.read(1)
        self.SET_FREQ = 100
        self.SET_AMP = 101
        self.START_SYNCHRO = 102
        self.STOP_SYNCHRO = 103

    def close(self):
        self.ser.close()

    def write_int8(self,v):
        self.ser.write((v).to_bytes(1,byteorder='big'))

    def write_int16(self,v):
        v = numpy.int16(v)
        char1 = int((v & 0xFF00) >> 8)
        char2 = int((v & 0x00FF))
        self.ser.write((char1).to_bytes(1,byteorder='big'))
        self.ser.write((char2).to_bytes(1,byteorder='big'))
        
    def write_int32(self,v):
        v = numpy.int32(v)
        char1 = int((v & 0xFF000000) >> 24)
        char2 = int((v & 0x00FF0000) >> 16)
        char3 = int((v & 0x0000FF00) >> 8)
        char4 = int((v & 0x000000FF))
        self.ser.write((char1).to_bytes(1,byteorder='big'))
        self.ser.write((char2).to_bytes(1,byteorder='big'))
        self.ser.write((char3).to_bytes(1,byteorder='big'))
        self.ser.write((char4).to_bytes(1,byteorder='big'))

    def set_frequence(self,frequence):
        self.ser.write((self.SET_FREQ).to_bytes(1,byteorder='big'))
        self.write_int32(frequence*1000)
    def set_amplitude(self,amp,offset):
        self.ser.write((self.SET_AMP).to_bytes(1,byteorder='big'))
        self.write_int16(amp*100)
        self.write_int16(offset*100)
    def stop_synchro(self):
        self.write_int8(self.STOP_SYNCHRO)
    def start_synchro(self):
        self.write_int8(self.START_SYNCHRO)
    			 
    			 

Le code suivant configure la fréquence, l'amplitude et l'offset de la modulation MLI (la fréquence de hachage est codée en dur dans le programme arduino), puis exécute une boucle qui permet à l'utilisateur de démarrer puis de stopper l'émission du signal de synchronisation (ou de sortir avec le caractère S).

ard = Arduino("COM18")
ard.set_frequence(5)
ard.set_amplitude(0.5,0.0)
ard.stop_synchro()
synchro = False
while True:
    r = input("?")
    if r=="S":
        break
    if synchro:
        ard.stop_synchro()
        synchro = False
        print("stop synchro")
    else:
        ard.start_synchro()
        synchro = True
        print("start synchro")   			 
    			 

La configuration de la caméra est faite dans le logiciel Basler Pylon Viewer. Elle consiste à ajuster le cadre de l'image afin que seule la zone d'intérêt soit enregistrée (dans le cas présent, une taille de 600 par 600 pixels environ convient) puis à ajuster le temps de pose par image à une valeur inférieure à la période d'échantillonnage (5 ms). L'utilisation d'un éclairage à LED homogène et assez puissant est nécessaire pour permettre un temps de pose aussi court. Il faut aussi configurer un déclenchement externe de l'acquisition de chaque image, qui est effectué par le signal de synchronisation délivré par l'arduino. La vidéo est enregistrée dans un fichier AVI avec une compression IYUV. Le nombre d'images de la vidéo doit être choisi en fonction de la durée voulue. Par exemple, un enregistrement de 1000 images permet d'obtenir une durée de 1000/200=5 s, puisque la fréquence d'acquisition est de 200 images par seconde.

2.d. Test de l'onduleur et étalonnage

Une sonde à effet Hall est placée sur l'axe des bobines (sans la boussole) afin de mesurer la composante B_z du champ magnétique, soit au centre (z=0), soit en un point correspondant à peu près à l'extrémité du barreau aimanté (z=30 mm).

Voici les enregistrements du champ magnétique pour une tension d'alimentation de 12,0 V et une amplitude de modulation MLI de 0,5 (sans offset), soit la moitié de l'amplitude maximale. L'acquisition est déclenchée avec le signal de synchronisation.

Le champ magnétique est nécessairement en phase avec le courant électrique dans les bobines, donc en phase avec le signal de modulation. Sur ces enregistrements, on constate que le signal du champ magnétique est légèrement en retard par rapport au signal MLI de commande, très probablement à cause du retard introduit par l'amplificateur de la sonde à effet Hall. Pour cette raison, la sonde à effet Hall ne peut être utilisée pour effectuer une mesure du champ magnétique synchronisée avec la prise de vue vidéo. Ces courbes permettent néanmoins de vérifier la période de variation du champ (ici 200 ms) et fournissent son amplitude. La variation du champ entre le centre du barreau et son extrémité est très faible. Pour cette amplitude de modulation de 0,5, la valeur efficace de B_z est de 1,3 mT au centre et de 1,4 mT à l'extrémité. On retiendra donc une valeur moyenne de champ efficace de 1,35 mT, soit une amplitude de 1,9 mT pour une amplitude de modulation de 0,5. Cette valeur est suffisante pour cette expérience mais il est possible de l'augmenter en augmentant la tension d'alimentation. Par exemple, une tension de 24 V devrait en principe conduire à un courant deux fois plus grand donc un champ magnétique deux fois plus grand.

2.e. Traitement vidéo : tracé du champ magnétique

Ce traitement consiste à ajouter sur la vidéo une flèche montrant le champ magnétique, d'afficher le temps et la valeur du champ magnétique et de produuire un fichier vidéo avec compression. On utilise le fait que la première image de la vidéo coïncide avec le début d'un cycle de modulation en sinus, c'est-à-dire que le champ est nul à cet instant. On doit aussi tenir compte de l'orientation de la caméra par rapport à la direction du champ lorsque le courant est positif. On commence par appliquer un courant positif constant (en choissant une amplitude nulle et un offset positif), ce qui a pour effet d'aligner le barreau avec l'axe des bobines. La caméra est alors orientée afin que le côté vertical de l'image soit parallèle à cette direction, tout en s'assurant que son axe reste bien coïncidant avec l'axe de rotation du barreau. On peut alors repérer sur l'image la direction du champ magnétique lorsque le rapport cyclique du signal MLI est positif. Ce réglage est délicat et nécessite donc un support de caméra flexible et facilement réglable.

Le script python utilise la bibliothèque de traitement d'images OpenCV. La vidéo produite est compressée avec l'encodage XVID (implémentation libre de la norme MPEG-4).

import cv2 as cv
import numpy

nom = 'B-sinus-12V-amp0,5-4,5Hz-200fps-2'
video = cv.VideoCapture('%s.avi'%nom)

freq=4.5
amp=0.5
offset=0
echelle_B = 1.9/0.5

cv.namedWindow("frame")
def onMouse(event,x,y,flags,param):
    global i0,j0
    if event==cv.EVENT_LBUTTONUP:
        i0=x
        j0=y
        
cv.setMouseCallback('frame',onMouse)
success,frame=video.read()
cv.imshow('frame',frame)
cv.waitKey()
(h,w,c) = frame.shape

fourcc = cv.VideoWriter_fourcc(*'XVID')
out = cv.VideoWriter('%s-champB.mp4'%nom, fourcc, 30.0, (w,h))

n=0
scale=400
t=0
nombre_images = 500

while success and n<nombre_images:
    n += 1
    B=offset+amp*numpy.sin(2*numpy.pi*freq*t)
    t += 1.0/200
    j1 = int(j0-B*scale)
    s = 1
    if B<0 : s=-1
    frame = cv.line(frame,(i0,j0),(i0,j1),(1,0,0),2)
    frame = cv.line(frame,(i0,j1),(i0-10,j1+s*10),(1,0,0),2)
    frame = cv.line(frame,(i0,j1),(i0+10,j1+s*10),(1,0,0),2)
    frame = cv.putText(frame,"t = %0.2f s"%t,(30,40),cv.FONT_HERSHEY_SIMPLEX,1,(1,0,0))
    Bz = B*echelle_B
    frame = cv.putText(frame,"Bz = %0.2f mT"%Bz,(30,80),cv.FONT_HERSHEY_SIMPLEX,1,(1,0,0))
    out.write(frame)
    cv.imshow('frame',frame)
    if cv.waitKey(10) == ord('q'):
        break
    success,frame=video.read()
    
video.release()
out.release()   	      
    			      

Au démarrage, il apparaît une fenêtre avec la première image. Il faut alors cliquer à la souris (bouton gauche) sur le centre de la boussole pour marquer le point origine du vecteur B, puis appuyer sur Entrée pour débuter la parcours de la vidéo. La boucle parcourt un nombre d'images déterminé.

Le fichier vidéo MP4 résultant (avec l'encodage XVID) peut être converti en fichier WEBM pour le web (encodage vp9) avec ffmpeg au moyen de la commande suivante :

ffmpeg -i video.mp4 -b:v 500k video.webm

L'option -b:v 500k permet de préciser le débit pour la vidéo produite, ici 500 kbits/s.

L'encodage vp9 permet d'obtenir une très bonne qualité pour un débit relativement faible, dans le but de diffuser la vidéo sur le web. Ce fichier vidéo peut être directement inséré dans une page HTML (via l'élément HTML5 video).

Il est possible de convertir seulement une partie de la vidéo, par exemple les images 0 à 500 de la manière suivante :

ffmpeg -i video.mp4 -vf select='between(n,\0,\500)' -b:v 500k video.webm

2.f. Traitement vidéo : extraction de l'angle

Voici à titre d'exemple une vidéo montrant le mouvement d'oscillation de la boussole dans un champ magnétique constant. Le rapport cyclique du signal MLI est fixé à 0,5, ce qui correspond à un champ magnétique de 1,9 mT. La vidéo est lue à 30 images par secondes, soit un ralenti d'un facteur 200/30.

On souhaite mettre en place un traitement d'image permettant d'obtenir l'angle entre le barreau et la direction du champ magnétique (parallèle au côté vertical de l'image). La couleur rouge très marquée d'une moitié du barreau nous oriente vers une extraction à partir de la couleur, par la méthode de rétroprojection d'histogramme (voir Suivi d'un objet coloré dans une vidéo).

L'algorithme comporte les étapes suivantes :

• Sélection d'un rectangle délimitant la région rouge du barreau.
• Conversion de l'image en HSV et calcul de l'histogramme HS de la région sélectionnée.
• Début de la boucle de traitement de la vidéo.
• Pour chaque image de la vidéo, conversion en HSV et rétroprojection de l'histogramme sur l'image, suivi d'un seuillage pour extraire les pixels rouges.
• Pour obtenir un ensemble de pixels connexes, dilatation puis érosion.
• Extraction de la liste des pixels de la région, par la technique du coloriage.
• Sélection des P pixels les plus éloignés du centre et calcul de leur barycentre.
• Calcul de l'angle à partir des positions du barycentre et du centre.

import cv2 as cv
import numpy
import matplotlib.pyplot as plt
import sys
sys.setrecursionlimit(10000)
                 

Après extraction des pixels rouges, on obtient une image en niveaux de gris dont tous les pixels sont noirs (valeur 0) sauf les pixels détectés qui sont blancs (valeur 255). Les pixels associés à la zone rouge du barreau sont supposés connexes. La fonction suivante effectue une extraction de la liste des pixels de l'objet par la technique du coloriage, consistant à colorier chaque pixel détecté avec la valeur 100 puis à explorer récursivement ses voisins. Le coloriage se fait en partant d'un pixel appartenant à l'objet (ici la partie rouge du barreau) et permet de sélectionner uniquement les pixels appartenant à cet objet, au cas où des pixels situés en dehors de l'objet auraient été sélectionnés par le seuillage.

def test_pixel_objet(A,M,N,i,j,liste_pixels):
    liste_pixels.append([i,j])
    A[i,j]=100
    voisins=[(i+1,j),(i-1,j),(i,j-1),(i,j+1)]
    for (k,l) in voisins:
        if k>=0 and k<M and l>=0 and l<N:
            if A[k,l]==255:
                test_pixel_objet(A,M,N,k,l,liste_pixels)                     
                   

La fonction suivante recherche le pixel blanc le plus proche d'un pixel donné. Elle permet de rechercher un pixel de l'objet à partir d'une position qui est proche de l'objet sans être dans l'objet (avant d'appliquer le coloriage).

def recherche_voisinage(A,M,N,i0,j0):
    if A[i0,j0]==255:
        return([i0,j0])
    q=0
    while True:
        q+=1
        for i in range(i0-q,i0+q+1):
            j=j0-q
            if i>=0 and i<M and j>=0 and j<N and A[i,j]==255:
                return([i,j])
            j=j0+q
            if i>=0 and i<M and j>=0 and j<N and A[i,j]==255:
                return([i,j])
        for j in range(j0-q+1,j0+q):
            i=i0-q
            if i>=0 and i<M and j>=0 and j<N and A[i,j]==255:
                return([i,j])
            i=i0+q
            if i>=0 and i<M and j>=0 and j<N and A[i,j]==255:
                return([i,j])                  
                  

La fonction suivante calcule le barycentre d'un liste de pixels :

def barycentre(liste_pixels):
    N=len(liste_pixels)
    xG=0.0
    yG=0.0
    for pixel in liste_pixels:
        xG += pixel[1]
        yG += pixel[0]
    xG /= N
    yG /= N
    return ([xG,yG])                 
                  

La fonction suivante extrait d'une liste de pixels les P pixels les plus éloignés du centre (dont les coordonnées sont données par ic,jc), puis renvoie leur barycentre.

def barycentre_points_extremite(liste_pixels,ic,jc,P):
    N=len(liste_pixels)
    liste_pixels_distance = []
    for k in range(N):
        i = liste_pixels[k][0]
        j = liste_pixels[k][1]
        d2 = (i-ic)**2+(j-jc)**2
        liste_pixels_distance.append((i,j,d2))
    liste_pixels_distance = sorted(liste_pixels_distance, key=lambda point:point[2])
    liste_pixels_extrem = liste_pixels_distance[N-P:N]
    return barycentre(liste_pixels_extrem)                 
                  

On commence par afficher les 50 premières de la vidéo afin de permettre à l'utilisateur de sélectionner (une ou plusieurs fois) un rectangle inscrit dans la zone rouge. Les rectangles sélectionnés sont convertis au format de couleur HSV et stockés dans une liste. Dans le cas présent, la sélection d'un seul rectangle (coïncidant avec la partie rouge du barreau) suffit.

nom = 'B-constant-12V-amp0,5-200fps-2'
nom_video = "%s.avi"%nom
techant = 1/200
video = cv.VideoCapture(nom_video)
success,frame=video.read()
liste_hsv_region = []
n=0
while success and n<50:
    r=cv.selectROI(frame)
    if r!=(0,0,0,0):
        print(r)
        region = frame[int(r[1]):int(r[1]+r[3]), int(r[0]):int(r[0]+r[2])]
        hsv_region = cv.cvtColor(region,cv.COLOR_BGR2HSV)
        liste_hsv_region.append(hsv_region)
    success,frame=video.read()
    n += 1          
                  

En deuxième étape, une image de la vidéo est affichée et l'utilisateur doit cliquer deux fois avec le bouton gauche de la souris, une première fois pour pointer le centre du barreau (c.a.d. l'axe de rotation) et une seconde fois pour sélectionner un pixel situé dans la zone rouge, après quoi il ferme la fenêtre en appuyant sur une touche quelconque.

click = 0
def onMouse(event,x,y,flags,param):
    global i0,j0,ic,jc,click
    if event==cv.EVENT_LBUTTONUP:
        if click==0:
            ic=y
            jc=x
            print(ic,jc)
            click = 1
        else :
            i0=y
            j0=x
            click = 2
            print(i0,j0)
cv.namedWindow("frame")       
cv.setMouseCallback('frame',onMouse)
video = cv.VideoCapture(nom_video)
success,frame=video.read()
(h,w,c) = frame.shape
cv.imshow('frame',frame)
cv.waitKey()                    
                    

L'histogramme HS (Hue-Saturation) est calculé à partir de la liste des rectangles selectionnés précédemment, puis affiché.

hist1 = cv.calcHist(liste_hsv_region,[0,1],None,[30,30],[0,180,0,255])
cv.normalize(hist1, hist1, alpha=0, beta=255, norm_type=cv.NORM_MINMAX)
plt.matshow(hist1,extent=[0,180,0,255],cmap=plt.cm.gray)
plt.xlabel("Hue")
plt.ylabel("Sat")
plt.show()                  
                    

L'étape suivante consiste à tester et à régler l'extraction de la zone de couleur rouge. La rétroprojection de l'histogramme (réalisée par la fonction calcBackProject) consiste à créer une image en niveaux de gris dont la valeur de chaque pixel est égale à celle de l'histogramme HS pour les valeurs (H,S) du pixel correspondant sur l'image en couleur. Cette opération est suivie d'un seuillage qui produit une image où tous les pixels sont noirs sauf les pixels rouges qui sont blancs. À ce stade, les pixels sélectionnés ne forment pas une région connexe, ce qui rend impossible l'application du coloriage. On procède donc à une dilatation suivie d'une érosion.

hsv = cv.cvtColor(frame,cv.COLOR_BGR2HSV)
back_proj = cv.calcBackProject([hsv],[0,1],hist1,[0,180,0,255],scale=1)
cv.imshow('frame',back_proj)
r,img = cv.threshold(back_proj,100,255,cv.THRESH_BINARY)
element = cv.getStructuringElement(cv.MORPH_RECT,(20,20))
img = cv.dilate(img,element)
img = cv.erode(img,element)
cv.imshow('frame',img)                   
                    

Voici enfin la boucle de traitement de la vidéo. La boucle s'arrête lorsque le nombre d'images indiqué est atteint ou lorsque l'utilisateur appuie sur q. On enregistre une nouvelle vidéo contenant les pixels coloriés, le barycentre des points de l'extrémité et le segment reliant ce point au centre de rotation. Le temps et l'angle sont aussi affichés.

fourcc = cv.VideoWriter_fourcc(*'XVID')
out = cv.VideoWriter('%s-angle.mp4'%nom, fourcc, 30.0, (w,h))
n= 0
t=0
liste_theta = []
while success and n<950 :
    n += 1
    success,frame=video.read()
    (M,N,c)=frame.shape
    hsv = cv.cvtColor(frame,cv.COLOR_BGR2HSV)
    back_proj = cv.calcBackProject([hsv],[0,1],hist1,[0,180,0,255],scale=1)
    r,img = cv.threshold(back_proj,100,255,cv.THRESH_BINARY)
    element = cv.getStructuringElement(cv.MORPH_RECT,(20,20))
    img = cv.dilate(img,element)
    img = cv.erode(img,element)
    [i0,j0] = recherche_voisinage(img,M,N,i0,j0)
    liste_pixels=[]
    test_pixel_objet(img,M,N,i0,j0,liste_pixels) # extraction d'une région de pixels connexes
    
    P=500
    (xG,yG) = barycentre_points_extremite(liste_pixels,ic,jc,P)
    i2 = yG
    j2 = xG
    X = j2-jc
    Y = ic-i2
    Z = Y+1j*X
    a = numpy.angle(Z)
    liste_theta.append(a)
    
    cv.circle(img,(int(j2),int(i2)),5,(255,0,0),2)
    cv.line(img,(int(j2),int(i2)),(jc,ic),(255,0,0),2)
    cv.putText(img,"a = %0.3f rad"%a,(40,550),cv.FONT_HERSHEY_SIMPLEX,1,(255,0,0))
    t += techant
    cv.putText(img,"t = %0.2f s"%t,(40,500),cv.FONT_HERSHEY_SIMPLEX,1,(255,0,0))
    
    frame[:,:,0] = img
    frame[:,:,1] = img
    frame[:,:,2] = img
    out.write(frame)
    cv.imshow('frame',img)
    if cv.waitKey(1) == ord('q'):
        break
    
video.release()
out.release()          
			

Pour finir, on trace l'angle en fonction du temps puis on enregistre les valeurs dans un fichier texte.

theta = numpy.array(liste_theta)
N = len(theta)
t = numpy.arange(N)*techant
plt.figure()
plt.plot(t,theta,"ko-")
plt.xlabel("t (s)")
plt.ylabel("theta (rad)")
plt.grid()
nom_fichier = "%s-angle.txt"%nom
numpy.savetxt(nom_fichier,numpy.array([t,theta]).T,delimiter="\t",fmt="%.4e",header="t(s)\t theta(rad)")
plt.show()			
			

Voici le tracé de l'angle en fonction du temps :

[t,a] = numpy.loadtxt('B-constant-12V-amp0,5-200fps-2-angle.txt',unpack=True,skiprows=1)
figure(figsize=(12,6))
plot(t,a,'-')
xlabel(r"$t\ (\rm s)$",fontsize=18)
ylabel(r"$\theta\ (\rm rad)$",fontsize=18)
grid()

			

On constate que l'angle moyen n'est pas nul, sans doute parce que le champ magnétique n'est pas tout à fait aligné avec le bord vertical de l'image (ou à cause de l'incertitude de pointé du centre).

Voici la vidéo générée par l'analyse :

Une analyse spectrale permet de déterminer la fréquence de ces petites oscillations pour un champ statique de 1,9 mT :

import numpy.fft
import scipy.signal 

def spectre(t,u):
    N=len(u)
    te=t[1]-t[0]
    zeros=numpy.zeros(6*N)
    U = numpy.concatenate((u*scipy.signal.blackman(N),zeros))
    NN=len(U)
    spectre = numpy.absolute(numpy.fft.fft(U))*2.0/N/0.42
    freq = numpy.arange(NN)*1.0/(NN*te)
    return (freq,spectre)

(freq,sp) = spectre(t,a)
figure()
plot(freq,sp)
xlim(0,10)
grid()
            

La fréquence d'oscillation (c'est-à-dire ${\omega }_0/\left(2\pi \right)$) est donc 4,5 Hz pour un champ de 1,9 mT. On s'attend donc à une première instabilité pour une excitation paramétrique de fréquence 9,0 Hz, une deuxième pour une fréquence de 4,5 Hz, une troisième pour une fréquence de 2,25 Hz.

3.a. Résultats expérimentaux pour n=1

Pour un champ constant obtenu avec un rapport cyclique de 0,5 (soit un champ magnétique 1,9 mT), nous obtenons une fréquence de petites oscillations de 4,5 Hz. On commence par chercher l'instabilité paramétrique pour une fréquence d'oscillation du champ magnétique de 9,0 Hz, avec une composante statique de 1,9 mT. La boussole est initialement dans sa position d'équilibre relative à la partie statique du champ. Si on ne la touche pas, elle reste à l'équilibre jusqu'à une amplitude d'oscillation de l'ordre de 0,2 (h=0,4). Si au contraire on la dévie très légèrement de sa position d'équilibre (à la main), l'instabilité paramétrique se déclenche pour une amplitude de l'ordre de 0,03 (h 0,06). Voici la vidéo obtenue pour une amplitude de 0,09 (h=0,18) :

Voici l'angle en fonction du temps :

[t,a] = numpy.loadtxt('B-sinus-12V-amp0,09-offset0,5-9Hz-200fps-1-angle.txt',unpack=True,skiprows=1)
figure(figsize=(12,6))
plot(t,a)
xlabel(r"$t\ (\rm s)$",fontsize=18)
ylabel(r"$\theta\ (\rm rad)$",fontsize=18)
ylim(-0.9,0.9)
title("A = 0,18   f = 9Hz")
grid()

			

Voici le démarrage de l'instabilité, avec le tracé du champ magnétique :

h=0.18
f=9.0
Bz=1.9*(1+h*numpy.sin(2*numpy.pi*f*t))
figure(figsize=(12,6))
subplot(211)
plot(t,Bz)
ylabel(r"$B_z\ (\rm mT)$",fontsize=18)
xlim(1,4)
ylim(0,3)
grid()
subplot(212)   
plot(t,a)
xlabel(r"$t\ (\rm s)$",fontsize=18)
ylabel(r"$\theta\ (\rm rad)$",fontsize=18)
xlim(1,4)
ylim(-0.2,0.2)
grid()

			

Les oscillations qui apparaissent sous l'effet de l'excitation paramétrique se font à la fréquence propre (4,5 Hz), conformément à ce que prévoit la théorie pour un oscillateur harmonique.

Faisons une analyse spectrale :

(freq,sp) = spectre(t,a)
figure()
plot(freq,sp)
ylabel(r"$A\ (\rm rad)$",fontsize=18)
xlabel(r"$f\ (\rm Hz)$",fontsize=18)
xlim(0,15)
grid()
            

Cette analyse confirme que les oscillations se font bien à la fréquence propre. La présence d'une seconde raie de fréquence un peu plus basse est reliée au phénomènre de battement observé sur la représentation temporelle.

Voici le tracé de l'angle en fonction du temps pour une amplitude de 0,03, correspondant à h=0,06.

[t,a] = numpy.loadtxt('B-sinus-12V-amp0,03-offset0,5-9Hz-200fps-1-angle.txt',unpack=True,skiprows=1)
figure(figsize=(12,6))
plot(t,a)
xlabel(r"$t\ (\rm s)$",fontsize=18)
ylabel(r"$\theta\ (\rm rad)$",fontsize=18)
ylim(-0.9,0.9)
title("h = 0,06   f = 9Hz")
grid()

			

Pour cette valeur de h plus petite, la croissance des oscillations est plus lente et les battements sont plus lents.

3.b. Résultats expérimentaux pour n=2

La fréquence d'excitation est égale à la fréquence propre, soit 4,5 Hz. Dans ce cas, le seuil de déclenchement de l'instabilité est plus élevé. Voici le résultat pour h=0,2 :

[t,a] = numpy.loadtxt('B-sinus-12V-amp0,2-offset0,5-4,5Hz-200fps-1-angle.txt',unpack=True,skiprows=1)
figure(figsize=(12,6))
plot(t,a)
xlabel(r"$t\ (\rm s)$",fontsize=18)
ylabel(r"$\theta\ (\rm rad)$",fontsize=18)
grid()
title("h = 0,4  f=4,5 Hz")
ylim(-0.9,0.9)

			

Pour une même amplitude d'excitation paramétrique, l'amplitude des oscillations est plus faible, la croissance et les modulations sont plus lentes.

h=0.18
f=4.5
Bz=1.9*(1+h*numpy.sin(2*numpy.pi*f*t))
figure(figsize=(12,6))
subplot(211)
plot(t,Bz)
ylabel(r"$B_z\ (\rm mT)$",fontsize=18)
xlim(0,2)
ylim(0,3)
grid()
subplot(212)   
plot(t,a)
xlabel(r"$t\ (\rm s)$",fontsize=18)
ylabel(r"$\theta\ (\rm rad)$",fontsize=18)
xlim(0,2)
grid()

			

Voici le spectre

(freq,sp) = spectre(t,a)
figure()
plot(freq,sp)
ylabel(r"$A\ (\rm rad)$",fontsize=18)
xlabel(r"$f\ (\rm Hz)$",fontsize=18)
xlim(0,15)
grid()
            

4.a. Oscillateur harmonique

Voyons tout d'abord le cas d'une excitation paramétrique d'un oscillateur harmonique, qui correspond au cas où l'angle reste petit. L'équation différentielle considérée est :

Le coefficient K permet de représenter un frottement de type fluide. Voici le calcul pour une excitation à une fréquence double de la fréquence propre :

from scipy.integrate import odeint   			

w2=(2*numpy.pi)**2
omega = 2*numpy.pi*2
k=2e-2
h=0.2
def systeme(y,t):
    return [y[1],-w2*(1+h*numpy.cos(omega*t))*y[0]-k*y[1]]

Yi=[1e-3,0]
t=numpy.linspace(0,30,1000)
Y=odeint(systeme,Yi,t,rtol=1e-7,atol=1e-7)
x=Y[:,0]
v=Y[:,1]

figure()
plot(t,x)
xlabel('t')
ylabel('x')
grid()


    			

(freq,sp) = spectre(t,x)
figure()
plot(freq,sp)
xlim(0,5)
xlabel('f')
ylabel('A')
grid()   			
    			

Pour un oscillateur linéaire, l'instabilité se traduit par des oscillations à la fréquence propre (soit la moitié de la fréquence d'excitation). L'amplitude de ces oscillations croît indéfiniment.

Voyons le cas de l'équation du pendule :

w2=(2*numpy.pi)**2
omega = 2*numpy.pi*2
k=2e-2
h=0.18
def systeme(y,t):
    return [y[1],-w2*(1+h*numpy.cos(omega*t))*numpy.sin(y[0])-k*y[1]]

Yi=[1e-3,0]
t=numpy.linspace(0,100,1000)
Y=odeint(systeme,Yi,t,rtol=1e-7,atol=1e-7)
x=Y[:,0]
v=Y[:,1]

figure()
plot(t,x)
xlabel('t')
ylabel('x')
grid()


    			

(freq,sp) = spectre(t,x)
figure()
plot(freq,sp)
xlim(0,5)
xlabel('f')
ylabel('A')
grid()   			
    			

La non-linéarité (provenant du sinus) a pour effet de limiter l'amplitude des oscillations et de faire apparaître des battements de l'amplitude, conformément à ce qui est observé expérimentalement.

Voici, pour la même fréquence d'excitation, le calcul avec une amplitude d'oscillation (paramètre h) plus petite :

w2=(2*numpy.pi)**2
omega = 2*numpy.pi*2
k=2e-2
h=0.06

Yi=[1e-3,0]
t=numpy.linspace(0,100,1000)
Y=odeint(systeme,Yi,t,rtol=1e-7,atol=1e-7)
x=Y[:,0]
v=Y[:,1]

figure()
plot(t,x)
xlabel('t')
ylabel('x')
grid()


    			

La vitesse de croissance des oscillations et celle des battements sont plus faibles, conformément aux observations expérimentales.

Voici le cas d'une excitation à la fréquence propre (n=2) :

w2=(2*numpy.pi)**2
omega = 2*numpy.pi
k=2e-2
h=0.2

Yi=[1e-3,0]
t=numpy.linspace(0,500,10000)
Y=odeint(systeme,Yi,t,rtol=1e-7,atol=1e-7)
x=Y[:,0]
v=Y[:,1]

figure()
plot(t,x)
xlabel('t')
ylabel('x')
grid()


    			

(freq,sp) = spectre(t,x)
figure()
plot(freq,sp)
xlim(0,5)
xlabel('f')
ylabel('A')
grid()   			
    			

Pour une même valeur de h, la croissance des oscillations est beaucoup plus lente et l'amplitude maximale est plus faible. La période des battements est aussi beaucoup plus grande. C'est bien ce qu'on observe expérimentalement.

Cette simulation numérique confirme qualitativement les observations expérimentales mais la croissance de l'amplitude des oscillations est plus rapide dans la réalité (pour une même valeur de h).