APPLICATIONS OF MATLAB ON ENGINEERING PROBLEMS

12/22/2006

肌電訊號（EMG）之頻域分析

肌電訊號分析可以分為兩個主要部份，一個即是時域
（time domain）另一則是頻域（frequency domain）。

這裡介紹物理治療會使用去計算頻率變化的方法，利
用傅立葉轉換的方式求得頻率和其強度的關係圖。

按我連結

12/15/2006

林冠宏
R95631022
之前以有人做回應
我想表中所對應的值,有時會不齊,所以讓可以找到一個表中對應較近的值,是否會比較好,就可以有所依據,去做選擇
及修改yes和no大小寫都轉為大寫
程式流程:
% gear.m
% Design of a gear
efc=0.86;g=9.81;% efficiecy of gear
gcof=[10 0.176 0.201
11 0.192 0.226
12 0.210 0.125
13 0.223 0.261
14 0.236 0.276
15 0.215 0.289
16 0.254 0.295
17 0.264 0.302
18 0.270 0.308
19 0.276 0.314
20 0.283 0.320
21 0.289 0.327
23 0.295 0.333
25 0.305 0.339
27 0.311 0.349
30 0.317 0.358
34 0.327 0.371
38 0.333 0.383
43 0.339 0.396
50 0.346 0.408
60 0.355 0.421
75 0.361 0.434
100 0.368 0.446
150 0.374 0.459
300 0.383 0.471
999 0.390 0.484];
gmax=length(gcof);
while 1

while 1
GT=input('請輸入齒輪之齒數(-)[預設值為30]:
','s'); %輸入齒輪之齒數

if isempty(GT), %檢查是否有輸入,沒有則令為30
GT='30';
end

gteeth=str2num(GT);
[N,I]=find(gcof(:,1)==gteeth);

if isempty(N) %找表中較近的值

while isempty(N)
gteeth=gteeth+1;
[N,I]=find(gcof(:,1)==gteeth);
end

disp([' 表中較近的值為 ' num2str(gteeth)]);
on = upper(input('是否使用(Y/N) [Y]
','s'));

if on ~= 'Y'
gteeth=0;
[N,I]=find(gcof(:,1)==gteeth);
end

end

if isempty(N), %檢查齒數是否在表中可以查到
disp('齒數不符合，請重新輸入！');
else
break;
end
end
while 1

GP=input('請選擇壓力角[14.5度(1)/20度(2)][預設值為1]:
');%選擇壓力角
if
isempty(GP),%檢查是否有輸入,沒有則令為1
GP=1;
end

if GP==1 |
GP==2,%如果是1或是就跳出
break
else

disp('請選１或２，重新輸入！');
end
end

GF=input('請輸入齒輪之作用力(kg)[預設值為200]:
','s');%輸入齒輪之作用力
if isempty(GF),
GF='200';end;gf=str2num(GF);%檢查是否有輸入,沒有則令為200

GD=input('請輸入齒輪節徑(mm)[預設值為200]:
','s');%輸入齒輪節徑
if isempty(GD),
GD='200';end;gd=str2num(GD);%檢查是否有輸入,沒有則令為200

GB=input('請輸入齒厚(mm)[預設值為100]: ','s');
%輸入齒厚
if isempty(GB),
GB='100';end;gb=str2num(GB);%檢查是否有輸入,沒有則令為100

RPM=input('請輸入齒輪迴轉速(mm)[預設值為1200]: ','s');
%輸入齒輪迴轉速
if isempty(RPM),
RPM='1200';end;rpm=str2num(RPM);
%檢查是否有輸入,沒有則令為1200

%計算抗彎應力

gv=rpm*pi.*gd/60/1000;
gfv=3./(3+gv);
gm=gd/gteeth;
gy=gcof(N,GP+1);

gst=gf./gfv./gb./gm./gy;

%輸出資料
disp('齒輪資料：')
disp(['齒數：'
num2str(gteeth) ' 節徑：' num2str(gd)...
'mm. 齒厚：'
num2str(gb) 'mm.'])
disp(['轉速：'
num2str(rpm) 'rpm. 作用力：' num2str(gf) ' KG.'])

disp(['齒輪之抗彎應力為: ' num2str(gst) ' kg/m^2.
<=====']);

%是否還要計算

A=upper(input('繼續執行嗎？(Y/N) [Y] ','s'));

if ~isempty(A) |
A~='Y'
break;
end
end

結果(1):
請輸入齒輪之齒數(-)[預設值為30]: 40
表中較近的值為 43
是否使用(Y/N) [Y] y
請選擇壓力角[14.5度(1)/20度(2)][預設值為1]:
請輸入齒輪之作用力(kg)[預設值為200]:
請輸入齒輪節徑(mm)[預設值為200]:
請輸入齒厚(mm)[預設值為100]:
請輸入齒輪迴轉速(mm)[預設值為1200]:
齒輪資料：
齒數：43 節徑：200mm. 齒厚：100mm.
轉速：1200rpm. 作用力：200 KG.
齒輪之抗彎應力為: 6.5817 kg/m^2. <=====
繼續執行嗎？(Y/N) [Y]
結果(2):
請輸入齒輪之齒數(-)[預設值為30]: 40
表中較近的值為 43
是否使用(Y/N) [Y] n
齒數不符合，請重新輸入！
請輸入齒輪之齒數(-)[預設值為30]: 30
請選擇壓力角[14.5度(1)/20度(2)][預設值為1]:
請輸入齒輪之作用力(kg)[預設值為200]:
請輸入齒輪節徑(mm)[預設值為200]:
請輸入齒厚(mm)[預設值為100]:
請輸入齒輪迴轉速(mm)[預設值為1200]:
齒輪資料：
齒數：30 節徑：200mm. 齒厚：100mm.
轉速：1200rpm. 作用力：200 KG.
齒輪之抗彎應力為: 4.9105 kg/m^2. <=====
繼續執行嗎？(Y/N) [Y] n

_______________________________________
YM - 離線訊息
就算你沒有上網，你的朋友仍可以留下訊息給你，當你上網時就能立即看到，任何說話都冇走失。
http://messenger.yahoo.com.hk

吳子青：資料的前處理

吳子青：資料的前處理

前言：
這週開始應該會進入統計的課程，我想對統計大家都有一定的認識，從分佈、主成份分析和變方分析等，都是常用的方法。而且就算所研究的主題和統計一點關係都沒有，在論文和期刊上也是會要求利用討統計方法驗證結果的可靠度。

其實撇開統計上的理論方法，統計上最重要的是資料的品質，有好的資料品質分析出來的假設檢定才會具有意義。所以在進入統計前我們先要思考的問題是，我的資料夠好嗎？要如何改善？

我的問題是：
如果遇到資料有遺漏的處理方法

內容與方法
我寫了兩個小程式解決遇到資料有遺漏的處理方法，因為要分析的資料通常是千筆以上，不太可能手動自己算或用Excel慢慢拉，在分析數據上遇到資料有遺漏是常常有的狀況，但程式也只會回給你一個Error ，甚麼都不會和你說。

假設今天所要處理的資料有n 筆，每一筆資料共有p個變數，所以資料矩陣就是一個 n x p的矩陣。

處理方法一：
將有遺漏資料的那一筆刪除，因為不同資料長度或是不同的樣本作統計分析是無意義的，例如想知道一天之內日溫的變化，資料取禮拜一到禮拜日 (n=7)，變數為每小時的日溫(p=12)，假如溫度計禮拜三的中午忘記紀錄，這樣把禮拜三的11筆資料和禮拜四的12筆資料作比較是無意義的，所以處理方法一就是把禮拜三刪除。

處理方法二：
如果資料本身並不夠多，如果用處理方法一去處理就如同落井下石，所以在熟悉資料的特定或是合理的假設下(重要)，可以對資料進行補遺的動作，在這裡使用的是最基本的內插法，前提是資料變動性不大或者變動是線性的。
以上面的例子解釋就是把禮拜三11點的日溫和13點的日溫作內插，估計12點的日溫。

程式流程：
在這裡要處理的資料設定為從Excel載入到Matlab的數據，所以當資料有遺漏處會顯示 “ NaN “ ，所以在這裡以指令 “isnan”來找出數據遺漏的地方，在進行處理。
Isnan(A) 的結果為和原矩陣同樣維度的矩陣，而A裡面NaN的元素對應到a是1，而其他的地方為0

處理方法一：
Step 1：資料矩陣為 A，a= isnan(A) 為判斷矩陣
Step2 : 從第一筆資料開始跑迴圈，然後慢慢的累加結果，而當遇到a矩陣裏出現1時直接跳過
Step3: 所以最後加總的結果就是刪除缺失資料的處理後數據

程式內容：


function [out]=de(A);
a=isnan(A(:,5));
[mm,nn]=size(A);
B=zeros(1,nn);
for i=1:mm
    if a(i)==0
        B=[B;A(i,:)];
    end
end
out=B(2:end,:);

處理方法二
Step 1：資料矩陣為 A，a= isnan(A) 為判斷矩陣
Step2：搜尋每一筆資料中從0到1的位置，代表開始數具有缺失，之後利用while指令找出1總共延伸有幾個，代表中間有幾個數據要作內插。
Step3：將中間需補遺的地方用前一個和後一個資料進行內插的計算在矩陣B
Step4：跳過已搜尋的位置對全部的資料筆數進行補遺的動作
Step5：將原矩陣有資料的地方填補在矩陣B上面，完成處理
Step6：決定外插處理的方式

程式內容：


function [out]=inc(N12);

[mm,nn]=size(N12);
a=isnan(N12);
out=zeros(mm,nn);

for k=5:nn
    switch sum(a(:,k))
        case mm
        otherwise        
m=1;
s=1;
while a(s,k)==1
    s=s+1;
end

for i=s:mm
    if i>m
        if a(i,k)==1
            n=i;
            w=n+1;
                while a(w,k)==1
                    if w"<"mm
                         w=w+1;   
                    else
                        m=w;
                        A(n:m-1,1)=0; %(data尾內插修正)
                        break;
                    end
                end
            m=w;
            if m~=mm
                for q=n:m-1
                    out(q,k)=N12(n-1,k)+(N12(m,k)-N12(n-1,k))/(m-n+1)*(q-n+1);
                end
            end
        end
    end
end
clear i w;
    end
end
for i=1:mm
    for j=1:nn
        if a(i,j)==0
            out(i,j)=out(i,j)+N12(i,j);
        end
    end
end

實際案例：

圖一原始資料為Excel格式中間有很多遺漏的資料
https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiBq50IiJDdqRdSucYkKME6gY8MpoJ1wNVsi0TFdZrwYi5ZHYwkBjRYU5f_gM99GRgPcL11kQyKeoI_m4xfOXbVVd6GiWM-944mjs4V5WORr4iYr8wxl45RwoxbUEPSr9uh76Ox9pvhwzY/s1600-h/11.JPG

圖二用Matlab把資料載入遺漏的資料會顯示NaN
https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjzzVZR6OmV7OJRGy_4qHHuswfY1JJ9K_uQ00a_XCaN7oqiAvrIOaNVfHgCnr1OczAdaJvO6p4XJbgzn7QbhvrfwEWE6wqJYuIEeb66QXN8GhkZd_9bREe84mFoxeyaPY63QHX0oSdHxSo/s1600-h/12.JPG

圖三用處理方法一處理 (合併資料) 注意看年份的排序而且從第五行開始處理
https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgwesaPRU24mIoEXEKMeNQamu9Uj4vpHnAzh1DAAOIzTUOHxn67lXQpXshQ0yrAYHQ56nig1ZfKZfFjxBy7fUiwLDSIyWlKMjNjLQbVOhLkB2DZWXuKXnaf6S82X2SEUFUTa-tuG97PD9M/s1600-h/13.JPG

圖四用處理方法二處理 (內插補遺) 從第五行開始處理
https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgZhCmxWC0zKIs_vfZHDdYuxyz72N1FyHtRuWTpB3C1AyWuGQpRUd5om6fDCdWuOfSrHIjMqgwvoTi50JvuRmZ5MRo0nKGV4sX3Ja9S0ZY_YoTQ3ioxmp9uO5fH9F5Vifhws3FxeanpZM8/s1600-h/14.JPG

林書如：人體受外力牽拉後引發之姿勢反應_肌肉活化情形

林翊展--水分蒸發潛熱

在濕氣特性的章節中,根據老師所提供的公式(Brooker, 1967)

1.當 273.16K ≦ 亁球溫度(Tdb) ≦338.72K 時

hfg = 2,502,535.259 – 2,385.76424 (Tdb- 273.16)

2.當 338.72K ≦ 亁球溫度(Tdb) ≦533.16K 時

hfg = (7,329,155,978,000 – 15,995,964.08 Tdb2)^1/2

3.當 255.38K ≦ 亁球溫度(Tdb) ≦273.16K 時

hsg =-2,839,683.144 – 212.56384(Tdb-255.38)

但是根據老師所給的程式,似乎無法跑出結果,所以另行撰寫了一個程式,可以直接輸入攝氏溫度,得到該溫度下相對應的水分蒸發潛熱.

程式碼如下:

####################################################


function [h]=hfg(Tc)
% function [h]=hfg(Tc)
% find the latent heat, kJ/kg, Tc dry-bulb temp in C, in column matrix
T=Tc+273.15;
if T>=273.16 & T<338.72
   h=2502535.259-2385.76424.*(T-273.16);
elseif T>=338.72 & T<533.16            
   h=(7329155978000-15995964.08.*T.^2).^1/2;       
elseif T>=255.38 & T<273.16
   h=-2839683.144-212.56384.*(T-255.38);
else
   disp('overrange')
end

####################################################

此functuion命名為"hfg"
執行結果如下:

指令執行結果一

如果以矩陣輸入,結果如下

指令執行結果二

再加入兩條指令,畫成曲線圖(攝氏1~10度)

>>plot(hfg([1:10]'));
>>xlabel('Temperatures, C');
>>ylabel('Evaporation heat, J/kg');

曲線圖

利用moving window的方式分析EMG

　　分析肌電訊號的時候，我們會以一個類似窗口
的方式去對一個時間片段作處理，然後這個窗口會
隨著訊號的時間點一直掃過，所經之處及分析該窗
口內的訊號，並且經由指定的運算方式得到答案。

　　裡面的運算，常見的有整流（負翻正）平均，
方均根，或是一階的濾波，還有就是積分。

一般窗口的作法是用for loop然後可以指定要
重複多少點。以下便以簡單的程式碼方式完成。

按此連結到此文章

12/14/2006

詹弘彥：齒輪設計

正齒輪簡介：
齒輪常裝於軸上以傳動動力，它的種類包括正齒輪（spur gear）、螺旋齒輪、蝸輪及斜齒輪等，正齒輪為最便宜的一種，它的齒是與軸心平行，兩正齒輪相互作用時，兩齒輪之軸在同一平面且平行。

題旨：設計一齒輪是否適用？？
抗彎強度定義：
齒輪之容許負荷，為互相咬合而傳達動力之齒輪間，基於齒根部所受之抗彎應力產生之咬合節圓上之容許圓周力。


基本參數:
齒輪之抗彎應力gst=gf./gfv./gb./gm./gy 
 gf：齒輪之作用力 
 gfv=3./(3+gv)：圓周速度所生齒輪材料之使用應力(kg/cm2) 低速度(V=10m/s以下) 
  齒輪迴轉速rpm
  gv=rpm*pi.*gd/60/1000
 gb：齒厚 
 gm=gd/gteeth：模數 
  齒輪節徑gd
  齒輪之齒數gteeth
 gy：路易斯形因子 
  gcof=[齒輪數 壓力角14.5度之路易斯形因子 壓力角20度之路易斯形因子]
  壓力角GP

=======================程式內容=====================


% gear.m
% Design of a gear
efc=0.86;g=9.81;% efficiency of gear
gcof=[10 0.176 0.201;11 0.192 0.226;12 0.210 0.125;13 0.223 0.261;...
   14 0.236 0.276;15 0.215 0.289;16 0.254 0.295;17 0.264 0.302;...
   18 0.270 0.308;19 0.276 0.314;20 0.283 0.320;21 0.289 0.327;...
   23 0.295 0.333;25 0.305 0.339;27 0.311 0.349;30 0.317 0.358;...
   34 0.327 0.371;38 0.333 0.383;43 0.339 0.396;50 0.346 0.408;...
   60 0.355 0.421;75 0.361 0.434;100 0.368 0.446;150 0.374 0.459;...
   300 0.383 0.471;999 0.390 0.484];
gmax=length(gcof);
while 1
  while 1
      GT=input('請輸入齒輪之齒數(-)[30]: ','s');
      if isempty(GT), GT='30';end;gteeth=str2num(GT);
      [N,I]=find(gcof(:,1)==gteeth);
      if isempty(N),
          disp('齒數不符合，請重新輸入！');
      else
          break;
      end
  end
  while 1
     GP=input('請選擇壓力角[14.5度(1)/20度(2)]: ');
     if isempty(GP), GP=1;end
     if GP==1 ¦ GP==2,
         break
     else
          disp('請選１或２，重新輸入！');
     end
  end 
  GF=input('請輸入齒輪之作用力(kg)[200]: ','s');
  if isempty(GF), GF='200';end;gf=str2num(GF);
  GD=input('請輸入齒輪節徑(mm)[200]: ','s');
  if isempty(GD), GD='200';end;gd=str2num(GD);
  GB=input('請輸入齒厚(mm)[100]: ','s');
  if isempty(GB), GB='100';end;gb=str2num(GB);
  RPM=input('請輸入齒輪迴轉速(rpm)[1200]: ','s');
  if isempty(RPM), RPM='1200';end;rpm=str2num(RPM);
  gv=rpm*pi.*gd/60/1000;
  gfv=3./(3+gv);
  gm=gd/gteeth;
  gy=gcof(N,GP+1);
  gst=gf./gfv./gb./gm./gy;
  disp('齒輪資料：')
  disp(['齒數：' num2str(gteeth) ' 節徑：' num2str(gd)...
      'mm.  齒厚：'  num2str(gb) 'mm.'])
  disp(['轉速：' num2str(rpm) 'rpm.  作用力：' num2str(gf) ' kg.'])
  disp(['齒輪之抗彎應力為: ' num2str(gst) ' kg/m^2. <=====']);
  A=input('繼續執行嗎？(Y/N) [Y] ','s');
  if ~isempty(A) ¦ A~='Y', break;end
end

12/12/2006

吳子青：熱流分析之通用情況推廣

前言：
上禮拜上課老師講解了10-2 範例五：簡易之熱流分析
並探討了2x2 和 3x3 的範例
我就想寫一個n x m 維通用Case的程式
因為探討在Steady state穩態熱流情況基本原理並不難

除了把探討的維度改成n x m 外，在這裡也改成可以自由選擇熱流輸入
和輸出的位置，也可以多點的輸入輸出。

注意的地方是Control Volume 控制體積仍然維持正方形，原因在上課講解的時候有解釋過了，而且這裡探討的是Steady state。

分析與流程：
這個問題是示範反矩陣解聯立的方法
所以重要的就是在於係數矩陣 (也就是聯立的係數) 怎麼寫
而C矩陣只和輸入輸出熱流項有關

首先介紹函式的用法
function [T]=heat_platen(n,m,ti,to,input,output)
輸入項 n 和m 代表熱流場的維度
Ti: 輸入溫度
To: 輸出溫度
Input：熱流輸入的位置座標，以矩陣表示並可以多選。
Output：熱流輸出的位置座標

輸出項： T 溫度矩陣
並在輸入項加入了預設值的設計
也就是如果輸入沒有填熱流流出和流入的位置
程式將會預設成從( 1, 1 ) 傳入而傳到 ( n , m ) 熱流流出，也就是範例五的樣式

接下來簡略的說明想法
把2x2的係數矩陣和3x3的係數矩陣列出來

2x2 case :

3 -1 -1 0
-1 2 0 -1
-1 0 2 -1
-1 -1 0 3

3x3 case

3 -1 0 -1 0 0 0 0 0

-1 3 0 -1 -1 0 0 0 0

0 -1 2 0 0 -1 0 0 0

-1 0 0 3 -1 0 -1 0 0

0 -1 0 -1 4 -1 0 -1 0

0 0 -1 0 -1 3 0 0 -1

0 0 0 -1 0 0 2 -1 0

0 0 0 0 -1 0 -1 3 -1

0 0 0 0 0 -1 0 -1 3

可以觀察到有一定的規則性
在這裡分成矩陣的對角線元素的探討 A 和其餘部份的探討 B
然後最後再組合起來，最後在加上輸入輸出項

(1) 對角線元素為相鄰的區塊數目，如果是在四個角落為2，在矩形的四個邊
為3，而其他的區域為4，由這樣的規則成矩陣A
(程式 line 27 到 line 46 )

(2) 如果探討的控制體積是第k塊且是位在中間區域 (即不是四角及邊線)，
則在k+1、k-1、k+m、k-m的位置都要填 -1，而在邊線部份或四個角就只要填部份三個或兩的位置為 -1 ，並要特別注意i= 1 or n 或 j= 1 or m 的情況
由這樣的規則形成矩陣B
(程式 line 51 到 line 77)

(3) 加上輸入和輸出的部份，在輸入和輸出位置對應的對角線位置分別 +1 ，而在矩陣C 填上輸入溫度Ti 和輸出溫度To

(4) 最後利用反矩陣計算聯立方程式的解 T= inv(A)*C T為溫度矩陣
也是最終 Steady state 的溫度分佈，並繪製成3D圖

程式內容


function [T]=heat_platen(n,m,ti,to,input,output)
% Prog calculating heat transfer through a plate by 3x3.
% Inputs:
%   ti,to : input & output temperature, C
%    n,m  : n x m mesh (Heat flow region)
%   input : Heat flow input location [x,y] x=1~n,y=1~m
%   output: Heat flow output location [x,y] x=1~n,y=1~m
% Outputs:
%   temp:temperatures at each layer,C
% Example:
%  T=heat_platen(20,-10,20,20,[1 1],[20 20])

if nargin==5
    outx=n;outy=m;
    inx=input(:,1);
    iny=input(:,2);
elseif nargin==4
    outx=n;outy=m;
    inx=1;iny=1;
else
inx=input(:,1);
iny=input(:,2);
outx=output(:,1);
outy=output(:,2);
end

A=zeros(n,m);
B=zeros(n,m);
C=zeros(n*m,1);
C((inx-1)*m+iny,1)=ti;
C((outx-1)*m+outy,1)=to;

    for i=1:n
       for j=1:m           
          if i==1|i==n                 
                 A(i,j)=3;
             elseif j==1j==m
                 A(i,j)=3;
          else
                 A(i,j)=4;            
          end
          if (i==1i==n)&(j==1j==m)
                A(i,j)=2; 
          end
       end
    end
    B(inx,iny)=1;
    B(outx,outy)=1;
    A=A+B;
  
    TT=zeros(m*n,m*n);
    for k=1:n*m
        for i=1:n
            for j=1:m
                if k==(i-1)*m+j
                    TT(k,k)=A(i,j);
                end
                if i==1
                    TT((i-1)*m+j,i*m+j)=-1;
                elseif i==n
                    TT((i-1)*m+j,(i-2)*m+j)=-1;
                else
                    TT((i-1)*m+j,i*m+j)=-1;
                    TT((i-1)*m+j,(i-2)*m+j)=-1;
                end
                 if j==1
                    TT((i-1)*m+j,(i-1)*m+j+1)=-1;
                elseif j==m
                    TT((i-1)*m+j,(i-1)*m+j-1)=-1;
                 else
                     TT((i-1)*m+j,(i-1)*m+j+1)=-1;
                     TT((i-1)*m+j,(i-1)*m+j-1)=-1;
                 end
               
            end
        end
    end
    T=reshape(inv(TT)*C,n,m);
    mesh (T); figure(gcf)

結果

Case 1 改變熱流場Size
圖一5x5 Case
圖二10x10 Case
圖三50x50 Case

(輸入溫度 Ti = 20度在 (1,1)
輸出溫度 To = -10 度在 (n,m))

Case2 熱流入流出位置的改變
圖四 Input = [ 1 1] Output = [1 20]
圖五 Input = [ 1 1] Output = [10 10]
圖六 Input = [ 10 10] Output = [1 1]

(輸入溫度 Ti = 20度輸出溫度 To = -10 度
Size= 20 x 20 )

Case3 多個熱流出入口的改變
圖七 Input = [ 1 1; 20 1] Output = [1 20]
圖八 Input = [ 1 1; 20 20] Output = [1 20 ; 20 1]
圖九 Input = [ 1 1 ; 1 20 ; 20 1 ; 20 20] Output = [10 10]
圖十 Input = [ 10 10] Output = [ 1 1 ; 1 20 ; 20 1 ; 20 20]

12/08/2006

林書如：學生期中考成績分佈情形

因為脫離物理這門學科太久了
所以老師網頁上的題目都讓我覺得很陌生....
因此，我是以網頁內跟統計較相關的程式"度數分配程式"做了一些修正以及增加一些小指令，以敘述統計方式呈現學生期中考成績分佈情形。

12/01/2006

吳子青: 影像處理及後製

http://r95622024.blogspot.com/2006/11/blog-post_2921.htm

前言:
常常在攝影的人都知道
用單眼數位拍出來的照片通常需要 “後製”
也就是對照片作修改

因為常常用到像photoimpact或photoshop等的編輯軟體
突然想到可不可以利用
Matlab的強大功能作類似的處理

經過摸索之後找到可以把照片載入matlab的方法
並進行幾項簡單的修改

介紹:

在這裡我選擇了幾項簡單的功能並製作了一個選單
基本的功能包括
1. 觀看原圖 View photo
2. 顏色濾鏡功能 RGB filter
3. 反白 Counter image
4. 彩色照片轉黑白照 Color ro B/W

功能說明：
1. View photo
觀看指定的照片並有放大縮小的功能

2. RGB filter
就像photoshop的濾鏡功能，分成R、G、B 濾鏡
點選之後會輸入三個數值，r、g、b 分別是0~1的數值
1 代表顏色原始呈現 0 代表完全被遮住。

3. Counter image
也就是反白功能，能將明暗對調

4. Color ro B/W
把彩色照片轉換成黑白照片

使用方法：
把想處理的照片 xxx.xxx (如 w.jpg) 放在matlab的 Current Directory。
輸入指令 : Photo
打入 xxx. xx (記得打副檔名)
就進入選單模式

範例：
選這張照片進行處理 “ w.jpg “

5. (1) View photo
會顯示原照片

 http://photos1.blogger.com/x/blogger2/2907/859880749988154/1600/963015/w.jpg

6. (2) RGB filter

R=1
G=0
B=0

http://photos1.blogger.com/x/blogger2/2907/859880749988154/1600/36517/r.jpg

R=0
G=1
B=0

http://photos1.blogger.com/x/blogger2/2907/859880749988154/1600/388570/g.jpg
R=0.5
G=0.5
B=0.5
http://photos1.blogger.com/x/blogger2/2907/859880749988154/1600/500246/05.jpg

R=0
G=1
B=1
http://photos1.blogger.com/x/blogger2/2907/859880749988154/1600/859043/bg.jpg

7. Counter image

http://photos1.blogger.com/x/blogger2/2907/859880749988154/1600/872523/count.jpg

8. Color ro B/W
http://photos1.blogger.com/x/blogger2/2907/859880749988154/1600/673361/gra.jpg

結語：
目前只是最基本的修改功能
未來計畫要加入的包括進接的修改功能及目前富士最新的臉部辨識功能

程式碼

file=input('select your photo (xxx.xxx)\n:','s');
w=imread(file);
k=menu(' Matphoto v.1 by Ching','View photo','RGB filter','Counter image','Color ro B/W');
switch k
case 1
imview(w);
case 2
r=input('Input your R value(0~1)\n:');
g=input('Input your G value(0~1)\n:');
b=input('Input your B value(0~1)\n:');
w2(:,:,1)=w(:,:,1)*r;
w2(:,:,2)=w(:,:,2)*g;
w2(:,:,3)=w(:,:,3)*b;
imshow(w2);
case 3
w2=256-w;
imshow(w2);
case 4
w2=(w(:,:,1)+w(:,:,2)+w(:,:,3))/3;
imshow(w2);
end

熱傳模組建立（續）

繼上次的穩態一般邊界
這次加入穩態所有邊界情形

11/28/2006

小程式月曆產生器


把程式加入介面


  GUIDE程式說明
    首先在command window下打guide就會出現下圖

  接下來選Blank GUI(Default)，就會開啟新的window，如下圖

接下來你會看到左邊有很多控制物件，
你只要以滑鼠點選其中一項
並把它拖到你想要的位置就可以.
    月曆這個程式我有用到[push button]二個(run 和run)、
九個[text](如萬年曆、月、年和藍底的框都是)，
而其它的空白框則都是[edit](可以key in數字的)。
如下圖
(小格子是中文字，因為我7版的不能用，所以我就用6.5版的
，在此要注意，應該會有些指令在6.5不能用)。


而要是想更改其顏色，字體等…則在其物件點擊二下則會出現下圖(這是以紅色的run為例子)，
或是mark其物件，選View，再選Property Inspector也可以
  
當你選擇好你要用的物件，即可存檔。此時你會存成二個檔
一個叫 xxx.fig的檔，另一個是xxx.m的檔。接下來則是程式編輯的部份了,
點選View，再選M-file Editor即可開始編輯程式。
當你一開啟M-file你就會看到密密麻麻的字，但在這個程式會用到的只有幾個而已。
      大部份我們都是在function 物件名稱_callback(hObject, eventdata, handles)
下編輯其程式，如下圖：
這個程式的[edit物件名稱]空白處的程式碼為
function  edit物件名稱_Callback(hObject, eventdata, handles)
% hObject    handle to edit8 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
newstrval=get(hObject,'string');  % 截取輸入的字串
newval=str2double(newstrval);  % 變成數字
handles.edit物件名稱=newval;  % 存儲起來
guidata(hObject,handles);
% Hints: get(hObject,'String') returns contents of edit8 as text
%        str2double(get(hObject,'String')) returns contents of edit8 as a double
以下為pushbutton的圖
  這個程式的[pushbutton物件名稱]的程式碼為  
Function  pushbutton物件名稱_Callback(hObject, eventdata, handles)  
% hObject    handle to pushbutton3 (see GCBO)  
% eventdata  reserved - to be defined in a future version of MATLAB  
% handles    structure with handles and user data (see GUIDATA)  
Y=handles.edit6; % 截取年的值
M=handles.edit7; % 截取月的值
K=handles.edit8; % 截取日的值  
% 以下為判斷你輸入的日子是星期幾  
if rem(Y,100)==0 & rem(Y,400)==0 | rem(Y,4)==0
    d1=[31 29 31 30 31 30 31 31 30 31 30 31];
else
    d1=[31 28 31 30 31 30 31 31 30 31 30 31];
end
s=Y-1;
c=d1(1:(M-1));
c1=sum(c)+K;
ss=s+fix(s/4)-fix(s/100)+fix(s/400)+c1;
dd=rem(ss,7)+1;
di=['日','一','二','三','四','五','六']; 
ggo=sprintf('你查詢的%d年%d月%d日是星期%s',Y,M,K,di(dd)); 
set(handles.text12,'string',ggo);  % 把它show出來在藍色框框中如下圖


另一個[pushbutton物件名稱]的程式碼為
% --- Executes on button press in pushbutton4.
function  pushbutton4_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
Y=handles.edit9; % 讀取年的值
M=handles.edit10; % 讀取用的值
ans=[ ];
if rem(Y,100)==0 & rem(Y,400)==0 | rem(Y,4)==0
    d1=[31 29 31 30 31 30 31 31 30 31 30 31];
else
    d1=[31 28 31 30 31 30 31 31 30 31 30 31];
end
s=Y-1;
c=d1(1:(M-1));
% 以下為建月曆的運算式
for k=1:d1(M)
    c1=sum(c)+k;
    ss=s+fix(s/4)-fix(s/100)+fix(s/400)+c1;
    ans(k)=rem(ss,7);        
end
ans(ans==0)=7;
pp=ans(1);
ans1=[ ];
k1=0;
for n=1:6;
    for m1=pp:7
        ans1(n,m1)=k1;
        k1=k1+1;    
    end
    pp=1;
end
ans1(ans1>d1(M))=0; % 結果
 % 以下就是把所出來的值給代入小格子
(這個東東，還可以修改)  
set(handles.e1,'string','日');
set(handles.e2,'string','一');
set(handles.e3,'string','二');
set(handles.e4,'string','三');
set(handles.e5,'string','四');
set(handles.e6,'string','五');
set(handles.e7,'string','六');  
set(handles.e8,'string',num2str(ans1(1,1)));
set(handles.e9,'string',num2str(ans1(1,2)));
set(handles.e10,'string',num2str(ans1(1,3)));
set(handles.e11,'string',num2str(ans1(1,4)));
set(handles.e12,'string',num2str(ans1(1,5)));
set(handles.e13,'string',num2str(ans1(1,6)));
set(handles.e14,'string',num2str(ans1(1,7)));  
set(handles.e15,'string',num2str(ans1(2,1)));
set(handles.e16,'string',num2str(ans1(2,2)));
set(handles.e17,'string',num2str(ans1(2,3)));
set(handles.e18,'string',num2str(ans1(2,4)));
set(handles.e19,'string',num2str(ans1(2,5)));
set(handles.e20,'string',num2str(ans1(2,6)));
set(handles.e21,'string',num2str(ans1(2,7)));  
set(handles.e22,'string',num2str(ans1(3,1)));
set(handles.e23,'string',num2str(ans1(3,2)));
set(handles.e24,'string',num2str(ans1(3,3)));
set(handles.e25,'string',num2str(ans1(3,4)));
set(handles.e26,'string',num2str(ans1(3,5)));
set(handles.e27,'string',num2str(ans1(3,6)));
set(handles.e28,'string',num2str(ans1(3,7))); 
set(handles.e29,'string',num2str(ans1(4,1)));
set(handles.e30,'string',num2str(ans1(4,2)));
set(handles.e31,'string',num2str(ans1(4,3)));
set(handles.e32,'string',num2str(ans1(4,4)));
set(handles.e33,'string',num2str(ans1(4,5)));
set(handles.e34,'string',num2str(ans1(4,6)));
set(handles.e35,'string',num2str(ans1(4,7)));  
set(handles.e36,'string',num2str(ans1(5,1)));
set(handles.e37,'string',num2str(ans1(5,2)));
set(handles.e38,'string',num2str(ans1(5,3)));
set(handles.e39,'string',num2str(ans1(5,4)));
set(handles.e40,'string',num2str(ans1(5,5)));
set(handles.e41,'string',num2str(ans1(5,6)));
set(handles.e42,'string',num2str(ans1(5,7)));  
set(handles.e43,'string',num2str(ans1(6,1)));
set(handles.e44,'string',num2str(ans1(6,2)));
set(handles.e45,'string',num2str(ans1(6,3)));
set(handles.e46,'string',num2str(ans1(6,4)));
set(handles.e47,'string',num2str(ans1(6,5)));
set(handles.e48,'string',num2str(ans1(6,6)));
set(handles.e49,'string',num2str(ans1(6,7)));    
  下面為顯示結果圖

等程式都key in完，就可以按run，即可執行程式
月曆產生器檔案下載

11/24/2006

熱傳模組建立

this a test

11/13/2006

簡單的濾波器製作

　　對於擷取肌電圖（EMG）時，擷取到的資料對於我們而言
並不是都是有意義的，有時候常常會收取到一些雜訊，例如直
流電所產生的60Hz的雜訊，或是其實我們收取EMG的訊號以取
樣頻率1000為例，大多介於20-500之間。
　　因此濾波對於我們肌電訊號的分析格外重要。這裡撰寫了
一個非常簡單的濾波功能。有興趣的話可以點我看看。或是按
標題也可以連結。

11/11/2006

請問您有關MATLAB的問題楊淑媛

老師您好

抱歉打擾您了

想請問您有關"比較資料作累加動作:的問題

data2是k*3的矩陣

現在的比較方式是

當第三行的其中一個數比另一個數小

其對應的第二行資料卻比另一個數對應的第二行資料大時

則p+1

相反的

當第三行的其中一個數比另一個數大

其對應的第二行資料也比另一個數對應的第二行資料大時

則m+1

附上的是我寫的程式

基本的debug是沒錯

但它跑很久

(雖然說k=1000)

因此我怕這是一個無窮迴圈

請老師幫我看一下

謝謝老師

抱歉打擾您

麻煩您了

學生淑媛

11/03/2006

黃世榮:經緯儀的定位

前言:早前出差去使用經緯儀，希望能用兩種角度去作定位
問題:如何在只有兩組經緯儀的角度下作空間點群的定位，在建立公式以及程式後，發覺程式當時剛學技巧有需要改進且似乎這樣的計算是有錯的，希望大家幫忙想
程式流程:
clear
format long
matrix=xlsread('diow.xls','sheet1','A2:L3'); %% read the data
[m,n]=size(matrix); %% get row and column
L=input('Input the L you get==> ')
%%request the distance between two measure point

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% to build the alpha's array
for i=1:m
A(i)=matrix(i,1)
%% to build the alpha's array
A2(i)=matrix(i,2)
%% to build the alpha's minute array
A3(i)=matrix(i,3)
%% to build the alpha's second array

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% to build the theta's array
S(i)=matrix(i,4)
%% to build the theta's array
S2(i)=matrix(i,5)
%% to build the theta's minute array
S3(i)=matrix(i,6)
%% to build the theta's second array

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% to build the Beta's array
B(i)=matrix(i,7)
%% to build the Beta's array
B2(i)=matrix(i,8)
%% to build the Beta's minute array
B3(i)=matrix(i,9)
%% to build the Beta's second array

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% to build the PHI's array

F(i)=matrix(i,10)
%% to build the PHI's array
F2(i)=matrix(i,11)
%% to build the PHI's minute array
F3(i)=matrix(i,12)
%% to build the PHI's second array
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% to make their unit coherence
for i=1:m
A(i)=A(i)+(A2(i)/60)+(A3(i)/3600)
S(i)=S(i)+(S2(i)/60)+(S3(i)/3600)
B(i)=B(i)+(B2(i)/60)+(B3(i)/3600)
F(i)=F(i)+(F2(i)/60)+(F3(i)/3600)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%find that limit it or not is no matter
for i=1:m
% if B(i)<90
C(i,1)=L*sin(B(i)/180*pi)/sin((B(i)-A(i))/180*pi)
% to build X's array (X,Y are not the real coordinate)
C(i,2)=L*sin(A(i)/180*pi)/sin((B(i)-A(i))/180*pi)
% to build Y's array
% end
% if B(i)>90
% C(i,1)=L*sin(B(i)/180*pi)/sin((B(i)+A(i))/180*pi)
% to build X's array
% C(i,2)=L*sin(A(i)/180*pi)/sin((B(i)+A(i))/180*pi)
% to build Y's array
% end
end

%% to build coordinates' matrix
for i=1:m
H(i,1)=C(i,1)*sin(A(i)/180*pi)
H(i,2)=C(i,1)*cos(A(i)/180*pi)
H(i,3)=C(i,1)*tan(S(i)/180*pi)
H(i,4)=C(i,2)*sin(B(i)/180*pi)
H(i,5)=C(i,2)*cos(B(i)/180*pi)+L
H(i,6)=C(i,2)*tan(F(i)/180*pi)
end

%output the data to excel
xlswrite('test',H,'sheet1','A2');

for i=1:m
plot3(H(i,1),H(i,2),H(i,3),'o');
plot3(H(i,4),H(i,5),H(i,6),'o');
hold on
end

吳子青: 試作搜尋引擎( I )

前言：
本來構想是希望作一個可以將老師部落格上面
關於Matlab程式碼斷行問題解決的程式

利用copy指令將部落格的程式碼複製下來
會發現變成一整行老師在上課的時候也不太方便
構想:
因為學到cell() 本來希望將整個文件按照字元存成
cell, 再去對矩陣作處理, 把斷行點找出來之後
重新輸出成正確的文件
但是...
找不到可以斷行的依據

所以轉換方向作一個網頁或是文件的搜尋器 (google..??)
目前進度緩慢，所以分階段進行

問題
如何將字元讀出並進行搜尋動作？

----------------------------------------------------------
答案:
題旨與分析
將文件讀入並將文件裡面的每個字元存入記憶體宣告成矩陣
(或者文件太大以動態搜尋方式不先將整個文件全部存入記憶體）
對矩陣作搜尋動作, 找出關鍵字對應的位置
將搜尋到的字元或位置進行後置處理
(列出一整行,或者標示不同顏色？）
作出使用者介面

程式流程
因為尚未完成所以先簡述如下

/tra.m/
[fid,message]=fopen('test2.txt','r');
mydata=fread(fid);
if fid==-1
disp(message);
end
fclose all;

%如此就將文件文字全部以二進位儲存方式
%儲存成矩陣mydata

/dis.m/
disp(char(mydata'))

/run.m/
[m,n]=find(mydata==keyword);
%找出矩陣中的keyword的二進位元,在mydata中作搜尋
找出對應的位置
A={mydata(m,n-2) mydata(m,n-1) mydata(m,n) mydata(m,n+1) mydata(m,n+2) };
disp(char(A));
%將關鍵字元的二進位碼轉換, 並顯示前兩字後兩字

執行結果
test2.txt 內容 ghuodjgpobinhindpohinjestpohdbvuisneptogindbuodb
關鍵字p
結果 :
jgpob , ndpoh , stpoh , nepto

討論
這目前只是初步的概念 ,只是針對字串中的字母搜尋關鍵字母並作處理
未來要克服的包括關鍵字串的設定, 文件檔案讀取更general, 及加入
if 等規則等等等

吳堉光：EMG之Onset時間點分析

物理治療研究中，肌電圖（electromyography）是非常重要的研究工具之一，以下是我利用所教的while loop找出肌電訊號開始的時間點
EMG時間點

11/02/2006

林書如：您的體重過重嗎

前言：
BMI指數可以反映出體重是否處於理想狀態。過去已有許多文獻指出BMI指數過高往往是造成糖尿病、高血壓，甚至是心血管疾病的危險因子。因此，適度地監控自己的BMI指數是否處於正常值內是很重要的。

問題：
試寫出一個函數，當輸入體重和身高的數值後，能計算出其身體質量指數BMI值。並且依據其BMI值判斷其是否有體重過重或過輕的問題。(18.5<=BMI<24時顯示"體重標準"，BMI>24時顯示"體重過重"，BMI<18.5時顯示"體重過輕")

指令：

一開始我寫了一個比較簡單敘述檔案(Script file)

function [A]=BMI(w,h)
% BMI(Body Mass Index) calculation
A=w/h^2; % A means BMI
if 18.5<=A & A<24
fprintf('BMI=%0.2f 體重標準\n',A)
elseif A<18.5
fprintf('BMI=%0.2f 體重過輕\n',A)
else
fprintf('BMI=%0.2f 體重過重\n',A)
end

模擬輸出結果：

>> BMI(60,1.5)
BMI=26.67 體重過重

ans =

26.6667

結果顯示，一位體重60公斤，身高150公分的人，其BMI指數為26.67，有體重過重的傾向

==============================================

但是，為了能夠加上交談的功能，讓使用者更能清楚地了解應該輸入的內容
例如，輸入之體重應以公斤為單位，輸入之身高應以公尺為單位
因此我又改寫成另一個函數檔(function file)，如下

% BMI2.m
% BMI(Body Mass Index) calculation
a=input('Please keyin your body weight? (kg)');
if isnumeric(a)
b=input('Please keyin your body height? (m)');
B=a/b^2; % B means BMI
if 18.5<=B & B<24
fprintf('BMI=%0.2f 體重標準\n',B)
elseif B<18.5
fprintf('BMI=%0.2f 體重過輕\n',B)
else
fprintf('BMI=%0.2f 體重過重\n',B)
end
end

模擬輸出結果：

>> BMI2
Please keyin your body weight? (kg)60
Please keyin your body height? (m)1.5
BMI=26.67 體重過重

結果同上，一位體重60公斤，身高150公分的人，其BMI指數為26.67，有體重過重的傾向

討論：

我是以簡單的if-elseif指令來撰寫這兩個小程式
程式中如果還有可以補強的地方(例如發現程式的bias，或是其實還有更簡單的寫法)，還麻煩大家提供意見 :p

10/25/2006

吳子青: 民意調查真的公正嗎

前言：
最近又到了選舉期了在電視媒體或報章雜誌常常可以看到
“ 某某參選人支持率多少”…
民意支持率各家都宣稱客觀公正，但其實是可以操作的
在這裡舉出一個以前碰過的問題讓大家思考一下　

問題：

若台北市某候選人的真實 “民意支持率”為15%，找出使調查樣本的
民意支持率超過0.18 的出現機率不大於 10% 的最小樣本數 N。

----------------------------------------------------------------------
答案：

題旨與分析
這問題的思考流程是這樣的
將全部台北市民當成一個族群，而調查的對象代表抽樣的樣本
“某候選人的真實民意支持率為15%” 這句話的意思就是代表族群的支持率
而調查樣本的民意支持率就是樣本的支持率。

將這問題視為一個二項式分佈的問題，也就是支持和不支持兩種情況
支持的機率 P=0.15 不支持的機率 Q=0.85
舉例說明
如果某一次的抽樣從台北市民隨機抽10個人 (Sample Size m=10)
10個人中支持候選的人有2個人，不支持的有8個人
這個事件出現的機率為

 http://photos1.blogger.com/blogger2/2907/859880749988154/1600/formula1.jpg

所以如果問題是 ”找出使調查樣本的民意支持率超過0.18”
如果抽樣人數維持10人，支持人數大於1.8的都符合條件
也就是支持人數 2 到10人分別算出的機率加總。

題目所求是 “出現機率不大於 10% 的最小樣本數 N。”所以從n=1
開始算到符合出現機率不大於10％為止。

程式流程
Step 1 : 宣告function p1(m,n) 為抽樣m個人，有n個人支持的機率

Step 2 : 計算A為抽樣m個人，支持率超過0.18的機率總和

Step 3 : 從m=1 開始增加用For 跑迴圈，決定停止條件為A < >程式內容


/ p1.m / (* 這行代表m-files p1.n 開始 *)

function output=p1(m,n)

if n>0.5*m
a=m-n;
output=prod((m-a+1):m)/prod(1:a)*(0.15)^n*(0.85)^(m-n);

else
output=prod((m-n+1):m)/prod(1:n)*(0.15)^n*(0.85)^(m-n);
end

/ pp1.m / (* 這行代表m-files pp1.n 開始 *)
clear;
max=400;
ppp1=zeros(max,2);

for m=1:max
k=ceil(m*0.18);
A=0;

if A<=0.1 fprintf('最小樣本數%d\t出現機率%d',m,A);
break
end
fprintf('最小樣本數%d\t出現機率%d\n',m,A);
clc;
end
plot (ppp1(1:max,2), 'DisplayName', 'ppp1(1:max,2)', 'YDataSource', 'ppp1(1:max,2)');
figure(gcf)
xlabel('樣本數目');
ylabel('出現機率');
/ Command windows / (* 這行代表command windows輸入 *)
pp1

執行結果
最小樣本數217 出現機率9.559103e-002
這是最後跑出來的結果
N=217為民意支持率超過0.18 的出現機率不大於 10% 的最小樣本數 N

這是程式結果

http://photos1.blogger.com/blogger2/2907/859880749988154/1600/hy2.jpg

這是圖

http://photos1.blogger.com/blogger2/2907/859880749988154/1600/hydro.jpg

討論
(1)因為老師上週的上課內容剛交到For 和If的用法，所以在這裡練習一下

(2) 程式中碰到最大的問題在於 prod () 連乘函數括號內的矩陣不能太大
從1＊2＊3……..＊171 ＝prod(1:171) 這樣結果會爆掉，可是在算二項式分佈
的機率時就一定會用到。

所以在這裡用了一個小技巧，在p1裡當 n > 0.5 m 時令n=m-n
因為在排列組合裡 C10取2=C10取8，10!/8!2! =10!/2!8!
在這裡可以避免連乘太多項，可以有效避免算到M=171時，
結果一定會出現NaN。但是當m 算到300多時
仍然會出現同樣的問題。或許會有更有效率的寫法。

訂閱：文章 (Atom)