Boussole dans un champ magnétique

from matplotlib.pyplot import *
import numpy
import math
import cmath
import numpy.fft
import scipy.signal
import StringIO
            

3.a. Calcul de l'angle

On commence par lire les données dans le fichier texte. Les virgules doivent être converties en points. On saute les deux premières lignes du fichier. Pour chaque colonne, on enlève les 5 premières valeurs, qui correspondent à un marqueur situé au centre de la boussole.

s = open('CIMG5043.txt').read().replace(',','.')
data = numpy.loadtxt(StringIO.StringIO(s),skiprows=2,unpack=True)
debut = range(5)
t=numpy.delete(data[0],debut)
x=numpy.delete(data[1],debut)
y=numpy.delete(data[2],debut)
figure(figsize=(6,6))
plot(x,y)
xlabel("x")
ylabel("y")
grid()
                

figA.pdf

Pour obtenir l'angle, le plus simple est de créer un tableau avec des nombres complexes et d'en extraire l'argument. Comme la position d'équilibre de l'anneau se trouve à gauche de l'image, on change le signe de x.

angle = numpy.angle(-x+1j*y)

figure()
plot(t,angle*180.0/math.pi)
xlabel("t (ms)")
ylabel("angle (deg)")
grid()
                 

figB.pdf

Voici un détail vers la fin du mouvement :

figure()
plot(t,angle*180.0/math.pi)
xlabel("t (ms)")
ylabel("angle (deg)")
grid()
axis([6e4,7e4,-10,10])
                 

figB1.pdf

On voit sur ce tracé la quantification de l'angle, qui vient de celui des coordonnées (x,y) qui varient par saut de 1 pixel. Par rapport à la courbe idéale qui devrait être obtenue en signal analogique, la numérisation a introduit un bruit de quantification.

3.b. Analyse spectrale

On commence par une analyse spectrale globale, sur toute la durée de la trajectoire.

fe=120.0
te=1.0/fe
N=y.size
T=1.0/fe*N
f = numpy.arange(angle.size)*1.0/T
spectre = numpy.abs(numpy.fft.fft(angle))

figure()
plot(f,spectre,'-o')
xlabel('f (Hz)')
ylabel('A')
grid()
axis([0,2,0,spectre.max()])
                

figC.pdf

Le spectre est assez complexe car d'une part le mouvement est pseudo-périodique (périodique amorti), d'autre part la pseudo-période évolue au cours du mouvement (elle diminue lorsque l'amplitude diminue).

Pour y voir plus clair, on fait deux analyses spectrales, chacune sur une moitié de la trajectoire.

N1 = int(N/2)
spectre1 = numpy.abs(numpy.fft.fft(angle[0:N1]))
f1 = numpy.arange(spectre1.size)*1.0/(N1*te)
spectre2 = numpy.abs(numpy.fft.fft(angle[N1:N]))
f2 = numpy.arange(spectre2.size)*1.0/((N-N1)*te)

figure()
plot(f1,spectre1,'-o',label='[0,T/2]')
plot(f2,spectre2,'-o',label='[T/2,T]')
xlabel('f (Hz)')
ylabel('A')
grid()
legend(loc='upper right')
axis([0,2,0,spectre.max()])
                 

figD.pdf

Les spectres sont deux fois moins précis puisque la durée a été divisée par deux. Sur le spectre de la seconde moitié de la trajectoire, on voit le maximum correspondant aux oscillations de petite amplitude, pour un angle inférieur à 30 degrés environ.

On en déduit la fréquence des oscillations de petite amplitude :

$$f_0=0,475\pm 0,025Hz\text{(2)}$$

D'après l'équation du mouvement, cette fréquence est :

$$f_0=\frac{1}{2\pi }\sqrt{\frac{mB}{J}}\text{(3)}$$

Au lieu où l'expérience est réalisée, le champ magnétique terrestre peut être estimé (Earth magnetic field) à

$$B=47,5\pm 0,5\mu T\text{(4)}$$

On en déduit :

$$\frac{m}{J}=(1,87\pm 0,14)\cdot 10^5{s}^2\cdot T^{-1}\text{(5)}$$

On peut alors utiliser cette valeur pour déterminer l'intensité d'un champ magnétique (par exemple celui d'un aimant), en faisant osciller cette boussole dans le champ.

3.c. Vitesse angulaire

La vitesse angulaire est obtenue avec un filtre dérivateur, c'est-à-dire un filtre RIF dont la réponse impulsionnelle est (1,-1) :

a=[te]
b=[1,-1]
omega = scipy.signal.lfilter(b,a,angle)
omega[0] = omega[1]

figure()
plot(t,omega)
xlabel("t (ms)")
ylabel("w (rad/s)")
grid()
                

figE.pdf

Voyons un détail :

figure()
plot(t,omega)
xlabel("t (ms)")
ylabel("w (rad/s)")
axis([6e4,6.2e4,-2,2])
grid()
                

figF.pdf

La vitesse n'a que trois valeurs possibles dans cet intervalle. Cela vient de la quantification de l'angle, qui varie par sauts. On voit donc que le bruit de quantification est très amplifié par la dérivation. D'une manière générale, la dérivation amplifie le bruit. Le signal obtenu ici est pratiquement inutilisable.

La solution est de lisser le signal avant d'effectuer la dérivation. Pour cela, on peut utiliser un filtre RIF passe-bas. Dans le cas présent, la fréquence d'échantillonnage (120 Hz) est très grande devant la fréquence du signal (environ 0,5 Hz). On peut donc utiliser un filtre passe-bas avec un indice de troncature assez élevé, ce qui donnera un très bon filtrage. Voici le résultat avec une fréquence de coupure de 5 Hz. On superpose le signal filtré avec le signal d'origine pour vérifier que l'amplitude n'est pas modifiée :

P=30
fc=5.0/120
h = scipy.signal.firwin(numtaps=2*P+1,cutoff=[fc],nyq=0.5,window='hamming')
angle_filtre = scipy.signal.convolve(angle,h,mode='same')
figure()
plot(t,angle*180.0/math.pi,'b')
plot(t,angle_filtre*180.0/math.pi,'r')
xlabel("t (ms)")
ylabel("angle (deg)")
grid()
                      

figG.pdf

Voici le même détail que plus haut :

figure()
plot(t,angle*180.0/math.pi,'b')
plot(t,angle_filtre*180.0/math.pi,'r')
xlabel("t (ms)")
ylabel("angle (deg)")
grid()
axis([6e4,7e4,-10,10])
                      

figH.pdf

Le lissage a bien fait disparaitre le bruit de quantification. On peut à présent effectuer la dérivation :

a=[te]
b=[1,-1]
omega_filtre = scipy.signal.lfilter(b,a,angle_filtre)

figure()
plot(t,omega_filtre)
xlabel("t (ms)")
ylabel("w (rad/s)")
axis([0,8e4,-8,8])
grid()
                      

figI.pdf

3.d. Portrait de phase

Pour finir, on trace la trajectoire dans le plan de phase.

figure()
debut=P
plot(angle_filtre[debut:N],omega_filtre[debut:N])
xlabel("angle (rad)")
ylabel("omega (rad/s)")
grid()
                

figJ.pdf