matlab課堂筆記 061215_內插法

內插法 interp1
當查表時，如果欲查的值不在表格內，可用interp1求得
舉例：

>> x=0:10;y=sin(x);[x' y']

ans =
         0         0
    1.0000    0.8415
    2.0000    0.9093
    3.0000    0.1411
    4.0000   -0.7568
    5.0000   -0.9589
    6.0000   -0.2794
    7.0000    0.6570
    8.0000    0.9894
    9.0000    0.4121
   10.0000   -0.5440

>> y1=interp1(x,y,4.5)

y1 =   
     -0.8579

也可以由y值反求x值

>> x1=interp1(y,x,0.8)
x1 = 
     2.3488

Topic: 人體受外力牽拉後引發之姿勢反應_肌肉活化情形

(Horak et al., 1992)






前言：
當人體在站立姿勢下接受一外在的平衡干擾時，下肢與軀幹肌肉的活化程度可以反映人體受到外在干擾後引發之姿勢反應。而當人體在站姿下受到外力干擾而產生前向的身體擺盪時，背側的肌肉主要負責活化以維持或回復平衡，避免跌倒發生。過去文獻指出，當站立於一塊會產生水平後向移動之力板上時，隨著干擾強度增加(力板移動速度加快或移動位移增加)，受試者肌肉活化程度亦會隨之增加。但，由腳底板所接受之干擾與從腰間承受之外在干擾對於受試者的姿勢干擾程度與形式是有所不同的，因此，當受試者由腰間接受干擾後，其肌肉活化情形是否亦會隨著干擾強度增加而增加則值得探討。


實驗目的：
探討人體在站姿下從腰部位置受到之前向牽拉干擾時背側肌肉之活化情形，且比較在不同牽拉強度干擾下，人體反應是否會隨之變化。



實驗流程：
在受試者肚臍位置穿戴一腰帶，將三牽拉軸之牽拉線固定於受試者腰帶上(如上圖)，由三牽拉軸產生之合力提供為干擾受試者站立平衡之牽拉力量。受試者在受到牽拉之後必須在不跨步的前提下盡可能地維持平衡不跌倒。牽拉強度共有4種 (固定牽拉速度為9cm/s，牽拉位移則依序增加：3,5,7,9cm)，每個牽拉強度情境各進行3次測試。以肌電圖所量測到的肌肉活化情形作為受試者受牽拉後的姿勢反應依據，受測肌肉包括：腓腸肌 (MG)、膕旁肌(HAM)、臀中肌(GM)，以及脊旁肌(PARA)。

	牽拉位移 (cm)	牽拉速度 (cm/s)
V1	3	9
V2	5	9
V3	7	9
V4	9	9

MATLAB程式目的：
1. 由既訂之檔案途徑擷取資料：一為原始檔案(包括時間與肌肉活化大小(以voltage表示))，一為選取後之onset時間(牽拉開始以及各肌肉開始活化時間)。
2.
功能一(mode 1)：由原始EMG資料(excel檔案)計算不同牽拉強度所引發之integral EMG大小(IEMG計算方式為，肌肉開始活化後150~500毫秒減去牽拉開始前150~500毫秒時間內之IEMG大小做標準化)。並計算各個牽拉強度下各受測肌肉IEMG之平均值與標準差，以excel格式輸出並儲存於既定之資料夾內。
功能二(mode 2)：畫出各個牽拉強度下EMG活化情形(各以其中一個具代表性的trial呈現)，並儲存於既定之資料夾內。
Flowchart


以下為部分完成之程式內容
程式內容：

% IEMG.m

% function:IEMG calculation and EMG figure for one trial of a subject

%          onset data were inputted by key-in

%          IEMG results and EMG figure were displayed on the command window% raw_data: 1_Time(ms) 2_Fy 3_Mx 4_MGr 5_HAMr 6_PARAr 7_GMr 8_Pullonset

%           range from (600ms prior to the pull onset) to (5000ms after the pull onset)

% baseline EMG:time interval between (500ms prior to the pull onset) to (150ms priot to the pull onset)

% IEMG:time interval between (150ms after the burst onset) to (500ms after the burst onset)

pullonset=input('please input the onset time of the perturbation (ms):');

MGronset=input('please input the onset time of the Medial Gastrocnemius muscle activation (ms):');

HAMronset=input('please input the onset time of the Hamstring muscle activation (ms):');

PARAronset=input('please input the onset time of the Paraspinal muscle activation (ms):');

GMronset=input('please input the onset time of the Gluteus Medius muscle activation (ms):');

raw_data=xlsread('D:\Physical Therapy\Graduate School\Thesis\Pilot study for Thesis\Data\ASCII output\Selected channels_Fy_Mx_rMG_rHAM_rPARA_rGM_pull onset\He01C1V1R1.xls','C19:J5618');

EMGonset=[MGronset,HAMronset,PARAronset,GMronset];

baseline=mean(raw_data(101:451,4:7));

for i=1:4

    rawIEMG(i)=(1/2*(raw_data((EMGonset(i)-pullonset+149),i+3)+raw_data((EMGonset(i)-pullonset+499),i+3))+sum(raw_data((EMGonset(i)-pullonset+150):(EMGonset(i)-pullonset+498),i+3))); % unit:V*ms

end

IEMG=rawIEMG-baseline;

IEMG_MGr=IEMG(1),IEMG_HAMr=IEMG(2),IEMG_PARAr=IEMG(3),IEMG_GMr=IEMG(4)



% 以下為繪圖指令

muscle_name = {'MG (V)' 'HAM (V)' 'PARA (V)' 'GM (V)'};

linecolor=['b' 'r' 'g' 'm'];

for i=1:4

    subplot(4,1,i),fig_EMG=plot((raw_data(:,1)),(raw_data(:,i+3)),linecolor(i));xlabel('Time (ms)');box off;hold on;onsetline=plot(pullonset*ones(100,1),0:max(raw_data(:,i+3))/99:max(raw_data(:,i+3)),'k');text(pullonset,max(raw_data(:,i+3)),'pull onset');

    ylabel_command = ['ylabel(''',muscle_name{i},''')'];

    eval(ylabel_command)

end

執行結果：

IEMG數值


EMG圖形


討論：
1. 實驗結果的蒐集會依實驗情境不同而讓受測者進行多個測試，所以結果大多會以多個檔案分別儲存(以測試名稱做區別)。程式的目的主要是想將實驗結果以最簡便的方式做運算，不須經過太多人為的介入，所以程式的重點將大部分著重在資料的處理部份，減少對話框方式的進行，讓使用者可以在同一時間內處理多筆資料。
2. 目前的程式內容只單純建立IEMG計算方式以及EMG活化之圖形呈現格式，且牽拉和肌肉活化onset時間是以對話方式呈現。之後待完成的部分還包括：
(1) 將onset時間輸入方式改為由程式直接從已建立之excel檔案內擷取所需資料。
(2) 繼續完成在同時間內讀取多個trial之測試結果，並將各個trial之IEMG計算值和每3個trial之平均值和標準差(代表各個牽拉強度)之執行結果以excel檔案形式輸出並儲存。
(3) 選出欲代表的EMG活化圖形之trial(以對話方式進行)，並將圖形輸出且儲存於既定資料夾內

matlab課堂筆記061207_稀疏矩陣

稀疏矩陣的應用 sparse精簡一般含有零元素甚多的矩陣，僅就非零元素之位置加以陳述。可以將一個含有許多零的元素矩陣改寫為位址表示，但其實質仍如原來之矩陣，功能相同 (其大小仍然與原來之Ｍ矩陣相同，而且可以直接與其他矩陣相互運算)。可使用full指令復原。
S = sparse(A) (A為矩陣之變數)
S = full(A)

稀疏函數 spfun
f = spfun(@fun,S)
對稀疏矩陣S之非零元素作函數(@fun)取代

MATLAB課堂作業：學生期中考成績 (敘述統計_資料分佈)

前言：
敘述統計(descriptive statistic)能夠初步呈現整體資料情形，透過直方圖圖形的表示方法讓人更能清楚地了解資料之分佈情況。以下，以"學生期中考成績"做為範例。

老師網頁程式如下：
原程式

程式指令：

% statistic.m
% Degree distribution
% Data Input
n=0;
while 1
  A=input('輸入重量資料(按enter輸入結束): '); % '請輸入學生期中考成績(按enter輸入結束):'
  if isempty(A), break;end
  n=n+1;Y(n)=A;
end
%
disp(['==>總計輸入 ' num2str(n) ' 筆資料.'])
disp(['最小值：' num2str(min(Y)) '      最大值：' num2str(max(Y))])
fprintf('平均值：%0.2f  標準差：%0.2f\n',mean(Y),std(Y))
% 加入平均值和標準差的計算，並且設定數值的格式為小數點後兩位，所以改用fprintf指令
while 1
    Rk=input('輸入初始點及級矩(以空格分開)','s');
    if ~isempty(Rk),
     rank=str2num(Rk);
     r1=rank(1);r2=rank(2);
     break; % 如果按沒有輸入直接enter的畫會重複出現問句
    end
end
D=zeros(fix((max(Y)-min(Y))/r2)+2,1); % fix指令為捨去小數部分，D矩陣為結果中的frequency部分。原本只須加1，但是因為結果中有<的部分所以加2
for i=1:n
  if Y(i)<r1,
      D(1)=D(1)+1;
  else
      J=fix((Y(i)-r1)/r2)+2;D(J)=D(J)+1;
  end
end
disp(['  <' num2str(r1)  '  ' num2str(D(1))]);
for i=2:length(D)
  disp([num2str(r1+(i-2)*r2) '-'...
      num2str(r1+((i-1)*r2))  '  ' num2str(D(i))])
end
% 以bar指令畫出資料分佈頻率之圖形 (長條圖)
H=zeros(fix((max(Y)-min(Y))/r2)+2,1);
H(1)=r1-(r2)/2;
for j=2:length(D)
    H(j)=r1+(j-2)*r2+(r2)/2;
end
bar(H,D)
xlabel('Value') % '期中考分數'
ylabel('Frequency') % '學生個數'
title('Distribution of the Weight')% 'T大xx系學生期中考成績分佈圖'

執行結果：
統計結果
應用：學生期中考成績

討論：
1. 加入平均值和標準差資料更能清楚呈現資料。
2. 此程式中，若輸入值皆為正值，而起始點輸入為0時，程式執行會出現錯誤，這一部分需再修改。

matlab課堂筆記 061201_第十章 聯立方程式解

行列式與反矩陣
AA^-1=I
A^-1A=I
Ｉ稱為單位矩陣(Identity matrix)，為一方矩陣，且其對角線元素值均為1。以eye指令可製作此矩陣。
矩陣Ａ之反矩陣(A^-1)以inv(A)表示。Ａ矩陣之行列式值需不得為零，否則其反矩陣不能存在 (以det(A)作測試，計算其行列式值，若det(A)=0，則矩陣A為奇異矩陣)。

線性聯立方程式矩陣解法
聯立線性方程式可以化成ＡＸ＝Ｃ之型式：
3x +4y =10
5x -2y =8

(1) 左除法求解：利用MATLAB求解時，只要用矩陣倒除即可，或稱為左除法。即X=A\C：
A=[3 4;5 -2];
C=[10;8];
X=A\C
X =
2 1
(2) 反矩陣法求解：A^-1AX=A^-1C 　或A^-1AX=IX=X=A^-1C
A=[3 4;5 -2];
C=[10;8];
X=inv(A)*C
X =
2.0000 1.0000
答案和用左除法所求得的解相同，但是左除法較常見。