... ... ... ... ... ...

三、一維 FFT 演算法

1. 將一維離散傅立葉轉換改寫如下：

where

2. 假設 ， n 為正整數，則

where

and

for u = 0,1,2,…,M-1 ,

如此，原本要計算 f(x) 的傅立葉轉換 F(u)，已經分解為計算偶數編號資料 f(2x) 的傅立葉轉換 Feven(u) 及奇數編號資料 f(2x+1) 的傅立葉轉換 Fodd(u) 的合成。

3. 利用 且 可推得

如此， F(u) 與 F(u+M) 在運算時只相差一個減號，也就是在計算 F(u) 時，所求出的 Feven(u) 及 Fodd(u) 與的乘積，只要用加法運算，即可得到F(u)，以減法運算，即可得到 F(u+M) ，這樣在算 F(u) u=0,1,2,...,M-1時，可順便算出 F(u) u=M,M+1,...,2M-1 這樣運算量就大幅減少了。

4. 同樣地，在計算 Feven(u) 及 Fodd(u) 時，也是利用相同的方法。因為總資料數 N ＝ 2^n ，在分割 n 次後，即可由原資料中某兩筆資料的加減運算分別算出第一層的兩筆的 Feven 及 Fodd ， 這樣就可以少掉Σ的運算，使運算量再減少。

5. 範例：假設我們要由 f(0),f(1),f(2),f(3),f(4),f(5),f(6),f(7)

算出其傅立葉轉換 F(0),F(1),F(2),F(3),F(4),F(5),F(6),F(7)

由演算法得知當 N = 8 時,

u = 0,1,2,3

又 Feven(u) 是 f( ) 偶數編號 f(0),f(2),f(4),f(6) 的傅立葉轉換。

為使符號表示方便，令 G(u) ＝ Feven(u) for u = 0,1,2,3

所以，可將問題視為由資料 f(0),f(2),f(4),f(6)

算出其傅立葉轉換 G(0),G(1),G(2),G(3)

由演算法得知當 N = 4 時,

u = 0,1

又 Geven(u) 是 f(0),f(2),f(4),f(6) 中偶數編號 f(0),f(4) 的傅立葉轉換。

由演算法得知當 N = 2 時,直接可得

Geven(0) = [ f(0) ＋ f(4) ] / 2 for

Geven(1) = [ f(0) － f(4) ] / 2

同樣地，Godd(u) 是 f(0),f(2),f(4),f(6) 中奇數編號 f(2),f(6) 的傅立葉

轉換。

由演算法得知當 N = 2 時,

Godd(0) = [ f(2) ＋ f(6) ] / 2

Godd(1) = [ f(2) － f(6) ] / 2

將所算出的 Geven(0),Geven(1),Godd(0),Godd(1) 代回，可求得

G(0),G(1),G(2),G(3) ，也就是 Feven(0),Feven(1),Feven(2),Feven(3) 。

同理， Fodd(u) 是 f() 奇數編號 f(1),f(3),f(5),f(7) 的傅立葉轉換，

也可由上述方式推得。再將算出之 Feven(u) 及 Fodd(u) 代回，

即可算出 F(0),F(1),F(2),F(3),F(4),F(5),F(6),F(7)

6. 在 FFT 的計算中，須要將原資料做適當的排列以方便計算。

就上面的例子來說， f(0),f(1),f(2),f(3),f(4),f(5),f(6),f(7)

在第一次的分割分為偶數及奇數排列，變為

f(0),f(2),f(4),f(6) f(1),f(3),f(5),f(7)

第二次再將分出的部分再分偶數及奇數排列，變為

f(0),f(4) f(2),f(6) f(1),f(5) f(3),f(7)

如果能將原 f( ) 一開始就以 f(0),f(4),f(2),f(6),f(1),f(5),f(3),f(7)

的方式排列，即可依序運算，加快速度。

這種排列方式，如果以二進位的方式來看，好像由最高位元累加，

進位時向低位元進位。這樣的方式正好與一般二進位累加的方式相反，

稱為位元倒置(bit-reversal)，如下：

0	1	2	3	4	5	6	7
000	001	010	011	100	101	110	111
↓	↓	↓	↓	↓	↓	↓	↓
000	100	010	110	001	101	011	111
0	4	2	6	1	5	3	7

一般數位訊號處理器(DSP)為了加速 FFT 的運算，用硬體實作出位元

倒置，使 FFT 的運算更方便、快速。

四、程式架構及流程

本程式是以 Borland C++ Builder(以下簡稱 BCB) 為開發環境 ... ... ... ... ... ...

在實作 FFT 方面：

利用 2D-FT 的平移特性，先將原圖形資料乘上 (-1)^(x+y) 使 2D-FT 的原點移到圖形的中心點。

算出位元倒置的排列順序，以便 FFT 的運算。

利用 2D-FT 的可分離性，先對 Y 方向作 1D-FFT，再對 X 方向作1D-FFT，以簡化運算。

以線性壓縮的方式顯示 FFT 的結果，或利用動態範圍影像壓縮將 FFT 的結果壓縮，使高頻的資訊更明顯。另外，還有只對 X 方向或只對 Y 方向作 FFT 再以動態範圍影像壓縮顯示結果。

由於 BCB 已經建好了複數物件，所以在使用時，比一般用 C 語言所寫的要直覺，例如，在 C 語言中，可能利用一個結構來定義複數：

struct complex { double r, i; };

在作 FFT 時須要使用到複數的相乘以及相加(減)，所以，在 C 語言的 FFT 中，可見到如下的程式碼：

Fodd_Wn.r = F[o].r*Wn.r - F[o].i*Wn.i; /* Fodd * Wn */

Fodd_Wn.i = F[o].r*Wn.i + F[o].i*Wn.r;

F[o].r = F[e].r - Fodd_W.r; /* Feven - Fodd*Wn */

F[o].i = F[e].i - Fodd_W.i;

F[e].r = F[e].r + Fodd_W.r; /* Feven + Fodd*Wn */

F[e].i = F[e].i + Fodd_W.i;

以 BCB 物件的寫法，在使用上與實數沒什麼分別：

Fodd_Wn = F[o] * Wn; // Fodd * Wn

F[o] = F[e] - Fodd_Wn; // Feven - Fodd*Wn

F[o] = F[e] + Fodd_Wn; // Feven + Fodd*Wn

詳細的程式原始碼請見附錄。

五、結果

以下影像皆為 256×256 像素 256 灰階。由左至右分別為原始影像 f(x,y)、傅立葉頻譜 |F(u,v)| 經線性壓縮及非線生壓縮D(u,v)＝log[1+|F(u,v)|]，而原始影像的上方則標示其中灰階255(白色)部分的幾何特性。因應|F(u,v)|和D(u,v)的動態範圍，將其最大值對應至灰階255(白色)，其餘則依此比例調整。

圖一、背景灰階為 0 前景灰階為 255 之 8×8 正方形

其 |F(u,v)| 最大值（請見六、結論與討論）＝16320

動態範圍影像壓縮 c 值＝26.29 壓縮後 |F(u,v)| 最大值＝9.70

... ... ... ... ... ...

圖三、背景灰階為 0 前景灰階為 255 之 8×8 正方形旋轉 45°

其 |F(u,v)| 最大值＝12460

動態範圍影像壓縮 c 值＝27.04 壓縮後 |F(u,v)| 最大值＝9.43

... ... ... ... ... ...

Image Processing HomeWork #1(節錄)	Implementing 2D Discrete Fourier Transform

Dec.15,1998	Lee, Dong-Liang http://dllee.ktop.com.tw	5 / 6

2004-08-14T09:47:22