スポーツのための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 最大値 min 最小値 sum 合計 std 標準偏差

②行列の結合 行列同士を結合するには、行列 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

黒

k

点線

：

プラス +

青色 b

鎖線

-- 四角

s

赤色 r

一点鎖線 -. 下三角 v

緑色 m

五角形 p

紫色 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 1 2 2×2行列 3×2行列 1 2 1 2 3 4 3 4 5 6 2 左右に並べる 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)

⑪関節角度の算出 関節角度を算出するには、「上部セグメント角度」と「下部セグメント角度」の ２つを求めて、その２つのパラメーターを用いて関節角度を算出する。 ベクトルの作り方 中心となる関節を起点としてベクトルを２つ作る。 終点－起点＝起点を中心としたベクトル 大転子－膝＝大腿部ベクトル 外果 －膝＝下腿部ベクトル （課題１ １列行列ごとに作る 課題２ ３列行列ごとに作る）

角度算出 セグメント角度を求めるには、[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 を読んでみる。 わからない数学の部分は、数学の本を借りたり買ったりして理解する。

⑰ダウンサンプリング ダウンサンプリング：地面反力が 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/○○); フォースプレートのデータにノイズがあるので取 り除くためにフィルターを作ります。フォースプレ ート 1080Hz、 1080÷●●すると、○○に入力する数字がわかりま す。 FF2=filtfilt(B,A,FF2); 位相ずれなしのバターワースフィルターをかける。 ①基本の時間軸を作る 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) １行目をあけたい場合（ヘッダーを作る場合）は、１，０と入力する