簡單筆記一下今天工作的成果：將 Matlab 的矩陣噴出成 Amira 的 .am 檔案，也就是 Amira 3D mesh。
這個成果是根據我手上有的 .am 檔案還有 Matlab Central 上的 LoadData_Amira.m 逆推的，特此致謝。
雖然從 LoadData_Amira.m 看起來 AM 檔案有支援 uint8 之外的其他數值格式，但是我剛才弄半天弄不太出來，所以為了簡便，就先強制輸出為 uint8 吧，科科！
剛才突然想到，之前一直有個很差的習慣，就是遇到2D的矩陣，就用很彆扭的方式去解決，例如要解2D矩陣的 max ，就寫
A = rand(6, 3);
max_num = max(max(A));
照理來說min和mean這類 會從 column 開始做起然後才做更高維度 的函數，我也都是用一樣的方法來解決。
可是俗話說的好， 夜路走多了會遇到鬼，如果遇到 std 我就沒轍了，如果你這樣寫你就準備出糗了
有在用 Matlab 的人可能有用過 polar這個畫極座標的繪圖函數，但是很令人崩擴的是 polar 這個函數本身功能非常有限，如果要多別的事情，都要靠 handler 還有 set搭配來完成，不像 plot早就已經做到萬能的地步，幾乎想要做什麼，都可以直接在 argument 裡面完成。
現在的問題是，假設我想要把兩個函數都畫在同一張 polar 上面，聰明的各位可能馬上就會想到可以透過 hold on,來實作，這是沒錯的，但是如果實際去做，會發現一個很糟糕的事情：polar的limit會被第一個畫的函數所限制住，例如假設我想要畫以下的函數：
最簡單的寫法大概是這樣：
close all;
N = 101;
theta = linspace(0, 2*pi, N);
y1 = abs(cos(theta));
y2 = 1.5 * abs(sin(theta));
figure(1),
polar(theta, y1, 'k');
hold on,
polar(theta, y2, 'r');
hold off;
但是畫出來卻會變這樣：
囧！因為 limit 已經被第一個polar決定了，所以如果後面再hold on上去的圖案更大，就只能悲劇。
Amira 是一個許多做顯微術或者生物影像的人都會用到的軟體，可以用它來具象化許多3D的data。
除了原生的AM檔案之外，有時候我們也會希望藉由 Amira來顯現來自 Matlab我們加工後或者計算後的檔案，因此就會需要由 Amira來打開 .mat檔案；照理來說 Amira 是支援的，但是我自己實際測試，卻會噴出這樣的錯誤訊息：
Reading minivirus_A.mat ...
opening connection to Matlab engine...
Couldn't open libeng.dll:
File does not exist.
Couldn't open libmat.dll:
File does not exist.
Couldn't open libmx.dll:
File does not exist.
Failed to connect to Matlab
Be sure Matlab is correctly installed, and the path is valid
Error opening file G:/minivirus_A.mat
網路上查半天，都找不到問題所在或者解法，但是自己嘗試去解決，卻發現解法超乎想像地簡單。
首先先試著去找 libeng.dll 、libmat.dll或libmx.dll這三個鬼東西在哪。答案很簡單，就在
C:\Program Files\MATLAB\R2012b\bin\win64
的下面。我想 Amira 抓不到的原因非常簡單，是因為 Amira 不知道 64-bit 的 Matlab 已經改變存放這些lib的位置了。
解決方法
解法非常簡單，就是改系統變數中的 PATH值。
首先對著 我的電腦 按右鍵，選 內容 ，接著選 進階系統設定 ，然後按 進階 這個頁籤，最後按最下面的按鈕 環境變數 。接著在下面的小框框中找到PATH，點選 編輯，然後在最後面加上
;C:\Program Files\MATLAB\R2012b\bin\win64
就大功告成，此時重開 Amira 後就可以開啟 .mat 檔案了。
之前就有想過要用解析解求一個均勻的圓球的質量投影，但是推半天沒推出來，只好先用MATLAB做出一個3D的volume，再用 sum(V,3)來求投影質量。我知道這樣超白癡的，而且要是球稍微大一點，電腦就記憶體不足當掉給你看了。
剛才突然... 就夢到正確的解析解做法了......
原理
已知圓球球面的公式是\
其中 為一常數，是球面的半徑。
現在我們將其推廣，把球面推廣成球體，則參數 要滿足
其中
所以要快速求出一個球體的質量投影，只要移項還有乘以2就好了，得到：
不過記得要把超出 support 的地方手動設為零，不然可是會出現 complex 的。請直接參考下面的實作程式碼。
MATLAB 實作
function Z = sphere_proj(DIA)
RAD = DIA/2; % radius
x = linspace(-RAD, RAD, DIA);
[X, Y] = meshgrid(x, x);
R = sqrt( X.^2 + Y.^2);
Z = 2 * sqrt( RAD^2 - R.^2 );
Z( R > RAD) = 0;
簡介
訓練目標
以1D為切入點，了解 EEMD 之基本原理，並以 Matlab 語言實作之。
Matlab
完成這份訓練講義之後，預期學習者將會應用到或者學到如何
- 創造並使用
function - 找出資料向量中的極值
- 使用
interp1內插出平滑的資料曲線 - 打造 EMD 引擎
- 將 white noise 加入資料之中
- 完成 EEMD 函數
習題
尋找極值函數 findExtrema
在這題當中，我們要學會寫出一個傳入資料向量 data，可以回傳資料中極值的 值 以及 位置 的函數。以下將分為四個小題進行。
學會如何自定函數
試找尋資料，學習如何創造一個自定函數，或查詢 Matlab 內建教學
doc function
或參考以下範例：
function y = numAdd(a, b)
y = a + b;
尋找最大值
試寫一個函數 findMax.m，輸入任意長度向量 data 可傳回兩組輸出：
- 最大值所在的位置
max_ind - 最大值的值
max_val
提示： 函數的第一行可能長成這個樣子：
function [max_ind, max_val] = findMax(data)
% coding here
請注意，一串資料可以有多個 local maxima。此題要求的是找出所有的 local maxima，而非僅是資料中的最大值，否則用內建函數 max() 即可達成。
尋找最小值
承上題，找出「所有」最小值 min_val 及其分別之所在位置 min_ind。
提示： 試想，是不是將上一題的傳入資料由 data 改為 -data 就可以達成目標了？
尋找所有極值
綜合前面三題所學，請完成一個可以找出所有極值，並分別回傳最大最小值的值與位置的函數 findExtrema.m
function [max_ind, max_val, min_ind, min_val] = findExtrema(data)
% coding here
注意： data(1) 與 data(end) 是否也須考慮為「極值」？如果是，其原因為何？
內插曲線
下面這兩題，我們將會利用前一大題所得到的資料，以及 Matlab 內建的 interp1函數，得到 EMD 所需要的「趨勢」曲線。
使用離散資料內插曲線
利用 interp1，將 max_ind 和 max_val內插成一條長度和取樣率皆與 data 相同之「連續」曲線 max_trend。
利用 interp1，將 min_ind 和 min_val內插成一條長度和取樣率皆與 data 相同之「連續」曲線 min_trend。
注意： max_trend或min_trend的頭尾，都要和 end 一樣（可以想成是 boundary conditions），否則會有些區域沒有定義。
求得「平均」曲線
將上小題的兩條曲線做點對點平均，得到 avg_trend。
打造 EMD 引擎
參考 paper 上的演算法，以及上面所練習的技巧與成品，應該已經可以完成 EMD 的引擎。其界面如下：
function modes = emd(data, nmode)
% coding here
其中 nmode 是所要求的模態數量。請注意，如果原本的 data 長度是 的話，則得到的輸出資料 modes的大小會是 。
打造 EEMD 引擎
如果前面的 emd.m 已經驗證完成，那麼接下來就是讓它進化成 EEMD 的時候了。比起 EMD ， EEMD 多了下面兩項特色
- 外加了額外的 white noise 進去
- 先跑出 N 組各別使用不同的 white noise 所解出來的結果，再一起取平均。
下面我們只拆成兩小題練習，不足之處請直接參閱 paper。
加入 white noise
假定使用者指定一個特定的 white noise 強度 wn_std，我們有兩種方式可以將等效強度的 white noise Wn 加到 data 裡面去，第一種是：
其中 指 data， 是求標準差的函數；或者第二種方法
試問：兩種方法何種較優？較劣者，有什麼致命的缺點？
請利用 Matlab 內建的 randn() 完成此題。
完成 EEMD 架構
參考 paper 上的演算法，以及上面所練習的技巧與成品，終於我們可以完成 EEMD 的全面架構。完成品的輸入資料除了 data以及nmode之外，還多了指定 ensemble 數量的 nens 與 white noise 的標準化強度 wn_std。其界面如下：
function modes = eemd(data, nmode, nens, wn_std)
% coding here
