Matlab 講習会 目的 Matlab を用いて VICONや Winanalyze の座標データー 地面反力の分析必要な項目について習得する 本やヘルプに掲載されている情報を 実際に使用できる形で整理する 講習会 1 回目 (4 時間 ) 1. 行列操作について理解する 2. 時間軸を作る 3. エクセルデーターを取り込む 4. テキストデーターを取り込む 5. グラフの作成 6.1つのグラフに複数のグラフを出す 7. 計算によって新しい座標軸を作る 講習会 2 回目 (6 時間 ) 8. フィルターをかける 9. 速度を算出する 10. 合成速度を算出する 11. 最大値 最小値を求める 12. 時系列データをまとめる 13. データを出力する 14. 図表を出力する 15. スティックピクチャーを書く & 軌跡を描く 16. セグメントを書く &axis equal 縦横比をそろえる 17. 画面に出力するスティックピクチャーの数を減らす 18. 関節角度を算出する 19. 補間作業 20. 地面反力のサンプリングを下げる 1
1. 行列操作について理解する 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 行 列 1 行目だけ取り出す A(1,:) 1 列目だけ取り出す A(:,1) 問題 Q1 3 行目だけ取り出す Q2 3 列目だけ取り出す 2~3 列目を取り出す A(:,2:3) 2~3 行列を取り出す A(2:3,2:3) 2
2. 時間軸を作る 時間軸を作る 1 秒間 120 コマであったら 1/120=0.0083 Time=0:0.0083:1 101 行のデータに,0.0083 刻みの時間軸を作るには DA=0:100 このままだと101の列行列になる 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 べき乗 2.^3 2の3 乗 power(2,3) 2の3 乗 ゼロの行列がほしい場合 ZZZ=zeros(10,1) 1の行列がほしい場合 OOO=Ones(1,10) 3
3. エクセルデーターを取り込む clear all [DATA HEAD]=xlsread('Matlab1.xls'); と入力する DATA には 数値データが入る HEAD には ヘッダーが入る ( 重複ヘッダーがないか注意 ) 取り込んだデータを マーカーごとに変数にする II=strmatch('LASI:X,HEAD) 探したい文字は で囲むようにする HEAD( ヘッダー ) のなかに LASI:X が何行目にあるか探してくる (4 行目 ) LASI=DATA(:,II) (Field=DATA(:,4) と入力しているのと同じ ) 行はすべてほしいので : と入力する X Y Z のすべてほしい場合は LASI=DATA(:,ll:II+2) (Field=DATA(:,4:6) と入力しているのと同じ ) VICON を用いているときは 3 行ごと取り込むようにしている 慣れない場合は1 行ごと取り込む ( 計算によっては1 行ごと取り込んだ方が計算しやすいことも ) 位置座標の単位を ミリからメーターに変える ( 線速度などを出力するとき必要 ) LASI =DATA(:,ll:II+2)/1000 4
4. テキストデーターを取り込む LANK=dlmread('Read1-1.txt',' t',1,0); ヘッダーの部分は取り込めないので 1 と入力して2 行目以降を取り込む KNEE=dlmread('Read1-2.txt',' t',1,1); と入力すると 2 列目以降が取り込まれる LASI=dlmread('Read1-3.txt',' t',1,0); 5. グラフの作成 位置座標をプロットする 時間軸を作る LL=length(KNEE) TimeA=0:0.004:0.004*(LL-1) figure(1) plot(timea,knee(:,1),timea,knee(:,2),timea,knee(:,3)) タイトルを入れる X 軸タイトル Y 軸タイトル title( Knee ) で文字を挟む xlabel( Time(sec) ) ylabel( mm ) 凡例 ( はんれい ) legend( X, Y, Z ) plot の順番と対応させて書く 5
表示の範囲を設定する axis([0 0.15 -Inf Inf]) ( ) の中に [ ] を入力する はじめの2つが Xの最小値 最大値 残りの2つがYの最小値 最大値 軸の範囲を設定しない場合は 最小値には-Inf 最大値には Inf と入力 プロットの線 マークの種類を変更する線点 ( プロットシンボル ) 色の順で配列する k 黒 b 青 r 赤 g 緑 c シアン m マゼンダ X バツ O まる +プラス * アスタリクス s 正方形 dダイアモンド v 三角形 - 実線 : 点線 ; 鎖線 -- 点線 マーカーサイズの設定方法 plot(x,y, -ok, MarKerSize,2) MarkerSize 数字を変えることによって マーカーのサイズが変化する plot(time,lasi,'ok',time,lasif,'-k','markersize',2) 6
6.1 つのグラフに複数のグラフを出す 上下に並べる 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) 4 分割する 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) 7
7. 計算によって新しい座標軸を作る 両肩峰を結んだ点を作る (A-B)/2+B CSHO=((RSHO-LSHO)/2)+LSHO; 8. フィルターをかける フィルターの作成 [b,a]=butter(2,10/60); 4 次のバターワースフィルター ( ローパスフィルター ) 遮断周波数 /( ナイキスト周波数分析に用いた周波数 /2) フィルターをかける [b,a]=butter(2,10/60); LASIf=filtfilt(b,a,LASI); 8
9. 速度を算出する 位置座標 角度データ フィルター 速度 角速度計算 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]; 10. 合成速度を算出する CGV= sqrt(power(comassv(:,3),2)+power(comassv(:,2),2),power(comassv(:,1),2)); 9
11. 最大値 最小値を求める [CGVmax CGVmaxTime]=max(CGV) CGVmax 最大値 CGVmaxTime 最大値出現時のフレーム [CGVmax CGVmaxTime]=min(CGV) 最小値を求める 12. 時系列データをまとめる EX=[COMass LANK LASI Force Field] 13. データを出力する wk1write('matlab1',ex,1,0) 1 行目をあけたい場合 ( ヘッダーを作る場合 ) は 1,0 と入力する 14. 図表を出力する saveas(figure(1),'matlab1','tif') bmp Windows ビットマップ jpg JPEG イメージ tif TIFF イメージ 10
15. スティックピクチャーを書く 図を重ねあわせるコマンド (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 11
16. セグメントを書く 結線したい A と B を X 座標では [(AX),(BX)] Y 座標では [(AY),(BY)] と入力する plot( X, Y ) plot( [ (AX),(BX) ], [ (AY),(BY) ] X 座標 Y 座標 下腿のセグメントを作る 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 axis equal Q3 肩峰をプロットする Q4 肩峰と骨盤を結線する 12
17. 画面に出力するスティックピクチャーの数を減らす 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 13
18. 関節角度を算出する 関節角度は 2つのセグメント角度を求め その2つの角度を加算減算することによって求める 膝であれば 地面に対する大腿部セグメント角度と下腿部セグメント角度を求め その2つを加算して 360 度から引くと関節角度が算出される 1 動くマーカーの座標を求める ( 原点に対して動くマーカー )-( 原点になるマーカー ) 2 動くマーカーの座標を 関数 atan2( アークタンジェンド ) に入力する atan(y 座標, X 座標 ) 14
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; 膝関節角度 =360-( 大腿部セグメント角度 - 下腿部セグメント角度 ) LKNEag=360-(LKNEag1-LKNEag2); 伸展位 180 度 問題 Q1 足関節角度を求める Q2 股関節角度を求める 15
19. 補間作業 規格化する (3 次スプライン補間 ) 規格化したい範囲を 100 等分する 今回はすべての範囲を 100% に規格化する MaxT=max(Time); 最終フレームの時間を求める Inter=MaxT/100; その時間を100 等分する STime=0:Inter:MaxT; 0~100% までの時間を作る LKNEags=spline(Time, LKNEag,STime); 規格化した結果 =spline( 元の時間軸, 規格化したい変数 規格化時の時間軸 ) 直線補間地面反力を補間するときはこちらの方がよい場合もある LKNEagL=interp1(Time,LKNEag,STime,'linear'); 20. 地面反力のサンプリングを下げる ForceZ=decimate(ForceZ,9); サンプリングを9 分の1にする 16