スポーツのための Matlab 講習会

スポーツのための Matlab 講習会
対象者：VICON 利用を考えている学部生・大学院生
卒論、修士論文、博士論文のバイオメカニクスデータを分析するにあたって、
多変量、多変数のパラメーターを分析するとなると、エクセルを用いて分析を
していては、何ヶ月もかかって効率が悪い。そこで、行列計算を行うことがで
きる Matlab というプログラム言語を利用することで、短時間に分析することが
でき、研究の他の部分に時間をまわすことができる。
Matlab に関する本は多くあるが、スポーツの分析に必要な方法がバラバラに
記述されており、独学すると習得するために莫大な時間かかる。そこで、効率
よく Matlab の使い方を講義して、実践的に使えるようにするのが本講習会の目
的である。
①行列の抽出
BMI＝体重（Kg）÷（身長（ｍ）×身長（ｍ））
RJ_index=跳躍高（ｍ）/接地時間（秒）
跳躍高(滞空時間から求める)＝（1/8）×9.81（重力加速度 m/s）×滞空時間２
A=[1 2 3；4 5 6；7 8 9]
A=1 2 3
4 5 6
7 8 9
縦が「行」、横が「列」 行と列を区別できるようにする
A=[1 2 3；;4 5 6；7 8 9]
A=[1 2 3；4 5 6；7 8 9]； 式の最後に；セミコロンをつけると画面表示され
ない
2 列目だけ抽出する QQ(：，２)
4 行目だけ抽出する QQ（４，：）
：コロンをつけると、すべての行、すべての列を抽出する。
2 列目～5 列目を抽出する QQ(：，２：７)
抽出したい変数（はじめの列：終わりの列）
4 行目～7 列目を抽出する QQ（４：７，：）
抽出したい変数（はじめの行：終わりの行）
mean 平均値 max 最大値
sum 合計 std 標準偏差
median 中央値 var 分散
min
最小値
②行列の結合
行列同士を結合するには、行列 A と行列 B の行の大きさ、列の大きさが同じで
ないと結合できない。
行列
行列
A＝[1 2 3 4 5]
B＝[5 6 7 8 9]
行列を横に結合するには、変数間にスペースを空けて並べてから、括弧で囲む
[
]。
C=[A B]
C=[1 2 3 4 5 5 6 7 8 9]
行列を縦に並べるには、並べたい変数の間にセミコロン（；）を入れて、括弧
で囲む[ ]。
C=[A；B]
D=１：２：１０
E=２：３：３３
問題
D と E を横に並べてください
E と D を縦に並べてください。
行行列を列行列に変更する
列行列を行行列に変更する。
D’
アポストロフィを付ける（Shift を押して 7 のキーを押す）
③エクセルデータの取り込み
[DATA HEAD]=xlsread('Rotation.xls');
「xlsread」はエクセルのファイルを読み込む関数です。
エクセルファイルを読み込ますには、
「エクセル 5.0/95 ブック」の形式にする。
エクセルの最新のファイル形式は読み込まないので注意。
ファイル名は ’ ’ で囲む。
④グラフの作成
figure(1) 複数のグラフを作成するときには必ず先頭に入力する。
グラフの表示
plot(Time,DATA(:,2),Time,DATA(:,4),Time,DATA(:,6))
plot(X,Y,,X2,Y2,X3,Y3)
X 軸データ、Y 軸データの順に入力していく
タイトル、X 軸、Y 軸の表示方法
title('Rotation')
xlabel('Time(sec)')
ylabel('Angle(deg)')
‘ ’で表示させたい文字を囲む
凡例の表示
legend('Lower','Upper','Hip')
plot （X1Y1
X2Y2
X3Y3）に対応させて入力する。
‘ ’で表示させたい文字を囲む
凡例の位置の設定
legend('Lower','Upper','Hip',-1)
-1
１
２
３
４
図表の外に凡例が表示される。
図表の右上に表示される。
図表の左上に表示される。
図表の左下に表示される。
図表の右下に表示される。
⑤グラフの線、プロットのサイズの設定
plot(Time,DATA(:,2) ，’-oｋ ’，Time,DATA(:,4)，’-vｋ ’，Time,DATA(:,6) ，’
-sｋ ’，)
X軸データ、Y軸データ，
’グラフの線の種類、マーカの種類、線・マーカーの色の設定’,
X と Y のデータのあとに入力する。
線種類
直線
点線
：
鎖線
-一点鎖線 -.
マーカの種類
丸
プラス
四角
下三角
五角形
o
+
s
v
p
色
黒
青色
赤色
緑色
紫色
k
b
r
m
m
X 軸、Y 軸の大きさを調整する。
axis([X 軸最小値 X 軸最大値 Y 軸最小値 Y 軸最大値])
特に設定しないときは、最小値には-Inf
最大値には Inf を入力する。
⑥行列の行の大きさをもとめる
LL=length(DATA)
「length」は行・列で最も長い「行」「列」を調べる関数です。
⑦行列の行と列の大きさを求める
size(DATA)
size は行の長さ、列の長さを同時に調べる関数です。
⑧ヘッダーの検索、データーの抽出
エクセルデータからヘッダーの位置を検索する。
II=strmatch('LASI:X',HEAD)
「strmatch」を使って、何列目に必要なデータがあるか検索する関数
カッコの中は（‘検索したい文字列’，検索したい文字がある行列）の順に入力
する。
5 列目に必要なデータがあれば、５と表示される。
5 列目に必要なデータがあることがわかったので、5 列目のデータを、すべての
行を抽出するように式を書く。
LASI=DATA(:,II)
VICON のデータの場合は、すべてのデータが X 軸、Y 軸、Z 軸の順に記述されて
いるので、「○○：X」の列番号がわかっていれば、
「○○：X」＋１した列番号が「○○：Y」のデータになり、
「○○：X」＋２した列番号が「○○：Z」のデータになる。
3 列まとめて抽出すると、
LASI=DATA(:,II:II+2)
のように書く必要がある。
VICON は座標データが「ミリ単位」で出力されるので、
メートル単位に単位変換をする必要がある。
II=strmatch('LASI:X',HEAD);
LASI=DATA2(:,II:II+2)/1000;
1000 で割ると、ミリ単位からメートル単位に変換される。
⑨時間軸の作成
II=strmatch('Frame',HEAD);
Frame=DATA2(:,II);
LL=lenght(Frame)
length 関数を用いて、行列の長さを調べる(行の長さを調べる)
Time=0:(1/120):(LL-1)＊(1/120);
120Hz の時間軸作成
（LL-1）をしないと、実際の行数より多くなってしまうのでマイナス１する。
Time2=0:(1/250):(LL-1)*(1/250);
250Hz の時間軸作成
⑩subplot 関数の使い方
１つのグラフに複数のグラフに出力する
上下に並べる
figure(1)
subplot(2,1,1)
plot(x,y,x1,y1)
subplot(2,1,2)
plot(x,y,x1,y1)
subplot(行、列、番号)を用いて、1 つの図表に複数の図表を表示する。
何個の図表を表示させたいか決めて、行と列の大きさを入力する。
２×１行列であれば
上部に表示させたければ１、下部の表示させたければ 2 と入力する。
2×１　行列
1
2
１×２　行列
1
2
2×2行列
1
3
3×2行列
1
3
5
2
4
6
左右に並べる
figure(2)
subplot(1,2,1)
plot(x,y,x1,y1)
subplot(1,2,2)
plot(x,y,x1,y1)
４分割する
figure(3)
subplot(2,2,1)
plot(x,y,x1,y1)
subplot(2,2,2)
plot(x,y,x1,y1)
subplot(2,2,3)
plot(x,y,x1,y1)
subplot(2,2,4)
plot(x,y,x1,y1)
2
4
⑪関節角度の算出
関節角度を算出するには、
「上部セグメント角度」と「下部セグメント角度」の
２つを求めて、その２つのパラメーターを用いて関節角度を算出する。
ベクトルの作り方
中心となる関節を起点としてベクトルを２つ作る。
終点－起点＝起点を中心としたベクトル
大転子－膝＝大腿部ベクトル
外果 －膝＝下腿部ベクトル
（課題１ １列行列ごとに作る
課題２
３列行列ごとに作る）
角度算出
セグメント角度を求めるには、[atan2]という関数を用いて角度を算出する。
角度(ラジアン単位)＝atan2(縦軸，横軸)；
縦軸、横軸を間違えると算出できないので注意。
角度単位変換
ラジアン単位だとわかりづらいので、Radian から Degree へ、角度単位を変換す
る。
角度=（算出した角度＊１８０）＊pi(円周率＝3.14)
%膝関節角度
UPLEG_V=RASI-RKNE;
KKUP=atan2(UPLEG_V(:,3),UPLEG_V(:,2));
UPLEG_ag=(KKUP*180)/pi;
DOLEG_V=RANK-RKNE;
KKDO=atan2(DOLEG_V(:,3),DOLEG_V(:,2));
DOLEG_ag=(KKDO*180)/pi;
RKNE_ag=UPLEG_ag-DOLEG_ag;
⑫フィルタリング
フィルタリングとは，角度、角速度算出にあたって生じるノイズを取り除くた
めに行われる。研究室によって使い方の考え方が異なるが、基本的には角度デ
ータにフィルターをかけたら、２回フィルターをかけない。論文では位相ずれ
なしの４次のバターワースローパスフィルターを、遮断周波数○○Hz で行った
という記述でフィルタリングについて書かれてある。
まず、バターワースフィルターの作成(設計)
[B,A]=butter(フィルターの次数,遮断周波数 / [取り込んだデータの周波数
/2]);
AG=filtfilt(B,A,[フィルターをかけたいデータ]);
filtfilt を使用するのであれば、4 次のフィルターであれば 2 と設定する。
[B,A]=butter(2,10/60);
AG=filtfilt(B,A,RKNE_ag);
⑬速度を算出する
位置座標、角度データ→フィルター→速度、角速度計算
3 点微分法（スポーツバイオメカニクス 阿江参照）
速度の計算方法①
COMassV1=(-3*COMassm(1,:)+4*COMassm(2,:)-COMassm(3,:))/(1/60);
COMassV2=(COMassm(3:LL,:)-COMassm(1:LL-2,:))/(1/60);
COMassV3=(COMassm(LL-2,:)-4*COMassm(LL-1,:)+3*COMassm(LL,:))/(1/60);
COMassV=[COMassV1;COMassV2;COMassV3];
速度の計算方法②
LL=length(x);
AGV=[(-3*x(1,:)+4*x(2,:)-x(3,:))/(1/60);...
(x(3:LL,:)-x(1:LL-2,:))/(1/60);...
(x(LL-2,:)-4*x(LL-1,:)+3*x(LL,:))/(1/60)];
２５０Hz の分析は
LL=length(x);
AGV=[(-3*x(1,:)+4*x(2,:)-x(3,:))/(1/125);...
(x(3:LL,:)-x(1:LL-2,:))/(1/125);...
(x(LL-2,:)-4*x(LL-1,:)+3*x(LL,:))/(1/125)];
⑭時間軸作成 応用編
膝の最大屈曲時を 0 秒にして、時間軸を作成する。
[RKneMin RKneMinF]=min(RKNE_AG);
Time=(0:(1/120):(1/120)*(RR-1))';
Frame=(1:RR)';
figure(7)
subplot(2,1,1)
plot(Time,RKNE_AG)
xlabel('Time(sec)')
subplot(2,1,2)
plot(Frame,RKNE_AG)
ylabel('Frame')
axis([1 Inf -Inf Inf])
RKneMinF
RR
Excel ファイル「時間軸」を参照
%わからなくなったら Excel で確認しながら作成する。
TimeA=0:(1/120):(1/120)*(RR-RKneMinF);
TimeB=-(1/120)*(RKneMinF-1):(1/120):-(1/120);
%0 は含まないので、行列の数をマイナス 1 する必要がある。
Time2=[TimeB TimeA]';
%2 つの行列を合成する。
⑮Matlab 情報の検索
Matlab プログラムをしてわからなくなったら？
１．ヘルプを見る。
２．Google で検索してみる。「Matlab ○○(捜したいキーワード)」
Yahoo は Matlab 関連の検索はあまりよくないので使っていない。
３．http://www.cybernet.co.jp/matlab/
サイバーネットシステムのサイトに行く
４．テクニカル FAQ で検索する。
http://www.cybernet.co.jp/matlab/support/techkwdb/
５．Mathwork 社（Matlab の開発をしている）のページに行く。
http://www.mathworks.com/matlabcentral/
File Exchange に行くと、人が作ったファイルを手に入れることができる。
６．Matlab の本を買って読んでみる。
マニア向け（BASIC や C 言語などの他のプログラムの発想を取り入れる）
⑯Matlab の勉強法
予備実験でもいいので実験を行って分析してみる。
講習会の資料を基にして、データーを分析してみる。
Matlab の本を一冊はじめから、終わりまで、式を書きながらやってみる。
David A.Winter. Biomechanics John Wiley & Sons,Inc を読んでみる。
わからない数学の部分は、数学の本を借りたり買ったりして理解する。
卒論のプログラムを作ってみる（はじめは 1 週間以上かかる）
⑰ダウンサンプリング
ダウンサンプリング：地面反力が 1080Hz で取り込まれていて、VICON のデー
タが 120Hz で取り込まれている場合、地面反力のデータを 9 分の 1 にして、VICON
のデータにあわせること。
サンプリングが小さいデータに、サンプリングが大きいデータを、データ数
を少なくしてあわせること。
規格化：0.13 秒でリバウンドジャンプする人と、0.22 秒でリバウンドジャンプ
をする人の各関節の動きを比較したいときに、時間の長さを同じにして比較す
る方法。
ダウンサンプリングの練習
clc close all clear all は入力しておく。
%データ取り込み
[DATA HEAD]=xlsread('○○'); エクセルのファイルを開く
FDATA=csvread('yano12.csv',10,0);
CSV ファイル（カンマ区切り）ファイ
ルを開く
10,0 は、11 行目、1 列目からデータを抽出することを示している
%ダウンサンプリング
FF2=FDATA(:,○○:○○);
る。
QQ=length(FF2);
フォースプレート２のデータが必要なので抽出す
FDATA をプロットしてください。
タイトル、Ｘ軸、Ｙ軸、凡例をつけてください。
axis([Xmin Xmax Ymin Ymax])
[B,A]=butter(4,40/○○);
FF2=filtfilt(B,A,FF2);
フォースプレートのデータにノイズがあるので取
り除くためにフィルターを作ります。フォースプレ
ート 1080Hz、
1080÷●●すると、○○に入力する数字がわかりま
す。
位相ずれなしのバターワースフィルターをかける。
①基本の時間軸を作る
Time=(0:(1/1080):(QQ-1)*(1/1080))'; サンプリングレート 1080Hz のデータ
Ｔｉｍｅ＝0：0.001：２
②変更したい時間軸を作る
TimeS=(0:(1/120):(QQ-1)*(1/1080))'; サンプリングレート 120Hz のデータ
ＴｉｍｅＳ＝0：0.008：2
③ダウンサンプリングする
出力=spline(①,ダウンサンプリングしたいデータを入力,②);
※1 行ごとしかダウンサンプリングできないので注意。
FX2=spline(Time,FF2(:,1),TimeS);
ＦＸ ＦＹ ＦＺごとを合成して出力する。
⑱規格化
MM=length(RANKLE_ExFx);
TimeA=(0:(1/120):(1/120)*(MM-1));
基本の時間軸
MaxA=max(TimeA); 時間の最後を捜してくる。
InterA=MaxA/100;
スタートから時間の最後までを 100 等分する。
※一般的には 100％に規格化することが多い。
ATime=(0:InterA:MaxA)';
新しい時間軸
=spline(TimeA【規格化前の時間軸】,規格化したいデータ,ATime【規格化後の
データ】);
ANKLE_SP=spline(TimeA,RANKLE_ExFx,ATime);
規格化（実際の場合）
①分析したい区間の時間軸を新しく作る。
②分析したい区間を 100 等分する。【上記の方法と同様】
③spline(①元の時間軸，データ（分析したい区間），②新しい時間軸)
出題：右膝の最大屈曲時から最後のフレームまでを規格化する。
[RKNE_max RKNE_maxF]=max(RKNE_ExFx); ①右膝の最大屈曲時のフレームを捜
す
BB=length(RKNE_ExFx); ②最終フレームを捜す
BL=BB-RKNE_maxF; ①スタートから②終了フレームまでに、何フレームある求
める。
TimeBB=(0:(1/120):BL*(1/120))'; 新しい時間軸を作る
新しい時間軸を 100 等分する
MaxB=max(TimeBB);
InterB=MaxB/100;
TimeCC=(0:InterB:MaxB)';
３次のスプライン関数を用いて規格化する。
ANKLE_SP=spline(TimeBB,RKNE_ExFx(RKNE_maxF:BB),TimeCC);
※規格化したい区間だけを計算にかける必要がある。
PTime=0:100; 規格化用の時間
直線補間
地面反力を補間するときはこちらの方がよい場合もある
SP_EX=interp1(Time,DATA,STime,'linear');
⑲グラフを出力する。
saveas(figure(1),'yano12','tif')
保存したいグラフが figure(2)であれば、figure(2)と入力する。
⑳データを出力する。
計算したデータを、エクセルファイルに出力する
PADATA=[A;B;C;D];
出力したいデータは［］カッコで囲む。
PAHead={'Data_A','DATA_B','DATA_C','DATA_D'};
ヘッダーをつけるときは、｛｝カッコで囲む。文字には’’で囲む。
xlswrite('yano12_KickX.xls',PADATA,'sheet1','B1');
xlswrite('yano12_KickX.xls',PAHead,'sheet1','A1');
※行列の形がそろっているか注意
計算したーデータを、ロータス 123 ファイルに出力する
wk1write('Matlab1',PADATA,1,0)
１行目をあけたい場合（ヘッダーを作る場合）は、１，０と入力する