Sous-programme de passage du centre de pression au centre de masse
Pour les plates-formes de stabilométrie.

Auteur : Bernard Gagey bernard.gagey@gmail.com

- Langage :

• Octave (compatible Matlab)
  

- Restrictions : Ce programme

• o Correspond aux mesures où l’approximation de l’individu par un « pendule inversé » est valable.
  o Le résultat n’est pas garanti dans les premières et dernières secondes (1à 3s)
  o Le programme contient un paramètre « Correction » qui vaut 1 + le moment d’inertie du corps par rapport au centre de gravité divisé par la masse. Par défaut de connaissance, cette valeur est approximée à 1,2 .
  

- Paramètres d’entrée

• o Vecteur à 1 dimension des coordonnées de la courbe du centre de pression sous la forme complexe x+iy
  o Fréquence d’échantillonnage
  o Taille de l’individu
  

- Sortie :

• Vecteur à 1 dimension des coordonnées de la courbe du centre de masse sous la forme complexe x+iy
  Remarque : Pour éviter les complexes, on peut travailler successivement avec les abscisses x, puis les ordonnées y : les complexes ne sont ici qu’un artifice permettant d’éviter ce doublet.
  

- Explication :

Le corps est supposé se comporter comme un pendule inversé, et satisfaire à l’équation différentielle :

P=M-k²M’’

Toutes les solutions de cette équation sont données, à partir d’une première solution M_ø par la formule :

M=M_ø +(ae^kt+be^-kt)

On recherche d’abord la solution M_ø qui serait en Ø au top précédent et suivant la mesure.
Pour cela on écrit l’équation différentielle en chaque point :

P(j)=M(j)-k²(M(j-1)+M(j+1)-2M(j))/dt²

Ayant supposé que M(Ø) et M(N+1)=Ø, ceci constitue un système de N équation linéaire à N inconnues qu’il suffit de résoudre, ou autrement dit :

Vecteur P= matrice A.Vecteur M
Vecteur M= (Matrice A)^-1.Vecteur P

La première étape du programme consiste à construire la matrice
La deuxième étape consiste a résoudre le système linéaire directement ou par inversion de la matrice
La 3ème étape consiste à trouver de bons coefficients a et b : Ils sont mathématiquement impossibles à trouver, mais physiquement les points M(Ø) et M(N+1) doivent au moins être sur la plateforme . La formule utilisé dans le programme a été trouvé empiriquement, mais on sait qu’elle n‘est pas, de loin, parfaite. En commentaire est mis une 2ème formule plus facile à comprendre. De toutes manières, au bout de 1 à 2s, l’importance de cette correction est négligeable.

% Début du sous-programme

function[g]=SpG(posi,freqech,taille);
%Programme de recherche de la courbe G
% parametres d'entrée:
% posi la courbe en coordonnées imaginaires
% freqch la frequence d'echantillonage (de la courbe totale) =nombre d'échantillons par seconde
% taille taille de l'individu corrigée si on ne prend pas le centre de gravité à 0.55
%retour
% g la dernière estimation de la courbe g qui satisfait à posil=g -kg''
%
gravite=9.8066;
H_gravite=.55;%_____________rapport hauteur centre de gravité / hauteur totale
Correction=1.2;% coefficient correctif de k Minimum 1 pour masse en G , 1,37 pour masse aux extrèmes 1.33 pour barre
ak=H_gravite*Correction/gravite;
coeffaccel=ak*taille
kacc=(1/coeffaccel)^.5;
Nbm=length(posi);
%centrer pour que le point (0,0) ne soit pas absurde
meanposi=mean(posi);
posi=posi-meanposi;
%calcul des coefficients de la matrice
co1=-freqech*freqech*coeffaccel;co2=-2*co1+1;
%creation de la matrice avec sa bonne diagonale
ma=co2*eye(Nbm,Nbm);
%remplissage des sus et sur-diagonale
for k=2:1:Nbm-1;ma(k,k-1)=co1;ma(k,k+1)=co1;endfor
ma(1,2)=co1;ma(Nbm,Nbm-1)=co1;
% Première méthode par inversion de matrice semble-t-il plus rapide sous Octave à cause de la forme de la matrice
na=inv(ma);
g=na*posi;
% 2ème méthode (en commentaire)par résolution de système linéaire qui devrait être plus rapide
%g=ma\posi;
%recherche debut et fin
% Première itération qui donnerait l'avant premeir point et le premier identique...
correc=zeros(Nbm,1);for jj=1:1:Nbm;correc(jj)=e^(-kacc*jj/freqech);endfor;
deb=g(1)/(1-correc(1));
fin=g(Nbm)/(1-correc(1));
%2ème itération
mul1=76;mul2=75;
deb=mul1*g(1)-mul2*g(2)+deb*(mul1*correc(1)-mul2*correc(2));
fin=mul1*g(Nbm)-mul2*g(Nbm-1)+fin*(mul1*correc(1)-mul2*correc(2));
% autre possibilité considérer que les points extrêmes de Vg sont ceux de Vp
%deb=(Vp(1)-Vg(1))/correc(1);fin=(Vp(Nbm)-Vg(Nbm))/correc(1);
%Effectuer l'addition de la correction et oter le centrage
g=meanposi+(g+deb*correc+flipud(fin*correc));
endfunction