Matlab 講習会
目的 Matlab を用いて、ＶＩＣＯＮや Winanalyze の座標データー、 地面反力の分析必要な項目について習得する。 本やヘルプに掲載されている情報を、実際に使用できる形で整理する。 講習会１回目（４時間） １．行列操作について理解する ２．時間軸を作る ３．エクセルデーターを取り込む ４．テキストデーターを取り込む ５．グラフの作成 ６．１つのグラフに複数のグラフを出す ７．計算によって新しい座標軸を作る 講習会２回目（６時間） ８．フィルターをかける ９．速度を算出する １０．合成速度を算出する １１．最大値、最小値を求める １２．時系列データをまとめる １３．データを出力する １４．図表を出力する １５．スティックピクチャーを書く＆軌跡を描く １６．セグメントを書く＆axis equal 縦横比をそろえる １７．画面に出力するスティックピクチャーの数を減らす １８．関節角度を算出する １９．補間作業 ２０．地面反力のサンプリングを下げる１． 行列操作について理解する a1=[1 2 3] 行行列 間にスペースを入れる a2=[1;2;3] 列行列 間にセミコロンを入れる b1=[1 2 3 4 5]；セミコロンを入れると非表示 A=[1 2 3;1 2 3;1 2 3] 1 2 3 1 2 3 1 2 3 →行 ↓列 １行目だけ取り出す A(1,:) １列目だけ取り出す A(:,1) 問題 Q1 ３行目だけ取り出す Q2 ３列目だけ取り出す ２～３列目を取り出す Ａ（:,2:3） 2～３行列を取り出す A(2:3,2:3)
２． 時間軸を作る 時間軸を作る １秒間 120 コマであったら 1/120=0.0083 Time=0:0.0083:1 101 行のデータに,0.0083 刻みの時間軸を作るには DA=0:100 このままだと１０１の列行列になる DA=(0:100)’ A.行の長さを調べる→ B.最終行までの時間を計算する→ C.時間軸を作る A LL=length(DA) B EndTime=LL*0.0083 C TIME=0:0.0083:EndTime 行が長くなってしまうので TIME=0:0.0083:EndTime-0.0083 のように書く 別解 EndTime=（(LL-1）*0.004） TIME=0:0.0083:EndTime べき乗 ２．＾３ ２の３乗 power(2,3) ２の３乗 ゼロの行列がほしい場合 ZZZ=zeros(10,1) １の行列がほしい場合 OOO=Ones(1,10)
３．エクセルデーターを取り込む clear all [DATA HEAD]=xlsread('Matlab1.xls'); と入力する DATA には、数値データが入る HEAD には、ヘッダーが入る（重複ヘッダーがないか注意） 取り込んだデータを、マーカーごとに変数にする II=strmatch('LASI:X’,HEAD) 探したい文字は、’ ’ で囲むようにする HEAＤ（ヘッダー）のなかに、’LASI:X’が何行目にあるか探してくる（４行目） LASI=DATA(:,II) （Field=DATA(:,4)と入力しているのと同じ） 行はすべてほしいので、：と入力する X Y Z のすべてほしい場合は LASI=DATA(:,ll:II+2) (Field=DATA(:,4:6)と入力しているのと同じ) VICON を用いているときは、３行ごと取り込むようにしている。慣れない場合は１行ごと 取り込む（計算によっては１行ごと取り込んだ方が計算しやすいことも） 位置座標の単位を、ミリからメーターに変える（線速度などを出力するとき必要） LASI =DATA(:,ll:II+2)/1000
４．テキストデーターを取り込む LANK=dlmread('Read1-1.txt','¥t',1,0); ヘッダーの部分は取り込めないので、「１」と入力して２行目以降を取り込む KNEE=dlmread('Read1-2.txt','¥t',1,1); と入力すると、２列目以降が取り込まれる LASI=dlmread('Read1-3.txt','¥t',1,0); ５．グラフの作成 位置座標をプロットする 時間軸を作る LL=length(KNEE) TimeA=0:0.004:0.004*(LL-1) figure(1) plot(TimeA,KNEE(:,1),TimeA,KNEE(:,2),TimeA,KNEE(:,3)) タイトルを入れる、Ｘ軸タイトル、Ｙ軸タイトル title(‘Knee’) ‘ ’で文字を挟む xlabel(‘Time(sec)’) ylabel(‘mm’) 凡例（はんれい） legend(‘X’,’Y’,’Z’) plot の順番と対応させて書く
表示の範囲を設定する
axis([0 0.15 -Inf Inf]) （ ）の中に、[ ] を入力する
はじめの２つが、Ｘの最小値、最大値、残りの２つがＹの最小値、最大値 軸の範囲を設定しない場合は、最小値には-Inf、最大値には Inf と入力 プロットの線、マークの種類を変更する 線 点（プロットシンボル） 色の順で配列する k 黒、b 青、r 赤、g 緑、c シアン、m マゼンダ X バツ Ｏ まる ＋プラス ＊アスタリクス s 正方形 d ダイアモンド v 三角形 - 実線、:点線、;鎖線、--点線 マーカーサイズの設定方法 plot(x,y,’-ok’,‘MarKerSize’,2) MarkerSize 数字を変えることによって、マーカーのサイズが変化する plot(Time,LASI,'ok',Time,LASIf,'-k','MarkerSize',2)
６．１つのグラフに複数のグラフを出す 上下に並べる figure(1) subplot(2,1,1) plot(x,y,x1,y1) subplot(2,1,2) plot(x,y,x1,y1) 左右に並べる 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)
７．計算によって新しい座標軸を作る 両肩峰を結んだ点を作る (A-B)/2+B CSHO=((RSHO-LSHO)/2)+LSHO; ８．フィルターをかける フィルターの作成 [b,a]=butter(２,10/60); ４次のバターワースフィルター（ローパスフィルター） 遮断周波数／（ナイキスト周波数分析に用いた周波数／２） フィルターをかける [b,a]=butter(２,10/60); LASIf=filtfilt(b,a,LASI);
９．速度を算出する 位置座標、角度データ→フィルター→速度、角速度計算 COMassM=COMass/1000; ％単位変換 メートルからセンチに変える COMassm=filtfilt(b,a,COMassM); ％フィルター 速度の計算方法 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]; １０．合成速度を算出する CGV= sqrt(power(COMassV(:,3),2)+power(COMassV(:,2),2),power(COMassV(:,1),2));
１１．最大値、最小値を求める [CGVmax CGVmaxTime]=max(CGV) CGVmax 最大値 CGVmaxTime 最大値出現時のフレーム [CGVmax CGVmaxTime]=min(CGV) 最小値を求める １２．時系列データをまとめる
EX=[COMass LANK LASI Force Field]
１３．データを出力する wk1write('Matlab1',EX,1,0) １行目をあけたい場合（ヘッダーを作る場合）は、１，０と入力する １４．図表を出力する saveas(figure(1),'Matlab1','tif') bmp Windows ビットマップ jpg JPEG イメージ tif TIFF イメージ
１５．スティックピクチャーを書く 図を重ねあわせるコマンド（hold on） figure(1) plot(Time,LKNE(:,3)) hold on plot(Time,LKNE(:,2)) 軌跡を描く 膝の軌跡をプロットする figure(1) hold on for I=1:LL plot(LKNE(I,2),LKNE(I,3),'ok') end 膝と骨盤の軌跡をプロットする figure(1) hold on for I=1:LL plot(LKNE(I,2),LKNE(I,3),'ok') plot(LP(I,2),LP(I,3),'hk') end figure(1) hold on for I=1:LL plot(LKNE(I,2),LKNE(I,3),'ok') plot(LP(I,2),LP(I,3),'hk') plot(LANK(I,2),LANK(I,3),'vk') end
１６．セグメントを書く
結線したい A とＢを、Ｘ座標では[(AX),(BX)]、Ｙ座標では[(AY),(BY)]と入力する。
plot( X , Y ) plot( [ (AX),(BX) ] , [ (AY),(BY) ] Ｘ座標 Ｙ座標 下腿のセグメントを作る for I=1:LL plot([LKNE(I,2),LANK(I,2)],[LKNE(I,3),LANK(I,3)],'k') end 大腿のセグメントを作る for I=1:LL plot([LKNE(I,2),LANK(I,2)],[LKNE(I,3),LANK(I,3)],'k') plot([LP(I,2),LKNE(I,2)],[LP(I,3),LKNE(I,3)],'k') end axis equal 縦横比をそろえる ‘--K’ ‘:K’ 結線の線の種類を変更することができる。 問題 Q1 つま先をプロットする Q2 外果からつま先を結線してみる for I=1:LL plot(LTOE(I,2),LTOE(I,3),'*k') plot([LANK(I,2),LTOE(I,2)],[LANK(I,3),LTOE(I,3)],'k') end
１７．画面に出力するスティックピクチャーの数を減らす for 文を変更する StickFrame=3 for I=1:StickFrame:LL figure(2) hold on StickFrame=2 for I=1:StickFrame:LL plot(LKNE(I,2),LKNE(I,3),'ok') plot(LANK(I,2),LANK(I,3),'vk') plot(LP(I,2),LP(I,3),'hk') plot(LTOE(I,2),LTOE(I,3),'*k') plot([LKNE(I,2),LANK(I,2)],[LKNE(I,3),LANK(I,3)],'k') plot([LP(I,2),LKNE(I,2)],[LP(I,3),LKNE(I,3)],'--k') plot([LANK(I,2),LTOE(I,2)],[LANK(I,3),LTOE(I,3)],'k') end axis equal
１８．関節角度を算出する 関節角度は、２つのセグメント角度を求め、その２つの角度を加算減算することによっ て求める。 膝であれば、地面に対する大腿部セグメント角度と下腿部セグメント角度を求め、 その２つを加算して、３６０度から引くと関節角度が算出される。 ①動くマーカーの座標を求める （原点に対して動くマーカー）－（原点になるマーカー） ②動くマーカーの座標を、関数 atan2（アークタンジェンド）に入力する atan(Y 座標 , X 座標)
UPPLeg=LP-LKNE; LOWLeg=LANK-LKNE; kk1=atan2(UPPLeg(:,3),UPPLeg(:,2)); LKNEag1=(kk1*180)/pi; %ラジアン→角度 変換 kk2=atan2(LOWLeg(:,3),LOWLeg(:,2)); LKNEag2=(kk2*180)/pi; 膝関節角度＝３６０－（大腿部セグメント角度－下腿部セグメント角度） LKNEag=360-(LKNEag1-LKNEag2); 伸展位１８０度 問題 Ｑ１ 足関節角度を求める Ｑ２ 股関節角度を求める
１９．補間作業 規格化する（３次スプライン補間） 規格化したい範囲を１００等分する 今回はすべての範囲を１００％に規格化する MaxT=max(Time); 最終フレームの時間を求める Inter=MaxT/100; その時間を１００等分する STime=0:Inter:MaxT; ０～１００％までの時間を作る LKNEags=spline(Time, LKNEag,STime); 規格化した結果=spline(元の時間軸,規格化したい変数、規格化時の時間軸) 直線補間 地面反力を補間するときはこちらの方がよい場合もある LKNEagL=interp1(Time,LKNEag,STime,'linear'); ２０．地面反力のサンプリングを下げる
