這原本是我三年前剛畢業的時候想要完成的計畫，但是一直都沒有時間完成。這次一來剛好利用這次要分析數據的機會，二來因為之前分析的數據實在太多了，多到我都記不得有寫過哪些程式，分析過哪些數據，所以藉此機會除了把它整合在一起外，順便將實驗室過去用到的數據整理起來。
檔案下載：
請按此 （密碼：同影印機）
以下為使用教學，我知道整個介面很醜，但是我實在很不喜歡修飾，等有空我再來修飾，功能可以用比較重要吧！
首先要先建立一個資料夾，這個程式的階層是這樣的：一個
主要資料夾（這邊我舉例的主資料夾名稱為
整合計畫test，裡面一定要自己設定一個叫data的資料夾），然後把同樣一隻腳的測試資料丟到data資料夾，這個程式會自動去尋找主資料夾裡的data資料夾，讀裡面的資料
主資料夾
- data
- xxxx
- xxxx
- xxxx
然後在matlab的command視窗打入main_program（記得上面的current directory要改成剛剛下載好的檔案存放的路徑）。
接下來就出現選單，可以選擇要計算的參數（還真多啊！）。
這裡我就示範全部勾選，一口氣用這兩筆做ITT的資料跑出全部數據！
接下來會出現一個視窗（我知道真的很醜！），可以填寫姓名，左右腳，以及腳長，主要是在最後輸出檔案的時候上面可以註明這個資料的主人，但是這僅限於第一次用，因為這邊的功能以及資料的寫入都是用附加的方式，所以只有在第一次跑的時候因為檔案尚未存在而會寫入這些基本資料（我知道我解釋的很難懂！）。
選擇剛剛自己建立的
主資料夾（此例就是整合計畫test），不要選到data資料夾！！
因為有些計算須要用到force onset，但是依我算了不算多的數據經驗發現，論文裡面所寫的用兩倍標準差當作開始根本就是一個理想，像我現在的例子就可以看到這個受試者居然在用力那瞬間力量突然往下掉（可能想利用stretch-shortening cycle的方式看可不可以力量大一點，所以先勾了一下），或者有些受試者在喊預備的時候就在那邊偷偷用力準備（這些不測風雲有時難以預測，有時候都不是當下可以發現到），所以這邊提供了輸了標準差的方式。
因為計算RFD以及RER須要換算成torque，所以都有勾選要計算這兩個參數的時候，就會出現要求輸入腳踝的長度（力臂）。
如果計算active EMD的話，會強迫重新輸入SD，至於原因則是
不能說的秘密！
中間會出現一大堆運算後的圖片，這些都可以在之後主資料夾裡面找到（圖中為，百分比MVC對於時間的變化率）。
如果有勾選curveRMS這個參數（主要是計算不同百分比MVC下，EMG RMS對時間變化的斜率），會發現突然停頓，這是因為計算RMS所需的時間較長。
這邊也會要求輸入RMS EMG的SD，原因跟force一樣。
劈哩啪啦跑完之後，可以發現主資料夾多了很多資料夾，以及有一個叫做
result.txt的檔案。
我們用一些文字編輯器或是notepad打開可以看到我們剛剛算的數據。
由於文字編輯器對於換格tab會因為字串長度不同可能無法對齊，導致看起來不是那麼好看，所以我們改用excel開啟看看。
由excel開啟之後，可以發現每個data都排排站站好。
而且每個參數都會把每次trial所計算的記錄下來（如果今天發生某次trial的值過於偏差，可以在此修改），以及這些trial的平均（有幾個檔案就用幾個平均，兩個檔案就兩個平均，三個就三個平均，一個當然就一個自己平均啦！真是太神奇啦！）！
然後我們打開資料夾看看，這邊是打開ITT資料夾，就可以看到這次所算的ITT圖以及結果。
當然如果要重算某樣參數也是可以，在不刪除result.txt的情況下，會自動把新算的數據直接附加到原來數據的最後，而資料夾裡面的圖片也自動被新的附蓋過去。當然如果想要重新跑的話，就把result.txt砍了就行了！
後記：
1. 這次的整合程式，主要的應用在我們做ITT的trial上面，也就是每次我都會跟受試者說要被電三次，我想如果都是照著這樣既定的protocol做的話，程式大多應該不會出現問題。
2. 我將原始檔案附在上面，希望大家使用上如果發生問題可以告知我（至少目前我的幾筆資料都是沒有問題），大家的反應才是可以把這個程式修改到好的原動力。
3. 這次有些程式我有改寫（變得比較精簡），然後整合程式撰寫我也是力求工整，這樣也比較利於之後如果有人想要增加新功能，新參數使用。
4. 有時間的話我再把他外觀弄得好看一點。
5. 最後，我該睡了！！
前言：改良自碩士班的原始版本，當初懵懵懂懂寫得很亂，現在以新觀念方法改造。
- Load File：開啟後選取圖片資料夾，資料夾內圖片會以[sort_name] = sort_file(path)排列，目的是將連續圖片排序，若是一般則以sort的方式（排列字母大小等）。但是當初因為自己的失誤把檔名寫成1.jpg, 2.jpg ….. 10.jpg, 11.jpg…..，經過排列後會變成1.jpg, 10.jpg, 11.jpg, ……, 19.jpg, 2.jpg, 21.jpg…..，因此才多寫若是檔名不一樣長，則以不同cell的方式將相同長度檔名放入排列（目前只有到十位數）。
- Sliderbar：會將sort過得file按照順序排列，提供一張張點選。
- position select：則是以sono_mousedown.m方式，將滑鼠點所點的點的座標顯示。
- Dist tool：以matlab內建的dist tool作用，不過目前沒有校正顯示的長度。
- Save Figure & value：將座標值存入table中，並且將點選好的圖片存入目前資料夾內，並且以[path new_filename] = changefilename(filename,add_word)改名（如：a_原始檔名）。改名之後，由於儲存圖檔的時候當初以plot方式所畫得十字記號無法存入，因此以[change_image] = draw_image(im,x,y)方式，將所選取的座標位置上下左右15個像素統統改成紅色（即變成十字）。
匯入座標資料於table，並且[cal_dist]=calcu_dist(ini_p,cal_p,dist_prop)，依目前深度計算和第一點的距離。 
- Export：將資料匯出，內建以目前資料夾為匯出處，以[file,path]=uiputfile('.txt','Save file as',eval('[handles.path default_file]'))完成想要匯出的資料夾，預設以dist.txt作為預設匯出檔名。
匯出後以sonofilemerge(path,sonofile,distfile)自動和sono.txt檔案結合成stiff.txt（個人需要）。
這個東西我實在沒有時間好好的把他製作一番。
electroGUI V.0.0001（毫無意義版）界面介紹
- 開啟（以autochoosetype.m）：目前內建可讀全部為數據的txt檔（以textread開啟），InstruNet（以readintranet.m開啟）所產生的.TXT檔（如：Ch1_Vin+.TXT），以及一般具有標頭檔案（以hdrload.m開啟）。
- 資料數據（以datacursor）：則是可以點選圖形，顯示原始x,y資料點，如下圖：
- 選取範圍（以get_x_point.m）：即點選圖形區域，目前規定要由左到右選取。如下圖：點選後按peak to peak即以所填入的放大倍率計算peak to peak amplitude。
若是資料選取順序錯誤，則會出現警告。
結論：
目前版本非常陽春而且bugs很多，待我有空有時間，在陸續修正並增加新功能，不過對於計算peak to peak amplitude目前來說應已經足夠。
當我們將肌電訊號轉換成頻譜（Power spectrum）後，我們的確可
以看到頻率和其強度的分佈，不過這代表什麼意義呢？
我們知道我們的肌肉由不同型態組成，簡單的說可以分為快肌
（fast twitch）跟慢肌（slow twitch），快肌肌肉纖維較粗，動作速
度快，不過較容易疲勞，慢肌動作速度慢，但是可以作用持續時間長，
就傳遞肌肉表面膜電位而言，快肌就像比較粗的電線，傳遞訊號的速度
較快，因此由電極片所接收到單位時間內的電位差訊號相對而言就會比
較多，所以頻率會落在比較高的位置，反觀慢肌，則為相反地現象。
而我們的肌肉，則由不同比例的快慢肌所組成，例如我們強而有力
的肱二頭肌（biceps），快肌的組成較多，所以也比較容易酸。
不過因為頻譜看出來的是各個頻率強度的變化，如果要相互比較的
話，則必須要像一個常態分佈圖形一樣，有平均數（mean），等數值，
比較常用來作為計算的為中位頻率（Median frequency）。
中位頻率（Median frequency），其定義為頻譜分析圖中，將其強
度由低頻積分到高頻，積分面積為全部面積一半時的頻率。因此就可以
中位頻率來作為訓練或是疲勞測試後，改變肌肉組成比例（訓練）或是
影響肌肉放電（疲勞）等改變觀察。
因此利用傅立葉轉換以及if指令來求得中位頻率，當積分到某一點
時，大於全部等積分的一半時的頻率即是中位頻率。
由於肌電訊號是時間性的，我們比較常做的方式會像之前所用的以
window的方式，每多少點計算一次，重複多少點，然而取樣頻率
（sampling rate）會影響最後每一個點所代表的頻率，所以也必須考慮
。
function [med_fq]=fft_i(data,n,op,sr)
% fft for n point
% op: the point of overlapping
% med_fq is the median frequency of the power spectrum
% sr is the sampling frequency
% calculate the times of the loop (-1 due to the initial loop is 0)
l=fix(length(data)/op-1);
fq=(1:n/2).*sr/n; % calculate the frequency of each point
for i = 0:(l-1)
int=1+op*i; % the initial value of the loop
f=fft(data(int:(int+511))); % fft for the data
mag=abs(f); % Amplitude of the fft
mag=mag(1:n/2);
for k=1:n/2
if sum(mag(1:k)) >= sum(mag)/2
med=k;
break
end
end
med_fq(i+1,1)=k*sr/n;
end
% calculate the final part of the sample
int=1+op*l;
f=fft(data(int:end)); % fft for the data
mag=abs(f); % Amplitude of the fft
mag=mag(1:n/2);
for k=1:n/2
if sum(mag(1:k)) >= sum(mag)/2
med=k;
break
end
end
med_fq(l+1,1)=k*sr/n;
首先必須先計算迴圈的次數，以所有的點數除上重複的點數（通常
會以window點數的一半）減一（種樹法則），然後就每多少點作一次傅
立葉轉換並且求得中位頻率，值得注意的是，因為常常所選取的資料數
並不能剛好夠最後一筆，舉例而言，也就是說原本應該要以512個點作計
算，但是迴圈到最後一次時可能只剩下500點，由於fft這個函數當不夠
點的時候會比剩下的點以零帶入，所以最後一個窗口則以另外方式計算
。
執行結果
>> data=textread('D:/EMG_fft.txt');
>> plot(data);
>> med_fq=fft_i(data,512,256,1000);
>> figure(2);plot(med_fq)
這是原始資料圖，資料可以
由此下載 經由計算後為
一開始很高可能是因為在還沒有任何訊號時就已經存在的高頻雜訊，
此時就要利用之前提得濾波器給予濾波。
肌電訊號分析可以分為兩個主要部份，一個即是時域
（time domain）另一則是頻域（frequency domain）。
大部分例如方均根植（root mean square）或是找尋
起始點（onset point），是屬於時域的部份；然而頻域的
部份會著重去看頻率和強度之間的關係（power spectrum）
，而以下所介紹的為利用傅立葉轉換（Fourier transform）
所做出的關係圖。
程式如下
data=textread('D:\EMG.txt');
fft_data=fft(data);
mag=abs(fft_data);
fq=(1:length(data)/2).*1200/length(data); % 每一點的所表示的頻率
mag=mag(1:length(data)/2); % 因為另一半為image
plot(fq,mag)
首先可以發現到當對傅立葉轉換後的圖形如下，每個
點都是有相位與長度大小的點。所做出為複數座標系的圖
。
不過對於分析強度和頻率而言，只要去看轉換後每個
點的強度即可，所以將傅立葉轉換後的值取絕對值，即是
強度。
作圖後可以發現到，會是一個以中線分隔兩邊成左右
對稱的圖形也表示只有在一半以下的點是有意義的。
至於為什麼只有一半是有意義而另一半會是跟前一半
對稱鏡射，這是由於我們取樣頻率的關係。
正確的取樣頻率必須至少是原來訊號的兩倍，這就是
所謂的Nyquist rate。即是Nyquist Theorem。
由下圖可以見到，如果取樣頻率較小，原始訊號頻率
很高，所取的訊號圖形和原始訊號就有差距，而變成較低
頻的訊號。
而每一個點所代表的頻率則為”取樣頻率／點數”，
由於只有一半的值是有意義的，所以須將頻率以及強度作
處理，如下：
fq=(1:length(data)/2).*1200/length(data); % 每一點的所表示的頻率
mag=mag(1:length(data)/2); % 因為另一半為image
因此做出來的圖形如下
一般我們在分析肌電訊號的時候簡單方式，讓訊號有如封包(envelope）一樣，或是利用這些方式作為分析比較。
常用的有，一階的低通濾波，方均根（rootmean square），或是平均的方式以及積分的方式。然而由於取樣頻率高，因此我們不可能把所選取的範圍全部一次作處理，而會利用窗(window）的方式，其意義是，每次選取特定的點作處理。
依據這樣的概念，利用for loop的方式，也就是每次選取固定的點作為分析使用。
舉例而言，下圖是一段EMG的訊號。
以下是所撰寫的程式碼
function [out point]=moving_avg(data,dur)
l=length(data);
k=menu('Select the function','Average','RMS');
switch k
case 1
for i=1:(length(data)-dur);
out(i,1)=mean(abs(data(i:i+dur,1)));
point(i,1)=(i+i+dur)/2;
end
case 2
for i=1:(length(data)-dur);
out(i,1)=sqrt(sum(data(i:i+dur,1).^2)/dur);
point(i,1)=(i+i+dur)/2;
end
end
使用的情形如下
目前以menu方式寫成兩種處理方式。point是說明每一次作處理後點應該落的位置（因為10作完後只有一個點）。
其結果和原始訊號圖畫在一起如下
首先我必須先聲明的是，我對於濾波的真正原理是
一知半解，但是他在我們物理治療中的訊號分析裡面，
卻是扮演著重要的角色。
什麼事butterworth的濾波，說真的我真的不知道，
反正數位濾波器就以這樣設計就好。
看到這裡如果你非常想打我的話，元智大學最佳化
實驗室有份文件大家可以參考看看。
按我看看 言歸正傳，反正濾波器有高通，低通，帶通，帶拒
四大類，於是乎我利用了matlab裡面的butter函數，以
及filter設計了這個含有四種功能的濾波器function，
內容很簡單，只不過利用了switch的方式選擇。
點我看看這個簡易的濾波功能 如果有興趣試試看的話，這裡提供了一筆資料可以
試試看
請以右鍵點選下載。
這分資料裡面，第一列跟最後一列分別都是force的
資料，中間幾個才是肌電圖，其取樣頻率為1000/s。
在我們物理治療研究中，不管是研究肌械延遲（electrical mechanical delay），或是病人對於肌肉張力，運動員之力量發展，都需要靠肌電訊號來告訴我們，其中有一個最重要的部份就是，去知道該塊肌肉什麼時候開始動作。
下載檔案（點右鍵下載此檔案） 其中，前三條分別為我們下肢小腿後側肌群。最後一條為實驗時我用數位類比轉換卡發出的類比訊號給儀器接收，以作為同步訊號。不過目前先不將同步訊號考慮在內。
我們必須要去找到開始的時間點，此時就要利用到統計學的概念，兩個母群若是差兩個標準差的話，就可以稱作是兩個不同的母群（否決虛無假設）。
不過很不幸的是，通常訊號可能沒有那麼漂亮，或是雜訊很多，可能需要到三個甚至四個標準差以上。
點此下載看上圖程式檔案 這個程式就是靠while+break所寫成的，當找到符合條件的點後，就結束迴圈，並且輸出結果。
