睡覺也是有方法的

不知道的人也許會認為，睡覺就是睡覺，一次睡8個小時跟分段睡8個小時並沒什麼不同。
但是事實上，一次睡足8小時跟分段睡滿8小時是絕對不同的。

正常人的睡眠有所謂的睡眠週期，每一個週期內又有不同的睡眠階段，
睡眠中的每一個階段，各有不同的用處。
分段睡眠會打亂睡眠週期，睡眠的品質可是會大打折扣的。

剛剛查資料時，找到一篇與睡眠品質有關的好文，並且還滿適合一般人閱讀的，
這篇文章的標題是「優質睡眠」，作者從睡眠的相關科學解釋，來教大家睡覺應該怎麼睡。
因為睡眠不只是單純的休息而已，還包括生理機能的修復、免疫系統的運作、甚至於 REM 期可以幫助記憶。睡不好不只是會勞累而已，對身體的長期影響也是很大，現代人應該要多了解喔。

優質睡眠︰http://mypaper.pchome.com.tw/news/kiske/3/1295571123/20070924140355/#centerFlag

The R Project for Statistical Computing

R 是著重於統計計算及繪圖的程式語言及環境，相當類似於Matlab。R 本身是一個 GNU project ，透過 自由軟體基金會 的 GNU General Public License 所發佈，並提供 Linux, MacOS X 及 Windows (95 and later) 三個平台的版本。因為 R 是以開放原始碼的方式發行，有許多社群或個人紛紛提供他們寫作的程式套件(package)。

R 可以描繪各種不同的統計圖，包括描繪3D圖形，官網中有放一些 截圖 可以參考。然而 R 本身並不提供 GUI(Graphical User Interface) 的環境，因此所有的工作都必須透過指令來完成。對非有 GUI 不可的使用者，恐怕會失望了 :D。

官網的文件雖然多，但是文件太多，反而讓人有點摸不清方向，只能以霧煞煞來形容(這是我的感覺啦，因為我也是新手)。A Quick, Painless Tutorial and Reference on the R Statistical Package 算是一個入門用的簡易手冊，不過提及的東西還是滿多的，新手學習應可從此文件著手。

Visual C++ 2005

手賤安裝了 Windows Vista ，結果發現 Visual C++ 6.0 不相容，目前是嘗試看看 Visual C++ 2005 能不能取代，不然要灌回 Windows XP 了。

以下列出目前遇到一些問題或麻煩事項及解決方式︰

1. 首先 Visual Studio 2005 必須先升級到 SP1 ，再另外安裝一個 Visual Studio 2005 SP1 Update for Vista

2. 開啟 Visual C++ 時，會出現視窗建議你以系統管理員身份啟動它。方法是在 Visual C++ 2005 的捷徑上按右鍵，再選 "以系統管理員身份執行" 。

3. 許多老舊函式由於有安全性問題，目前有許多替代的函式可使用，不過暫時實在不想去動它，所以在專案中加上一個 macro : _CRT_SECURE_NO_DEPRECATE ，以忽略這些煩人的警告(warnning)。

4. fstream 在讀、寫檔案時，如果遇到檔名內有中文字，在預設情況下會開啟失敗。這跟 locale 有關係，必須使用 setlocale() 來設定目前程式使用的語言。例如我的設定是︰

setlocale(LC_ALL, "Chinese");

我是把它放在 CWinApp::InitInstance() 函式內的第一行，目前讀、寫檔案看起來是都正常。

5.VC++ 6.0中好用的Resource Editor，我在 VC++ 2005還找不到怎麼叫出來，是該去買書回來看看了..

C++的 #pragma pack(n) - 位元組對齊 or 字節對齊

這東西的中文名稱似乎沒有統一，我也是第一次看到。

今天學弟在 C++Builder 下嘗試要把已有的 BITMAPINFO, 以及 Bitmap 的影像陣列儲存到 bmp 格式的檔案中，我看了看覺得，既然已經有 BITMAPINOF 跟 DIB 格式的影像陣列，那不如直接用 ofstream::write 寫進檔案裡。但是試了好久，卻一直無法輸出正確的檔案。

問了一下 Google 大神，馬上找到一篇

http://topic.csdn.net/t/20020402/07/615811.html

不過文中所說的方法必須要修改 C++Builder 內建的檔案，因為擔心會有後遺症，所以我決定自己寫 BITMAPINFO, BITMAPFILEHEADER, BITMAPINFOHEADER 來處理，但是對齊的問題我也會遇到。對齊的概念及會產生的問題請參考 字節對齊 以及 請問C 裡面的 #pragma 這兩篇文章。

32Bits 系統預設是 #pragma pack(4) ，而 BITMAPFILEHEADER 這個結構的長度為 14 Bytes，為了達到預設對齊的結果，編譯器會將 BITMAPFILEHEADER 的長度拉長成 16 Bytes (4 的倍數)，造成輸出的檔案的標頭檔會有問題，無法被讀圖軟體解析。為了避免這個問題，必須在宣告結構之前，將位元組對齊重設為 2 ︰

#pragma pack(2)

而在宣告結束後，重設為預設值︰

#pragma pack()

=> 在 pack() 的括號中不給值，即會重設回預設值。

如此一來，由於 14 為 2 的倍數，因此 BITMAPFILEHEADER 的結構就不會再受到影響。

Inverse of a matrix - 使用 LU 演算法

在 GSL 中已經實作了矩陣的反轉換函式

使用 GSL 定義的矩陣
首先應該知道怎麼在 GSL 中使用其定義的矩陣，這部份應參考 GSL Reference Manual 的 Vectors and Matrices 一節的 Matrices 小節。必須注意的是，GSL採用的是 C 語言的架構，即使用 struct 來實作矩陣 gsl_matrix ，再透過許多函式來操作之。

1. 宣告及配置記憶體給矩陣

/* 宣告及配置記憶體 */
gsl_matrix * m1 = gsl_matrix_alloc(m, n);

/* 宣告及配置記憶體，並將所有元素初始化為0 */
gsl_matrix * m2 = gsl_matrix_calloc(m, n);

　ps. gsl_matrix 宣告的矩陣，其元素的型態為 double ，若要使用其他型態，
　　 可使用 gsl_matrix_float, gsl_matrix_int, .. 等。

2.初始化矩陣

/* 設定矩陣 m1 內所有元素的值為 1.0 */
gsl_matrix_set_all (m1, 1.0);

/* 設定矩陣 m1 內所有元素的值為 0 */
gsl_matrix_set_zero (m1);

/* 將矩陣 m1 內對角線的元素值均設為 1 ，其餘為 0，即為單位矩陣 */
gsl_matrix_set_identity (gsl_matrix * m);

3.讀取及寫入矩陣元素

/* 讀出矩陣m1的第i列(row)第j行(column)的值 */
double v1 = gsl_matrix_get(m1, i, j);

/* 將矩陣m1的第i列(row)第j行(column)的值設為 5 */
gsl_matrix_set(m1, i, j, 5);

/* 直接傳回矩陣m1的第i列(row)第j行(column)元素的指標 */
double * ptr = gsl_matrix_ptr(m1, i, j);

　ps. 預設的情況下，這3個函式會檢查傳入的 i, j 值範圍，在大量運算的時候會拖慢系統執行速度，此時可以加上 preprocessor definition GSL_RANGE_CHECK_OFF來關掉此功能。

4.矩陣的操作

　GSL 在矩陣操作上提供了矩陣相加(gsl_matrix_add)、相減(gsl_matrix_sub)、乘上一個常數(gsl_matrix_scale)、加上一個常數(gsl_matrix_add_constant)等函式，比較令人不解的是，此處提供的 gsl_matrix_mul_elements 函式，並非正常的矩陣相乘運算，因此矩陣相乘的工作必須自行另外撰寫程式碼運作。這個程式碼很簡單，以下是我實作的矩陣相乘︰

#include <gsl/gsl_matrix.h>
#include <assert.h>

void gsl_matrix_mul_matrix(const gsl_matrix * a,
                      const gsl_matrix * b,
                      gsl_matrix * c)
{
assert(a->size2 == b->size1
   && a->size1 == c->size1
   && b->size2 == c->size2);

/* set all elements in c to zero */
gsl_matrix_set_zero(c);

/* calculate mul */
int i, j, k;
for(i = 0; i < a->size1; i++)
{
   for(j = 0; j < b->size2; j++)
   {
       for(k = 0; k < a->size2; k++)
       {
           gsl_matrix_set(c, i, j, gsl_matrix_get(c, i, j)
               + gsl_matrix_get(a, i, k)
               * gsl_matrix_get(b, k, j));
       }
   }
}
}

5.使用完陣列記得釋放記憶體

/* 釋放掉為矩陣 m1 配置的記憶體 */
gsl_matrix_free(m1);

計算矩陣的反轉換

矩陣的反轉換函式定義在 Linear Algebra 一節的 LU Decomposition 小節中。假設要計算反矩陣的是 n x n 的矩陣 A，則︰

1.計算 LU Decomposition

gsl_permutation * p = gsl_permutation_alloc(n);
int sign = 1;
gsl_linalg_LU_decomp(A, p, &sign);

　此時 LU Decomposition 的結果儲存在 A 及 p 中。

2.計算原矩陣 A 的反矩陣

gsl_matrix * inverse = gsl_matrix_alloc(n, n);
gsl_linalg_LU_invert(A, p, inverse);

　此時矩陣 inverse 即為原矩陣 A 的反矩陣。然而必須注意的是，經過 gsl_linalg_LU_decomp的計算，A中儲存的已經不再是原矩陣，故使用 gsl_linalg_LU_decomp 函式前應先使用gsl_matrix_memcpy 函式備份原矩陣。

gsl_matrix * A_c = gsl_matrix_alloc(n, n);
gsl_matrix_memcpy(A_c, A);

Reference
GSL Reference Manual - HTML
Matrix Inverse
Identity Matrix

GSL - Random Number Generation

Random Number Generation

剛剛學弟問到，所以看了一下怎麼使用 GSL 的 Random Number Generation 函式。

GSL 內建多個產生隨機變數的演算法，詳情請見 Random number generator algorithms，預設情況下，使用的 Generator 是 mt19937。

Random Number Generation 使用分成三個步驟︰

1. 初始化隨機變數產生器
　(1) 指定隨機變數產生器(Generator)，下面範例指定使用 taus 演算法

gsl_rng * r = gsl_rng_alloc (gsl_rng_taus);

　(2) 設定種子(Seed)，下面範例指定種子為 123

gsl_rng_set(r, 123);

2. 取得隨機數 - 重複執行以下函式，其傳回值即為隨機數
　(1) gsl_rng_get (const gsl_rng * r) : 取得隨機數，其隨機值範圍因不同的 Generator 而不同。
　(2) gsl_rng_uniform (const gsl_rng * r) : 取得 [0,1) 範圍內的浮點數，不包含 1。
　(3) gsl_rng_uniform_pos (const gsl_rng * r) : 取得 (0,1) 範圍內的浮點數，不包含 0 及 1。
　(4) gsl_rng_uniform_int (const gsl_rng * r, unsigned long int n) : 取得 [0, n-1] 之間的整數值。

範例︰
for (i = 0; i < n; i++)
{
　double u = gsl_rng_uniform (r);
　printf ("%.5f\n", u);
}

3. 釋放記憶體 - gsl_rng_alloc 會配置記憶體給 gsl_rng * r ，使用完必須釋放記憶體。

gsl_rng_free (gsl_rng * r)

C/C++ GNU Scientific Library(GSL) for Windows

The GNU Scientific Library (GSL) is a numerical library for C and C++ programmers. It is free software under the GNU General Public License.

The library provides a wide range of mathematical routines such as random number generators, special functions and least-squares fitting. There are over 1000 functions in total with an extensive test suite.

GNU Scientific Library (GSL) 是一個內含許多數值及科學運算函式的 C/C++ 函式庫，內含超過 1000 個以上的函式(例如︰數值微分、小波轉換、排序..等等) ，並且可以使用在各種作業系統平台。
在 GSL 在 Windows開發環境下的使用 中有簡單說明如何在 Windows 環境使用 GSL ，以下則是再詳細一點介紹如何在 Visual C++ 6.0 的 IDE (整合型視窗介面) 下使用︰

1. 首先到 http://gnuwin32.sourceforge.net/packages/gsl.htm 下載 Binaries、Developer files 。
2. 將 Binaries package 內 bin 子目錄下的 libgsl.dll, libgslcblas.dll 兩個檔案複製到 C:\Windows\System32 中。
3. 將 Developer files 解壓縮後，在 VC 的 IDE 中設定 include, lib 子目錄的路徑。
4. 由於此版本的 lib 子目錄中不含 .lib 檔，必須用以下指令產生。請先開啟 DOS 模式視窗(附屬應用程式/命令提式字元)，切換目錄到 lib 目錄下，並下達︰

lib /machine:i386 /def:libgsl.def
lib /machine:i386 /def:libgslcblas.def

即可產生 libgsl.lib, libgslcblas.lib 兩個檔案。
5. 在 VC 專案中加入 libgsl.lib, libgslcblas.lib 這 2 個 lib 到 link 參數中。
6. 在 VC 專案的 Preprocessor definitions 中加入 GSL_DLL。
6. 可用函式及說明請見 Reference Manual 。

ps. GSL是使用C語法寫成，故並沒有使用類別(class)，而是以結構(struct)及函式(function)組成，不過 GSL 仍可以用在 C++ 編繹器，並與 C++ 程式相容。

訊號內插 - Linear Interpolation 與 Sinc Interpolation

一串取樣間隔為 0.1 秒的時間訊號為例，如果想以訊號處理方法將訊號間隔增加到 0.025 秒，自然必須以內插加入新資料點。

Linear Interpolation

增加資料點最簡單的方式是使用 Linear Interpolation ，其方法是將新取樣點前後想鄰兩原取樣點連線，連線通過新取樣點時間位置時的值即設為新取樣點的值。如下圖所示，A、B為原取樣位置，SA及SB為原取樣位置的取樣值，N為新取樣位置，SN為新取樣位置的取樣值，則SN的計算公式如下︰

Sinc Interpolation

假設一訊號 f(t) 及其頻譜 F(ω) 如下圖︰

有一個 impulse train δT(t) ，其 pulse 之間的間隔為 T ，其訊號圖如下︰

若訊號 f(t) 以 T 的時間間隔做取樣，則取得取樣後的訊號 f(t) = f(t)．δT(t) ，其訊號及頻譜如下︰

從頻譜的角度來看，訊號的取樣過程會在頻譜上產生週期化的現象。因此若如上圖頻譜圖所示，在取樣後的頻譜上，使用一個頻寬與原訊號頻寬相同的低通濾波器(Lowpass filter)濾波，則又可將取樣後的訊號還原到原訊號。下圖為一理想低通濾波器的頻譜H(ω)及時域圖h(t)︰

則還原後訊號 fR(t) = f(t) * h(t) ，重建過程及重建後訊號如下圖︰

理想的低通濾波器在時域的型式為一 sinc 函數，然而 sinc 是一無限的函數(unbounded function)，只能以近似的方式模擬，因此實作時所還原得到的訊號與原訊號仍有差異。

若一取樣間隔為 T1 的離散訊號，欲以 T2 的取樣間隔重新取樣，可假想為將取樣間隔為 T1 的離散訊號以理想低通濾波器還原為原來的連續訊號，再使用 T2 取樣間隔重新取樣，又由於理想的低通濾波器在時域的型式為一 sinc 函數，因此一方法稱為 Sinc Interpolation。

若
f(k1T1)︰原取樣間隔為 T1 的離散訊號
f(k2T2)︰新取樣間隔為 T2 的離散訊號
T1︰原訊號取樣間隔
T2︰新取樣間隔
B︰sinc 函數(or 低通濾波器)的頻寬，T1 = 1/(2B)

則 Sinc Interpolation 的計算公式為︰

Sinc Interpolation 的 C 程式碼請參考這裡。

參考資料︰
B. P. Lathi, Signal Processing & Linear Systems, Berkeley Cambridge Press, 1998.

FFTW - 計算 FFT 的函式庫 (C 語言)

FFTW 是用 C 語言寫成的函式庫，用以計算 DFT 轉換。輸入的訊號可以是任意長度，可以是實數型式或是複數型式的訊號。當然，從它的名稱知道它是以 FFT 來實作。本來我已經有修改自 NUMERICAL RECIPES 的 fft 轉換程式，不過 NUMERICAL RECIPES 的授權限制滿大的，相較之下 FFTW 則是採用 GNU General Public License 授權，比較可以放心使用。

要在 VC 下使用 FFTW ，可以從這裡下載 zip 檔。

安裝︰

解壓縮之後，開啟 命令提示字元 ，轉換目錄到解壓縮的目錄，再執行以下三個指令︰

lib /machine:i386 /def:libfftw3-3.def
lib /machine:i386 /def:libfftw3f-3.def
lib /machine:i386 /def:libfftw3l-3.def

這會在目錄下建立三個 lib 檔。

將 libfftw3l-3.dll, libfftw3f-3.dll, libfftw3-3.dll 這3個 dll 檔複製到 system32 目錄。
在 VC 專案中指定 libfftw3l-3.lib, libfftw3f-3.lib, libfftw3-3.lib 這3個lib檔及 fftw3.h 檔所在的目錄。

在 VC 專案中加入 libfftw3l-3.lib, libfftw3f-3.lib, libfftw3-3.lib 這3個 lib。

使用︰
首先要引入 fftw3.h 這個標頭檔，如果要做的是一維 DFT 轉換，可參考以下的程式碼(取自fftw的線上說明文件)
#include <fftw3.h>
...
{
fftw_complex *in, *out;
fftw_plan p;
...
in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
...
fftw_execute(p); /* repeat as needed */
...
fftw_destroy_plan(p);
fftw_free(in); fftw_free(out);
}

fftw_complex 是 FFTW 自訂的複數型態，是一個2個元素的浮點陣列 double[2]，若 in 為 fftw_complex 所宣告的一個陣列

fftw_complex * in = new fftw_complex[10];

，則第一個元素的實部為 in[0][0], 虛部為 in[0][1]。
fftw_malloc 是類似於 malloc 的函數，之間的差別我還不太清楚，不過文件上有提到，使用其他開記憶體的方式也是可以的(例如: new)，使用 fftw_malloc 的好處是，開啟的陣列元素(包括 fftw_complex 的實虛部)會初始化為0(並非如此)。

fftw_plan_dft_1d - 從函式名稱大概可以知道，這個函式是用來產生一個準備的工作，它並不會真的執行 DFT 轉換，真正的 DFT 轉換要再用 fftw_exectue 來執行。

多維度號的轉換只需要將 fftw_plan_dft_1d 轉為以下任一個函式。

fftw_plan fftw_plan_dft_2d(int nx, int ny,

fftw_complex *in, fftw_complex *out,

int sign, unsigned flags);

fftw_plan fftw_plan_dft_3d(int nx, int ny, int nz,

fftw_complex *in, fftw_complex *out,

int sign, unsigned flags);

fftw_plan fftw_plan_dft(int rank, const int *n,

fftw_complex *in, fftw_complex *out,

int sign, unsigned flags);
輸出訊號格式︰
fftw_plan_dft_1d :
若輸入訊號的長度為 32 ，則輸出的FFT訊號從第17個元素(從0起算)開始重複。
out[17][0] == out[15][0], out[17][1] == out[15][1]
out[18][0] == out[14][0], out[18][1] == out[14][1]
                       .
                       .
out[31][0] == out[2][0], out[31][1] == out[2][1]


fftw_plan_dft_r2c_1d :
這個函式的傳入陣列為實數陣列，若輸入訊號長度為32，則輸出的FFT訊號會在
第17個元素之後均為0

Critticall home page

Critticall home page

A programming tool
根據說明，這似乎是一個可以幫你把演算法程式最佳化的工具，
例如一個 Bubble Sort 的演算法，Critticall 可以把他最佳化，而且結果比常用的 Quick Sort 還要快。

滿神奇的東西..