Labview and Matlab for PT

　　當我們將肌電訊號轉換成頻譜（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作完後只有一個點）。
　　其結果和原始訊號圖畫在一起如下

管路中之損失除管路本身外，其出口、入口及管徑擴張之損失均與動能有關，即V²/(2g)之係數。其程式內容如下

% pipefittings.m
% Local losses at pipe fittings
%
while 1
    k=menu('請輸入接合件之型式','驟然放大','驟然逸出','擴 大 管');
    Q=input('請輸入管路流量(m^3/s)[10]: ','s');
    if isempty(Q), Q='10';end
    D1=input('請輸入入口管徑(m)[1]: ','s');
    if isempty(D1), D1='1';end
    d1=str2num(D1);
    ke=str2num(Q)./(pi*d1.^2/4)/(2*9.81);
    switch k
        case 1
            HL=ke;
        case 2
            HL=ke/2;
        case 3
           D2=input('請輸入出口管徑(m)[2]: ','s');
           if isempty(D2), D2='2';end
           HL=(1-(d1./str2num(D2)).^2).*ke;
    end
    disp(['此接頭損失為 ' num2str(HL) ' (m).'])
    A=input('是否繼續執行(Y/N)[Y]:','s');
    if isempty(A), A='Y';end
    if upper(A)~='Y', break;end
    disp('請由選單選擇轉換單位。。。')
end

管路配件損頭

經過認養後

% pipefittings.m
% Local losses at pipe fittings
%
while 1
    k=menu('請輸入接合件之型式','驟然放大','驟然逸出','擴 大 管');
    switch k
        case 1
            prompt = {'請輸入管路流量(m^3/s)[10]','請輸入入口管徑(m)[1]:'};
            dlg_title = '您所選擇的是驟然放大';
            num_lines = 1;
            def = {'10','1'};
            answer = inputdlg(prompt,dlg_title,num_lines,def);
            Q=answer{1,1};d1=answer{2,1};
            ke=Q./(pi*d1.^2/4)/(2*9.81);
            HL=ke;
        case 2
            prompt = {'請輸入管路流量(m^3/s)[10]','請輸入入口管徑(m)[1]:'};
            dlg_title = '您所選擇的是驟然逸出';
            num_lines = 1;
            def = {'10','1'};
            answer = inputdlg(prompt,dlg_title,num_lines,def);
            Q=answer{1,1};d1=answer{2,1};
            ke=Q./(pi*d1.^2/4)/(2*9.81);
            HL=ke/2;
        case 3
            prompt = {'請輸入管路流量(m^3/s)[10]','請輸入入口管徑(m)[1]:','請輸入出口管徑(m)[2]:'};
            dlg_title = '您所選擇的是擴大管';
            num_lines = 1;
            def = {'10','1','2'};
            answer = inputdlg(prompt,dlg_title,num_lines,def);
            Q=answer{1,1};d1=answer{2,1};d2=answer{3,1};
            ke=Q./(pi*d1.^2/4)/(2*9.81);
            HL=(1-(d1./d2).^2).*ke;
    end
    figure('Position',[400 400 400 300],'Name','結果')
    uicontrol('style','text','position',[160 200 80 80],'string','此接頭損失為')
    uicontrol('style','text','position',[100 90 200 100],'string',HL)
    result=questdlg('是否繼續執行？','程式執行詢問','Yes','No','No');
    if strcmp(result,'No');
        break
    end
end

或是點此網頁

plot作圖時
如果新的plot圖時，會將前一次的洗掉，做出第二次的繪圖，舉例如下
第一次plot作圖

如果這時候在使用plot繪製另一個圖形的話，前一次的作圖就會消失
第二次plot作圖

若是要兩個圖顯示在一個figure的話，則必須用hold on指令
以hold on將兩張圖合併

如果要繪製三度空間的圖形的話，則必須要使用plot3的指令
plot3作圖

line作圖
若是以line的指令作圖時，連續使用line繪製不同函數時，則可以不必使用hold on的指令，即可繪製在同一個figure
以line繪製兩個函數

雖然說line可以給予三個變數繪製三度空間的圖，但是實際上使用後會發現繪製仍為2D平面的圖，並且是前面兩個變數的圖型
plot繪製三度空間圖

此時如果要會製成三度空間圖的話，則可以用plot3先繪製一點或是一個三度空間圖形
以plot3繪製三度空間點

最後在重新輸入line的函數，即可再三度空間繪圖
line指令三度空間繪圖

第一題
>> p=[133 0 122 0 0 1] % p=133x^5+122x^3+1

p =

133 0 122 0 0 1

>> q=[2 0 100 0 1] % q=2x^4+100x^2+1

q =

2 0 100 0 1

>> M=conv(p,q)

M =

Columns 1 through 9

266 0 13544 0 12333 2 122 100 0

Column 10

1

>> % p, q 兩多項式的乘積為 266x^9+13544x^7+12333x^5+2x^4+122x^3+100x^2+1

執行結果

第二題
>> p=[100 23 1 45];x=magic(5)

x =

17 24 1 8 15
23 5 7 14 16
4 6 13 20 22
10 12 19 21 3
11 18 25 2 9

以magic(5),x=每個元素帶入f(x)的答案為

>> f=polyval(p,x)

f =

498009 1395717 169 52725 342735
1228935 13125 35479 278967 415549
6817 22479 223645 809265 1075999
102355 176169 694267 936309 2955
135939 590715 1576945 939 74817

執行結果

第三題
>> p=[133 0 122 0 0 1] % p=133x^5+122x^3+1

p =

133 0 122 0 0 1

>> q=[2 0 100 0 1] % q=2x^4+100x^2+1

q =

2 0 100 0 1

>> [s,r]=deconv(p,q)

s =

66.5000 0


r =

1.0e+003 *

0 0 -6.5280 0 -0.0665 0.0010

>> % 商數s=66.5x 餘數r=e^3*(-6.528)x^3+e^3*(-0.0665)x+e^3*(0.001)

執行結果

第四題
>> p=[133 0 122 0 0 1] % p=133x^5+122x^3+1

p =

133 0 122 0 0 1

>> roots(p)

ans =

-0.0045 + 0.9578i
-0.0045 - 0.9578i
0.1039 + 0.1744i
0.1039 - 0.1744i
-0.1988

>> q=[2 0 100 0 1] % q=2x^4+100x^2+1

q =

2 0 100 0 1

>> roots(q)

ans =

0 + 7.0704i
0 - 7.0704i
0 + 0.1000i
0 - 0.1000i

執行結果

第五題
>> x=[1:9]; y=[1210, 1866, 2301, 2564, 2724, 2881, 2879, 2915, 3010];
>> p=polyfit(x,y,3)

p =

6.3047 -134.4603 994.3540 350.9127

>> f=@(a) (6.3047*a^3-134.4603*a^2+994.354*a+350.9127)

f =

@(a) (6.3047*a^3-134.4603*a^2+994.354*a+350.9127)

>> fplot(f,x)
>> hold on
>> plot(x,y)
>>

此題執行後，其polyfit的結果似乎只有x=1~2的時候比較符合，因為不管我限制多少範圍
做出來的圖ｘ座標都是從１～２
然後將其和原圖放在一起
執行結果

第一題
利用struct的方式撰寫

% student為整個struct的名稱   % 利用while的方式讓整個程式可以執行
student=struct('Name',0,'Age',0,'Sex',0,'Email',0);
ans='Y';
h=input('請問是否執行此資料程式？(Y/N) ', 's');
j=strcmp(upper(h),ans);
while j==1
    if student(1).Name==0                   %先判斷是否已經有存在資料庫
    i=0;
    else i=length(student);
    end
fprintf('目前資料庫有%d筆資料',i);
e=input('請問是否新增？（Y/N）','s');
f=strcmp(upper(e),ans);
if f==0                                     %判別起始回答時是否因為尚未建立而要去修改
    if i==0
        disp('目前資料庫裡面沒有資料,請新增')
    else
        g=input('請問要修改哪一筆資料？');
        i=g;
    end
    a=input('請輸入學生姓名（Name）　','s');
    b=input('請輸入學生年齡（Age）　','s');
    c=input('請輸入學生性別（Sex）　','s');
    d=input('請輸入學生電子信箱（Email）　','s');
    student(i)=struct('Name',{a},'Age',{b},'Sex',{c},'Email',{d});
elseif f==1
    i=i+1;
    a=input('請輸入學生姓名（Name）　','s');
    b=input('請輸入學生年齡（Age）　','s');
    c=input('請輸入學生性別（Sex）　','s');
    d=input('請輸入學生電子信箱（Email）　','s');
    student(i)=struct('Name',{a},'Age',{b},'Sex',{c},'Email',{d});
end
h=input('請問是否繼續執行此程式？(Y/N) ', 's');
j=strcmp(upper(h),ans);                    %決定是否要在繼續增加或修改資料
end 

執行結果1
執行結果2

第二題
利用Menu和Switch的方式撰寫

clear all
close all
k=menu('請選擇你喜歡的品牌','Apple','Microsoft','IBM','Acer','Asus');
switch k
    case 1
        disp('活著其實很好,在吃一口蘋果')
    case 2
        disp('準備迎接windows vista')
    case 3
        disp('無敵黑金剛,金剛不壞之身')
    case 4
        disp('為王建名加油')
    case 5
        disp('華碩品質,以卵擊石')
end

程式執行圖
執行結果圖

第四題
原始資料檔案
在command window下利用xlsread和xlswrite的方式修改資料
>> [num text raw]=xlsread('D:/grade.xls');
>> raw(2,:)

ans =

'Eric' [1] [180] [110]

>> raw(2,:)={'Eric',1,178,100}

raw =

'Name' 'Sex' 'height' 'weight'
'Eric' [ 1] [ 178] [ 100]
'Jane' [ 0] [ 160] [ 50]
'Jolin' [ 0] [ 158] [ 45]
'David' [ 1] [ 175] [ 54]

>> xlswrite('D:\grade.xls',raw)
程式執行圖
執行結果圖
確實在exel裡面修改了資料

第一題
自由落體：0~30s的時間內位置的高度的變化
>> t=(0:0.2:30);
>> height=freebody(t);
>> plot(t,height)
>> xlabel('time(s)');ylabel('height(cm)');
>> height
結果圖1 結果圖2

第二題
以半徑３和５並以圓心(0,0) and (10, 0) cm作圖
function drawcircles(r) % 主函數
% Draw circles with radii of r.
global MAXR i;
n=length(r);
MAXR=max(r);
hold on;
for i=1:n,
[X,Y]=circ(r(i));
plot(X,Y);
end
hold off;
axis equal;
end

function [x,y]=circ(rr) %次函數
% Calculate the points of a circle.
[xx,yy]=randxy;
theta=linspace(0,2*pi,60);
x=xx+rr*cos(theta);
y=yy+rr*sin(theta);
end

function [xx,yy]=randxy %次函數
global MAXR i
% locate the position of the center.
k=[0 10];l=[0 0];
xx=k(i);
yy=l(i);
end
作圖

第三題
以=3, b=5, 並且給予 y=[12 15 20] 然後以匿名函數作x²/a²+y²/b²
最後作圖在x=[3 4 5]的範圍下（應該是這樣吧，題意不是很懂）
a=3;b=5;
y=[12 15 20];
fplot(@(x) ((x.^2)./a^2+(y.^2)./b^2),[3 5])
結果圖

第四題
敘述檔為
Do=input('請輸入圓柱體的外半徑? ');
Di=input('請輸入圓柱體的內半徑? ');
height=input('請輸入圓柱體的高? ');
[a,b,c]=cylin(Do,Di,height);
fprintf('圓柱體的體積為%f,上表面積%f,側面面積%f\n',a,b,c)

呼叫function cylin
% function cylin %
function [volume,upper_face,side]=cylin(Do,Di,height);
volume=(Do./2-Di./2).^2*pi.*height; %體積
upper_face=pi*(Do./2-Di./2).^2; %上表面積
side=2*pi*(Do./2-Di./2).*height; %側面面積

結果在這

第五題
將程式改成呼叫外部的函數，並且設立global variable，即可將polyx移到原本的函數外面
而非原先的巢狀函數

function [rr_array]=nest_fun1(x,a)
%function to find sets of polynormials.
% a: set of constants, [A B C]
% x: variables in array
% Example: rr=nest_fun(2:10,[1 2 4;2 4 8])
n=size(a);
global A B C x
for i=1:n
  A=a(i,1);B=a(i,2);C=a(i,3);
  rr_array{1,i}=['A=',num2str(A),', B=',...
      num2str(B),', C=',num2str(C)];
  rr_array{2,i}=polyx(x);
end
end

function [r]=polyx(xx)
    global A B C x
    r=A.*x.^2 + B.*x +C;
end

結果在這

第六題

function pillar(Do,Di,height)
% Find the volume of a hollow pillar
% ro,ri:outside & inside diameters
if nargin<1 
   height=1; 
   Di=0;
   Do=1;
   volume=abs(Do.^2-Di.^2).*height*pi/4;
   fprintf('體積為%f（預設高為%d內面積%d外面積%d)\n',volume,height,Di,Do);
end
if nargin<2 & nargin>0
    height=1;
    Di=0;
    volume=abs(Do.^2-Di.^2).*height*pi/4;
    fprintf('體積為%f（預設高為%d內面積%d)\n',volume,height,Di);
end
if nargin<3 & nargin>1
    height=1;
    volume=abs(Do.^2-Di.^2).*height*pi/4;
    fprintf('體積為%f（預設高為%d)\n',volume,height);
end
if nargin==3
    volume=abs(Do.^2-Di.^2).*height*pi/4;
    fprintf('體積為%f\n',volume);
end

結果在這

　　話說現在的labview功能越來越強大了，
如果有看過Matlab簡易濾波器的話，看到這篇
之後應該會很想扁人。

　　看不清處嗎？那就按這邊下載


　　沒錯，就是這樣，開檔，濾波，畫圖。

　　首先我必須先聲明的是，我對於濾波的真正原理是
一知半解，但是他在我們物理治療中的訊號分析裡面，
卻是扮演著重要的角色。

　　什麼事butterworth的濾波，說真的我真的不知道，
反正數位濾波器就以這樣設計就好。

　　看到這裡如果你非常想打我的話，元智大學最佳化
實驗室有份文件大家可以參考看看。按我看看

　　言歸正傳，反正濾波器有高通，低通，帶通，帶拒
四大類，於是乎我利用了matlab裡面的butter函數，以
及filter設計了這個含有四種功能的濾波器function，
內容很簡單，只不過利用了switch的方式選擇。
點我看看這個簡易的濾波功能

　　如果有興趣試試看的話，這裡提供了一筆資料可以
試試看請以右鍵點選下載。

　　這分資料裡面，第一列跟最後一列分別都是force的
資料，中間幾個才是肌電圖，其取樣頻率為1000/s。

作業5.1 吳堉光
第一題：
利用陣列轉換函數（cat)依所指定的方向合併（concatenates）陣列A,B
程式碼：
A=magic(3);B=rand(3,3)*10;
cat(3,A,B)
或是直接
cat(3,magic(3),rand(3,3)*10)
執行結果

第二題：
將兩個陣列以水平或是垂直的方式合併，或是利用不同page的方式合併
水平合併的方式可以有兩種
A=[45 89 99; 12 34 55];B=[15 25 45; 65 50 30];
[A,B]
cat(2,A,B)
即可得到水平合併，點我看結果

垂直合併的方式也有兩種
A=[45 89 99; 12 34 55];B=[15 25 45; 65 50 30];
[A;B]
cat(1,A,B)
即可得到垂直合併，點我看結果

以不同page儲存的方式則要利用cat函數
A=[45 89 99; 12 34 55];B=[15 25 45; 65 50 30];
cat(3,A,B)
即可得到以不同page儲存，點我看結果

第三題
將A,B兩陣列以指定的位置放入page中
方式一
M(:,:,1)=A; M(:,:,2)=B; M(:,:,3)=A; M(:,:,4)=B;
執行結果
方式二
如果今天有大量的陣列放入而且是有規則性的放入的話，可以以for迴圈來執行
A=[45 89 99; 12 34 55];B=[15 25 45; 65 50 30];
for i=1:2
M(:,:,2*i-1)=A;
M(:,:,2*i)=B;
end

第四題
依題意看起來應該是要形成一個3 x 2的cell array
所以可以直接利用cell array的{}的方式
c={'Eric' [90 100]; 'Peter' [35 60]; 'Jan' [77 67]}
所形成的結果以cellplot的方式來表示
結果在這裡

第五題
將資料建構成一個struct的陣列
程式如下
patient(1).name='Philip';
patient(2).name='Peter';
patient(3).name='Roberts';
patient(4).name='Roe';
patient(1).Age=35;
patient(2).Age=45;
patient(3).Age=55;
patient(4).Age=60;
patient(1).salary=50000;
patient(2).salary=40000;
patient(3).salary=80000;
patient(4).salary=120000;
即可以建構出struct的陣列
不過這樣一個一個key實在挺累人的

第一題：
% 求兩個向量的外積的單位向量
% A=3i +5j +10k B=5i -6j +2k
A=[3 5 10];B=[5 -6 2];
m=cross(A,B);
r=sqrt(m*m');
u=m/r

第二題
由於在計算平均值，中位數時其syntax中（A,dim）當dim不寫或是1時，所計算的為每一行的平均值及中位數，而當填入為2時則為計算列
Mean&Median

而std其syntax
flag=0時為sample的標準差，而flag=1為population
STD

第三題
根據老師上所講的尤拉公式，在command windows下證明
證明結果

第四題
計算複數所成的陣列的絕對值,夾角,實部,虛部
絕對值,夾角
實部,虛部

第一題：
個人認為減少迴圈加快執行的方式為
n=(1:50).^2;
reshape(n,5,10)
這樣可以得到相同的答案

第二題
A=magic(5);
for n=1:5
for m=1:5
if A(n,m)>10
A(n,m)=NaN;
end
end
end
即利用迴圈一一的去檢查是否符合條件，符合的就替換成NaN

第三題
A=magic(5);
for n=1:5
for m=1:5
if isprime(A(n,m))
A(n,m)=0;
end
end
end
和第二題一樣，以迴圈的方式去檢測

第四題
將no_of_draw改成[100:100:1000]會發現到結果和no_of_draw的答案差不多，也就是只會執行100次，我們可以利用
a=100;b=(100:100:1000);
c=(a==b)
c=[1 0 0 0 0 0 0 0 0 0]
我想應該是這樣，所以程式裡面的條件
while n<=no_of_draw
因為第一個元素不符條件了，所以程式結束。
如果要一次全部呈現可以改成如下程式
function [B]=draw_number(no_of_draw)
% draw ball numbers within ndraw times
%
for i=1:length(no_of_draw)
C_i=zeros(1,5);
n=1;
while n<=no_of_draw(1,i)
ball=fix(rand*10);
if ball<2
C_i(1)=C_i(1)+1;
elseif ball<4,
C_i(2)=C_i(2)+1;
elseif ball<6,
C_i(3)=C_i(3)+1;
elseif ball<8,
C_i(4)=C_i(4)+1;
else
C_i(5)=C_i(5)+1;
end
n=n+1;
end
B(i,1:5)=C_i;
end
第五題
clear all
k=input('Please input a number to compare? (0-100)');
l=input('Please input the total number to do?');
a=round(k);
n=0;
num=0;
answer=0;
if a>=0 & a<=100
for i=1:l
while round(rand*100)~=a
n=n+1;
end
num=num+1/n;
end
answer=num/100;
fprintf('%d%s%d%s%f\n',l,'次抽中和',a,'相同的號碼機率為千分之',answer*1000);
else
fprintf('%s\n','你所輸入的號碼範圍號碼不對')
end
可以同時輸入號碼跟比對次數

　　在我們物理治療研究中，不管是研究肌械延遲（electrical mechanical delay），或是病人對於肌肉張力，運動員之力量發展，都需要靠肌電訊號來告訴我們，其中有一個最重要的部份就是，去知道該塊肌肉什麼時候開始動作。

下載檔案（點右鍵下載此檔案）
　　其中，前三條分別為我們下肢小腿後側肌群。最後一條為實驗時我用數位類比轉換卡發出的類比訊號給儀器接收，以作為同步訊號。不過目前先不將同步訊號考慮在內。

　　我們必須要去找到開始的時間點，此時就要利用到統計學的概念，兩個母群若是差兩個標準差的話，就可以稱作是兩個不同的母群（否決虛無假設）。
不過很不幸的是，通常訊號可能沒有那麼漂亮，或是雜訊很多，可能需要到三個甚至四個標準差以上。

點此下載看上圖
程式檔案

　　這個程式就是靠while+break所寫成的，當找到符合條件的點後，就結束迴圈，並且輸出結果。

　　之後都會上傳我的上課筆記，由於是邊上課邊用
writely寫的，所以並不會很有組織的一一紀錄上課
所教的東西，很多都是零碎的片段，或是一些小技巧
。

　　switch的用法，以及介紹如何在cell資料型態運用，特別是
利用page的方式，讓整個array多了三度的感覺

按我看上課筆記

　　主要以matlab中比較的關係句，以及如何在editor下除錯（筆記很少，看來我很不專心）

按我看上課筆記

　　對於我們大學四年都在醫學院的學生而言，姑且
不論工程，機電的問題，甚至連最基本微積分，聯立
方程式可能都不記得怎麼算，更不用說是寫程式了（
無法像PCman的作者如此厲害)。
　　
　　碩一上結束後的寒假，我毅然決然的決定來學寫
程式，到現在也已經有了半年的時間，剛開始從看似
容易的labview學起（因為看起來不用輸入程式語法）
，寒假大家都在出去玩，而我則一個人埋頭看書，自
己讀程式語言果然很痛苦，而我們醫學院這邊又沒有
提供這樣的資源讓我們可以去學習程式語言。
　　
　　即使是看了整個寒假，九九乘法表幾乎想了快一
個月才寫得出來，迴圈的用法到現在還是"哩哩辣辣"
，不過對於我而言，寫程式反而讓我回到過去高中的
時候，算數學物理時，可以忘卻時間，全心全意思考
一個問題，去嘗試，去失敗，去學習。

　　半年過去後的暑假，我拿起Matlab的書本，開始
接觸不敢碰觸的語法，直譯式的程式語言，可以快速
的知道自己錯在哪裡，慢慢的我將他利用在我的實驗
分析上面，Matlab和Labview相輔相成。

　　話說到總區修這堂課，當初心中有著千百個不安
（現在還是），對於我們這種整天都跟骨骼肌肉神經
為伍的人，工程的東西，程式的東西，都是讓我所恐
懼而不敢向前的。但是既然自己決定要學，就該不顧
一切的往前，希望到了這裡，能夠學會如何去思考程
式邏輯，藉由看大家寫程式的經驗，學習一些程式的
撰寫技巧，增進自己的實力。也希望各位先進，不吝
惜自己的實力。

　　個人願以野人獻曝的心情，且謙卑的向各位學習
。

　　

　　一提到MATLAB很多人可能就馬上想到張智星老師
，他的書雖然是入門款，不過剛開始看還真的有點小
難，不過市面上MATLAB的書琳瑯滿目，選擇對自己比
較容易理解的比較重要。

　　他的網頁裡面有許多東西，也有他上課的講義，
比較值得注意的音訊分析，裡面的一些概念以及濾波
的寫法，對於我們肌電訊號的處理我想應該是非常有
幫助。

　　這是我之前還蠻常瀏覽的網站，可以
參考許多高手對於一些問題的解決方式，
以及一些寫程式的技巧，不過因為我是初
學，很多討論的東西其實我都看不懂，不
然就是也不知道他們為何如此撰寫，不過
對於labview而言，是非長不錯並且可以參
考的地方。

這是馮老師上課內容的網站，介紹matlab的語
法與基本概念。

　　昨天詢問過老師，同意讓我將他的部落格放入
連結中，這樣可以直接連結到老師的matlab的部落
格。

　　由於東西都是放於網路免費空間上面，難保可
能會出現連結失效或是網路空間的意外，如有發生
連結失效的情形，請寫信給我：

我的信箱為happyeric1120@yahoo.com.tw

　　在EMG訊號中，不管是要作iEMG或是RMS，甚至式自己想要作
windows來分析，還是某一區域的最大最小值，都需要透過分割資
料應用，自己過去也想了一些方式，這份心得筆記是自己之前寫的
，大家不彷嘗試看看，相信一定有更好的方式，也歡迎大家提供。

點我看如何擷取某段資料

RMS常用在處理EMG的訊號上面，但是很少有去直接計算RMS的軟體，大多要靠自行寫程式，方法其實很多，在此只是給予一個我所寫的方式。

　　為了方便查詢以及分類，所以利用label的方式將文章一一分門。內容包含有

labview-訊號擷取
labview-訊號分析
matlab-訊號分析 等等

　　其他閒聊可以歸類在other。這樣應該會比較容易查找。