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

Transcription

1 スポーツのための Matlab 講習会 対象者 :VICON 利用を考えている学部生 大学院生 卒論 修士論文 博士論文のバイオメカニクスデータを分析するにあたって 多変量 多変数のパラメーターを分析するとなると エクセルを用いて分析をしていては 何ヶ月もかかって効率が悪い そこで 行列計算を行うことができる Matlab というプログラム言語を利用することで 短時間に分析することができ 研究の他の部分に時間をまわすことができる Matlab に関する本は多くあるが スポーツの分析に必要な方法がバラバラに記述されており 独学すると習得するために莫大な時間かかる そこで 効率よく Matlab の使い方を講義して 実践的に使えるようにするのが本講習会の目的である

2 1 行列の抽出 BMI= 体重 (Kg) ( 身長 (m) 身長 (m)) RJ_index= 跳躍高 (m)/ 接地時間 ( 秒 ) 2 跳躍高 ( 滞空時間から求める )=(1/8) 9.81( 重力加速度 m/s) 滞空時間 A=[1 2 3;4 5 6;7 8 9] A= 縦が 行 横が 列 行と列を区別できるようにする A=[1 2 3;;4 5 6;7 8 9] A=[1 2 3;4 5 6;7 8 9]; 式の最後に ; セミコロンをつけると画面表示されない 2 列目だけ抽出する QQ(:,2) 4 行目だけ抽出する QQ(4,:) : コロンをつけると すべての行 すべての列を抽出する 2 列目 ~5 列目を抽出する QQ(:,2:7) 抽出したい変数 ( はじめの列 : 終わりの列 ) 4 行目 ~7 列目を抽出する QQ(4:7,:) 抽出したい変数 ( はじめの行 : 終わりの行 ) mean 平均値 max 最大値 min 最小値 sum 合計 std 標準偏差 median 中央値 var 分散

3 2 行列の結合行列同士を結合するには 行列 A と行列 B の行の大きさ 列の大きさが同じでないと結合できない 行列 A=[ ] 行列 B=[ ] 行列を横に結合するには 変数間にスペースを空けて並べてから 括弧で囲む [ ] C=[A B] C=[ ] 行列を縦に並べるには 並べたい変数の間にセミコロン (;) を入れて 括弧で囲む [ ] C=[A;B] D=1:2:10 E=2:3:33 問題 D と E を横に並べてください E と D を縦に並べてください 行行列を列行列に変更する列行列を行行列に変更する D アポストロフィを付ける (Shift を押して 7 のキーを押す )

4 3エクセルデータの取り込み [DATA HEAD]=xlsread('Rotation.xls'); xlsread はエクセルのファイルを読み込む関数です エクセルファイルを読み込ますには エクセル 5.0/95 ブック の形式にする エクセルの最新のファイル形式は読み込まないので注意 ファイル名は で囲む 4グラフの作成 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 図表の外に凡例が表示される 1 図表の右上に表示される 2 図表の左上に表示される 3 図表の左下に表示される 4 図表の右下に表示される

5 5グラフの線 プロットのサイズの設定 plot(time,data(:,2), -ok,time,data(:,4), -vk,time,data(:,6), -sk,) X 軸データ Y 軸データ, グラフの線の種類 マーカの種類 線 マーカーの色の設定, X と Y のデータのあとに入力する マーカの種類 色 線種類 直線 - 丸 o 黒 k 点線 : プラス + 青色 b 鎖線 -- 四角 s 赤色 r 一点鎖線 -. 下三角 v 緑色 m 五角形 p 紫色 m X 軸 Y 軸の大きさを調整する axis([x 軸最小値 X 軸最大値 Y 軸最小値 Y 軸最大値 ]) 特に設定しないときは 最小値には-Inf 最大値には Inf を入力する 6 行列の行の大きさをもとめる LL=length(DATA) length は行 列で最も長い 行 列 を調べる関数です 7 行列の行と列の大きさを求める size(data) size は行の長さ 列の長さを同時に調べる関数です

6 8ヘッダーの検索 データーの抽出エクセルデータからヘッダーの位置を検索する II=strmatch('LASI:X',HEAD) strmatch を使って 何列目に必要なデータがあるか検索する関数 カッコの中は ( 検索したい文字列, 検索したい文字がある行列 ) の順に入力する 5 列目に必要なデータがあれば 5と表示される 5 列目に必要なデータがあることがわかったので 5 列目のデータを すべての行を抽出するように式を書く LASI=DATA(:,II) VICON のデータの場合は すべてのデータが X 軸 Y 軸 Z 軸の順に記述されているので :X の列番号がわかっていれば :X +1した列番号が :Y のデータになり :X +2した列番号が :Z のデータになる 3 列まとめて抽出すると LASI=DATA(:,II:II+2) のように書く必要がある VICON は座標データが ミリ単位 で出力されるので メートル単位に単位変換をする必要がある II=strmatch('LASI:X',HEAD); LASI=DATA2(:,II:II+2)/1000; 1000 で割ると ミリ単位からメートル単位に変換される

7 9 時間軸の作成 II=strmatch('Frame',HEAD); Frame=DATA2(:,II); LL=lenght(Frame) length 関数を用いて 行列の長さを調べる ( 行の長さを調べる ) Time=0:(1/120):(LL-1)*(1/120); 120Hz の時間軸作成 (LL-1) をしないと 実際の行数より多くなってしまうのでマイナス 1 する Time2=0:(1/250):(LL-1)*(1/250); 250Hz の時間軸作成 10subplot 関数の使い方 1つのグラフに複数のグラフに出力する 上下に並べる figure(1) subplot(2,1,1) plot(x,y,x1,y1) subplot(2,1,2) plot(x,y,x1,y1) subplot( 行 列 番号 ) を用いて 1 つの図表に複数の図表を表示する 何個の図表を表示させたいか決めて 行と列の大きさを入力する 2 1 行列であれば上部に表示させたければ1 下部の表示させたければ 2 と入力する

8 2 1 行列 1 2 行列 行列 3 2 行列 左右に並べる 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)

9 11 関節角度の算出関節角度を算出するには 上部セグメント角度 と 下部セグメント角度 の 2つを求めて その2つのパラメーターを用いて関節角度を算出する ベクトルの作り方中心となる関節を起点としてベクトルを2つ作る 終点 - 起点 = 起点を中心としたベクトル 大転子 - 膝 = 大腿部ベクトル外果 - 膝 = 下腿部ベクトル ( 課題 1 1 列行列ごとに作る課題 2 3 列行列ごとに作る )

10 角度算出セグメント角度を求めるには [atan2] という関数を用いて角度を算出する 角度 ( ラジアン単位 )=atan2( 縦軸, 横軸 ); 縦軸 横軸を間違えると算出できないので注意 角度単位変換ラジアン単位だとわかりづらいので Radian から Degree へ 角度単位を変換する 角度 =( 算出した角度 *180)*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; 12フィルタリングフィルタリングとは, 角度 角速度算出にあたって生じるノイズを取り除くために行われる 研究室によって使い方の考え方が異なるが 基本的には角度データにフィルターをかけたら 2 回フィルターをかけない 論文では位相ずれなしの4 次のバターワースローパスフィルターを 遮断周波数 Hz で行ったという記述でフィルタリングについて書かれてある まず バターワースフィルターの作成 ( 設計 ) [B,A]=butter( フィルターの次数, 遮断周波数 / [ 取り込んだデータの周波数 /2]); AG=filtfilt(B,A,[ フィルターをかけたいデータ ]); filtfilt を使用するのであれば 4 次のフィルターであれば 2 と設定する [B,A]=butter(2,10/60); AG=filtfilt(B,A,RKNE_ag);

11 13 速度を算出する位置座標 角度データ フィルター 速度 角速度計算 3 点微分法 ( スポーツバイオメカニクス阿江参照 ) 速度の計算方法 1 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]; 速度の計算方法 2 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)]; 250Hz の分析は 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)];

12 14 時間軸作成応用編膝の最大屈曲時を 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 つの行列を合成する

13 15Matlab 情報の検索 Matlab プログラムをしてわからなくなったら? 1. ヘルプを見る 2.Google で検索してみる Matlab ( 捜したいキーワード ) Yahoo は Matlab 関連の検索はあまりよくないので使っていない 3. サイバーネットシステムのサイトに行く 4. テクニカル FAQ で検索する 5.Mathwork 社 (Matlab の開発をしている ) のページに行く File Exchange に行くと 人が作ったファイルを手に入れることができる 6.Matlab の本を買って読んでみる マニア向け (BASIC や C 言語などの他のプログラムの発想を取り入れる ) 16Matlab の勉強法予備実験でもいいので実験を行って分析してみる 講習会の資料を基にして データーを分析してみる Matlab の本を一冊はじめから 終わりまで 式を書きながらやってみる David A.Winter. Biomechanics John Wiley & Sons,Inc を読んでみる わからない数学の部分は 数学の本を借りたり買ったりして理解する 卒論のプログラムを作ってみる ( はじめは 1 週間以上かかる )

14 17 ダウンサンプリング ダウンサンプリング : 地面反力が 1080Hz で取り込まれていて VICON のデータが 120Hz で取り込まれている場合 地面反力のデータを 9 分の 1 にして VICON のデータにあわせること サンプリングが小さいデータに サンプリングが大きいデータを データ数を少なくしてあわせること 規格化 :0.13 秒でリバウンドジャンプする人と 0.22 秒でリバウンドジャンプをする人の各関節の動きを比較したいときに 時間の長さを同じにして比較する方法

15 ダウンサンプリングの練習 clc close all clear all は入力しておく % データ取り込み [DATA HEAD]=xlsread(' '); エクセルのファイルを開く FDATA=csvread('yano12.csv',10,0); CSV ファイル ( カンマ区切り ) ファイルを開く 10,0 は 11 行目 1 列目からデータを抽出することを示している % ダウンサンプリング FF2=FDATA(:, : ); フォースプレート2のデータが必要なので抽出する QQ=length(FF2); FDATA をプロットしてください タイトル X 軸 Y 軸 凡例をつけてください axis([xmin Xmax Ymin Ymax]) [B,A]=butter(4,40/ ); フォースプレートのデータにノイズがあるので取り除くためにフィルターを作ります フォースプレート 1080Hz 1080 すると に入力する数字がわかります FF2=filtfilt(B,A,FF2); 位相ずれなしのバターワースフィルターをかける 1 基本の時間軸を作る Time=(0:(1/1080):(QQ-1)*(1/1080))'; サンプリングレート 1080Hz のデータ Time=0:0.001:2 2 変更したい時間軸を作る TimeS=(0:(1/120):(QQ-1)*(1/1080))'; サンプリングレート 120Hz のデータ TimeS=0:0.008:2 3ダウンサンプリングする出力 =spline(1, ダウンサンプリングしたいデータを入力,2); 1 行ごとしかダウンサンプリングできないので注意 FX2=spline(Time,FF2(:,1),TimeS); FX FY FZごとを合成して出力する

16 18 規格化 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); 規格化 ( 実際の場合 ) 1 分析したい区間の時間軸を新しく作る 2 分析したい区間を 100 等分する 上記の方法と同様 3spline(1 元の時間軸, データ ( 分析したい区間 ),2 新しい時間軸 ) 出題 : 右膝の最大屈曲時から最後のフレームまでを規格化する [RKNE_max RKNE_maxF]=max(RKNE_ExFx); 1 右膝の最大屈曲時のフレームを捜す BB=length(RKNE_ExFx); 2 最終フレームを捜す BL=BB-RKNE_maxF; 1スタートから2 終了フレームまでに 何フレームある求める TimeBB=(0:(1/120):BL*(1/120))'; 新しい時間軸を作る 新しい時間軸を 100 等分する MaxB=max(TimeBB); InterB=MaxB/100; TimeCC=(0:InterB:MaxB)';

17 3 次のスプライン関数を用いて規格化する ANKLE_SP=spline(TimeBB,RKNE_ExFx(RKNE_maxF:BB),TimeCC); 規格化したい区間だけを計算にかける必要がある PTime=0:100; 規格化用の時間直線補間地面反力を補間するときはこちらの方がよい場合もある SP_EX=interp1(Time,DATA,STime,'linear'); 19グラフを出力する saveas(figure(1),'yano12','tif') 保存したいグラフが figure(2) であれば figure(2) と入力する 20データを出力する 計算したデータを エクセルファイルに出力する 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) 1 行目をあけたい場合 ( ヘッダーを作る場合 ) は 1,0 と入力する