數字圖像處理上機之四：灰度圖 快速傅里葉變換 （ FFT IFFT 一維 二維 ）

1. 一維快速傅里葉變換的原理：

關于變量 X 的次數界為 n 多項式P(X)，

其系數表示法表示為

        P(X) = A0 * X^0 + A1 * X^1 + ... + An-1 * X^(n-1)

其點值表示法表示為

        n 個點值對組成的集合 { (X0,Y0), (X1,Y1), ..., (Xn-1,Yn-1) }，
        集合中所有 Xi 各不相同且 Yi = P(Xi)。
        顯然，點值表示不唯一。

定理：對于任意 n 個點值對組成的集合 { (X0,Y0), (X1,Y1), ..., (Xn-1,Yn-1) }，存在唯一的次數界為 n 的多項式 P(X)，滿足 Yi = P(Xi)，i = 0, 1, ... n-1 。

精心挑選 n 個點，可以使兩種表示相互轉化的算法的時間復雜度壓縮為 nlog(n)。
如果選擇“單位復根”作為求值點，則可以通過對系數向量做離散傅里葉變換（DFT），得到相應的點值表示，也可以通過對點值執行“逆DFT”運算，而獲得相應的系數向量。

n 次單位復根是滿足 W^n = 1 的復數 W ，n 次單位復根剛好有 n 個，它們是 e^(2*PI*i*k / n)，k = 0, 1, ..., n-1 。
Wn = e^(2*PI*i/n) 稱為主n次單位根，其它n次單位根都是它的冪。

引理：對任何整數 n>=0, k>=0, d>0, Wdn^dk = Wn^k 。

推論：對任意偶數 n>0, 有 Wn^(n/2) = W2 = -1 。

引理：如果 n>0 為偶數，n個n次單位復根的平方等于 n/2 個 n/2 次單位復根。

假定 n 為2的冪，則次數界為 n 的多項式 P(X) = A0 * X^0 + A1 * X^1 + ... + An-1 * X^(n-1) ，其系數向量為 A = (A0, A1, A2, ... An-1)，P(X)在 n 個單位復根處的值為 Yk = P(Wn^k)，向量 Y = (Y0, Y1, ... , Yn-1) 是系數向量 A 的離散傅里葉變換（DFT），寫作 Y = DFTn(A) 。

使用快速傅里葉變換(FFT)的方法，可以在 nlog(n) 的時間內計算出 DFTn(A)，主要利用了單位復根的性質。
FFT 用了分治策略，對 P(X) 定義兩個多項式
        P0(X) = A0 + A2 * X + A4 * X^2 + ... + An-2 * X^(n/2-1)    偶數下標
        P1(X) = A1 + A3 * X + A5 * X^2 + ... + An-1 * X^(n/2-1)    奇數下標
則
        P(X)  = P0(X^2) + X * P1(X^2)

由以上可以遞歸實現 FFT 。

使用FFT方法求“逆DFTn”，做法與以上類似，只是把 A 與 Y 的角色互換，用 Wn^(-1) 代替 Wn，并且將每個結果元素除以 n 。

考察FFT遞歸的樹形結構，可以將 A 中的元素按其在葉中出現的次序排序，然后從葉開始，一層層向根進行迭代計算。

以上學習自《算法導論》，具體實現見代碼。

2. 二維快速傅里葉變換的原理：

由其可分性，二維變換可以用二次一維變換實現，即先后在兩個“方向”上使用一維變換。

3. 數字圖像中的快速傅里葉變換：

以下對于灰度圖來討論。

將圖像看做二維函數，圖像灰度值為函數在相應XY處的函數值，對其進行二維快速傅里葉變換，得到一個復數矩陣，將此矩陣水平循環移動半寬，垂直循環移動半高。新圖像的大小同此矩陣，而圖像中像素的灰度為復數矩陣中相應位置的復數的模。新圖像還需進行圖像增強。

逆變換，是在變換后得到的復數矩陣上進行的二維傅里葉逆變換。

4. 注意：
一維快速傅里葉變換需要在長度為 2的冪 的序列上進行，若長度不足 2的冪，需在高位補零，湊足長度。其它類推。

1/**//*
2ProcessFftZ.h
3
4Copyright (C) 2011, coreBugZJ, all rights reserved.
5
6快速傅里葉變換FFT。
7*/
8
9
10#ifndef __PROCESSFFT_Z_H_INCLUDED__
11#define __PROCESSFFT_Z_H_INCLUDED__
12
13
14#include "TypeZ.h"
15#include "ClassImageZ.h"
16
17
18 /**//* 得到使 (1<<lgn) >= n 的最小的 lgn */
19PublicFuncZ R32 getUplgnZ( U32 *plgn, U32 n );
20
21/**//* 使用復數 CF64 */
22 /**//* 由 base 開始，步長為 step 的 (1<<lgn) 個元素按位翻轉 */
23PublicFuncZ R32 bitReverseZ( CF64 *base, U32 step, U32 lgn );
24 /**//* 釋放變換當中保存復數數據輸出所申請的內存 */
25 /**//* 此函數只為避免申請與釋放內存的方式不匹配 */
26PublicFuncZ R32 freeFftDataZ( CF64 **ppFftData );
27 /**//* 由 base 開始，步長為 step 的 (1<<lgn) 個元素構成的序列做 FFT */
28PublicFuncZ R32 doFftZ( CF64 *base, U32 step, U32 lgn );
29 /**//* 由 base 開始，步長為 step 的 (1<<lgn) 個元素構成的序列做 IFFT */
30PublicFuncZ R32 doIfftZ( CF64 *base, U32 step, U32 lgn );
31 /**//* 對寬為 (1<<lgw) 高為 (1<<lgh) 的矩陣做 2d FFT */
32 /**//* 點(x,y) 為 mat[ y * (1<<lgw) + x ] */
33 /**//* 先對每個 (mat + y * (1<<lgw)) 做 FFT */
34PublicFuncZ R32 doFft2dZ( CF64 *mat, U32 lgw, U32 lgh );
35 /**//* 對寬為 (1<<lgw) 高為 (1<<lgh) 的矩陣做 2d IFFT */
36 /**//* 點的定位同 2d FFT */
37 /**//* 2d IFFT 的順序與 2d FFT 相反 */
38PublicFuncZ R32 doIfft2dZ( CF64 *mat, U32 lgw, U32 lgh );
39 /**//* 對圖像做 FFT，若圖像大小不合適，會自動擴充零，但保證不修改原圖 */
40 /**//* 若 NULL == pImgOut 則無圖像輸出，否則此參數輸出 FFT 后的圖像，會銷毀此參數原來的圖像 */
41 /**//* 若 NULL == ppFftOut 則無數據輸出，否則此參數輸出 FFT 后的復數數據，會釋放此參數原來的內存 */
42 /**//* 若 NULL == ppFftOut 則忽略參數 plgw, plgh，否則此參數輸出 FFT 后的復數數據的規模 */
43 /**//* 輸出時，(*ppFftOut) 中的數據格式同 2d FFT */
44PublicFuncZ R32 doImageFftZ( cImageZ src, ImageZ *pImgOut, CF64 **ppFftOut, U32 *plgw, U32 *plgh );
45 /**//* 對 image FFT 后的數據 mat 做 2d IFFT 得到原始圖像 */
46 /**//* 及原始圖像 image FFT 前的數據 */
47 /**//* 會銷毀原來的 (*pImgOut) 圖像 */
48PublicFuncZ R32 doImageIfftZ( ImageZ *pImgOut, CF64 *mat, U32 lgw, U32 lgh );
49
50
51#endif /* __PROCESSFFT_Z_H_INCLUDED__ */
52

1/**//*
2ProcessFftZ.c
3
4Copyright (C) 2011, coreBugZJ, all rights reserved.
5
6快速傅里葉變換FFT。
7*/
8
9
10#include "stdafx.h"
11#include "ProcessFftZ.h"
12#include "MathZ.h"
13#include "ProcessGrayZ.h"
14
15#include <malloc.h>
16
17
18PublicFuncZ R32 getUplgnZ( U32 *plgn, U32 n ) {
19 if ( (NULL == plgn) || (0x80000000 < n) ) {
20 return RERR;
21 }
22 *plgn = 0;
23 while ( PWR2_U32_Z( *plgn ) < n ) {
24 ++(*plgn);
25 }
26 return ROK;
27}
28
29
30/**//* 使用復數 CF64 */
31
32PublicFuncZ R32 bitReverseZ( CF64 *base, U32 step, U32 lgn ) {
33 U32 i, j, k, v, n;
34 CF64 tmp;
35
36 if ( NULL == base ) {
37 return RERR;
38 }
39
40 n = PWR2_U32_Z( lgn );
41 for ( i = 0; i < n; ++i ) {
42 k = i;
43 j = 0;
44 for ( v = 0; v < lgn; ++v ) {
45 j = ( (j<<1) | (k&1) );
46 k >>= 1;
47 }
48 if ( i < j ) {
49 MOV_CF64_Z( tmp, base[ step * i ] );
50 MOV_CF64_Z( base[ step * i ], base[ step * j ] );
51 MOV_CF64_Z( base[ step * j ], tmp );
52 }
53 }
54 return ROK;
55}
56
57PublicFuncZ R32 freeFftDataZ( CF64 **ppFftData ) {
58 if ( NULL == ppFftData ) {
59 return RERR;
60 }
61
62 free( *ppFftData );
63 *ppFftData = NULL;
64 return ROK;
65}
66
67PublicFuncZ R32 doFftZ( CF64 *base, U32 step, U32 lgn ) {
68 U32 n, lgm, m, i, j;
69 CF64 w, wn, u, t;
70 F64 recos, imsin, tf;
71
72 if ( NULL == base ) {
73 return RERR;
74 }
75
76 bitReverseZ( base, step, lgn );
77 n = PWR2_U32_Z( lgn );
78 for ( lgm = 0; lgm < lgn; ++lgm ) {
79 m = PWR2_U32_Z( lgm );
80 DIV_F64_U32_Z( tf, PI_F64_Z, m );
81 COS_F64_Z( recos, tf );
82 SIN_F64_Z( imsin, tf );
83 MOV_CF64_F64_Z( wn, recos, imsin );
84 for ( i = 0; i < n; i += m ) {
85 MOV_CF64_U32_Z( w, 1, 0 );
86 for ( j = i+m; i < j; ++i ) {
87 MOV_CF64_Z( u, base[ step * i ] );
88 MUL_CF64_Z( t, w, base[ step * (i + m) ] );
89 ADD_CF64_Z( base[ step * i ], u, t );
90 SUB_CF64_Z( base[ step * (i + m) ], u, t );
91 MUL_CF64_Z( w, w, wn );
92 }
93 }
94 }
95
96 return ROK;
97}
98
99PublicFuncZ R32 doIfftZ( CF64 *base, U32 step, U32 lgn ) {
100 U32 n, lgm, m, i, j;
101 CF64 w, wn, u, t;
102 F64 recos, imsin, tf;
103
104 if ( NULL == base ) {
105 return RERR;
106 }
107
108 bitReverseZ( base, step, lgn );
109 n = PWR2_U32_Z( lgn );
110 for ( lgm = 0; lgm < lgn; ++lgm ) {
111 m = PWR2_U32_Z( lgm );
112 DIV_F64_U32_Z( tf, PI_F64_Z, m );
113 COS_F64_Z( recos, tf );
114 SIN_F64_Z( imsin, tf );
115 NEG_F64_Z( imsin, imsin );
116 MOV_CF64_F64_Z( wn, recos, imsin );
117 for ( i = 0; i < n; i += m ) {
118 MOV_CF64_U32_Z( w, 1, 0 );
119 for ( j = i+m; i < j; ++i ) {
120 MOV_CF64_Z( u, base[ step * i ] );
121 MUL_CF64_Z( t, w, base[ step * (i + m) ] );
122 ADD_CF64_Z( base[ step * i ], u, t );
123 SUB_CF64_Z( base[ step * (i + m) ], u, t );
124 MUL_CF64_Z( w, w, wn );
125 }
126 }
127 }
128 MOV_CF64_U32_Z( t, n, 0 );
129 for ( i = 0; i < n; ++i ) {
130 DIV_CF64_Z( base[ step * i ], base[ step * i ], t );
131 }
132
133 return ROK;
134}
135
136PublicFuncZ R32 doFft2dZ( CF64 *mat, U32 lgw, U32 lgh ) {
137 U32 x, y, w, h;
138
139 if ( NULL == mat ) {
140 return RERR;
141 }
142
143 w = PWR2_U32_Z( lgw );
144 h = PWR2_U32_Z( lgh );
145
146 for ( y = 0; y < h; ++y ) {
147 doFftZ( mat + y * w, 1, lgw );
148 }
149 for ( x = 0; x < w; ++x ) {
150 doFftZ( mat + x, w, lgh );
151 }
152
153 return ROK;
154}
155
156PublicFuncZ R32 doIfft2dZ( CF64 *mat, U32 lgw, U32 lgh ) {
157 U32 x, y, w, h;
158
159 if ( NULL == mat ) {
160 return RERR;
161 }
162
163 w = PWR2_U32_Z( lgw );
164 h = PWR2_U32_Z( lgh );
165
166 for ( x = 0; x < w; ++x ) {
167 doIfftZ( mat + x, w, lgh );
168 }
169 for ( y = 0; y < h; ++y ) {
170 doIfftZ( mat + y * w, 1, lgw );
171 }
172
173 return ROK;
174}
175
176PublicFuncZ R32 doImageFftZ( cImageZ src, ImageZ *pImgOut, CF64 **ppFftOut, U32 *plgw, U32 *plgh ) {
177 ImageZ img = NULL;
178 CF64 *mat = NULL;
179 B32 outImg = (NULL != pImgOut), outFft = (NULL != ppFftOut);
180 U32 width, height, lgw, w, lgh, h;
181 U32 x, y, i;
182 U32 tu;
183 U08 *ptrimg;
184 CF64 *ptrmat;
185 F64 tf, tf255;
186
187 if ( (! isImageValidZ(src)) || (GRAY_NUM_Z != src->colorNum) ) {
188 return RERR;
189 }
190 if ( (! outImg) && (! outFft) ) {
191 return ROK;
192 }
193 if ( outFft && ( (NULL == plgw) || (NULL == plgh) ) ) {
194 return RERR;
195 }
196
197 if ( outFft && (NULL != (*ppFftOut)) ) {
198 free( *ppFftOut );
199 *ppFftOut = NULL;
200 }
201
202 width = src->width;
203 height = src->height;
204 getUplgnZ( &lgw, width );
205 getUplgnZ( &lgh, height );
206 w = PWR2_U32_Z( lgw );
207 h = PWR2_U32_Z( lgh );
208
209 mat = (CF64*)malloc( sizeof(CF64) * w * h );
210 if ( NULL == mat ) {
211 return RERR;
212 }
213
214 for ( y = 0; y < height; ++y ) {
215 ptrimg = src->pPixel + y * src->linePitch;
216 ptrmat = mat + y * w;
217 for ( x = 0; x < width; ++x ) {
218 tu = *ptrimg++;
219 MOV_CF64_U32_Z( *ptrmat, tu, 0 );
220 ++ptrmat;
221 }
222 for ( x = width; x < w; ++x ) {
223 MOV_CF64_U32_Z( *ptrmat, 0, 0 );
224 ++ptrmat;
225 }
226 }
227 for ( y = height; y < h; ++y ) {
228 ptrmat = mat + y * w;
229 for ( x = 0; x < w; ++x ) {
230 MOV_CF64_U32_Z( *ptrmat, 0, 0 );
231 ++ptrmat;
232 }
233 }
234
235 if ( outImg && (NULL != (*pImgOut)) ) {
236 destroyImageZ( *pImgOut );
237 *pImgOut = NULL;
238 }
239
240 doFft2dZ( mat, lgw, lgh );
241
242 while ( outImg ) {
243 img = createImageZ( w, h, GRAY_NUM_Z );
244 if ( NULL == img ) {
245 break;
246 }
247 ptrimg = img->pPalette;
248 for ( i = 0; i < GRAY_NUM_Z; ++i ) {
249 ptrimg[ IMAGEZ_OFFSET_BLUE_Z ] =
250 ptrimg[ IMAGEZ_OFFSET_GREEN_Z ] =
251 ptrimg[ IMAGEZ_OFFSET_RED_Z ] = (U08)(i);
252 ptrimg[ IMAGEZ_OFFSET_ALPHA_Z ] = 0;
253 ptrimg += IMAGEZ_COLOR_SIZE_Z;
254 }
255 MOV_F64_U32_Z( tf255, 255 );
256 for ( y = 0; y < h; ++y ) {
257 ptrimg = img->pPixel + ((y+(h>>1))%h) * img->linePitch;
258 ptrmat = mat + y * w;
259 for ( x = 0; x < w; ++x ) {
260 M_CF64_Z( tf, *ptrmat );
261 ++ptrmat;
262 DIV_F64_U32_Z( tf, tf, 100 );
263 MIN_F64_Z( tf, tf, tf255 );
264 MOV_U32_F64_Z( tu, tf );
265 *(ptrimg + (x+(w>>1))%w ) = (U08)(tu);
266 }
267 }
268 break;
269 }
270 if ( outImg ) {
271 *pImgOut = img;
272 }
273
274 if ( outFft ) {
275 *ppFftOut = mat;
276 *plgw = lgw;
277 *plgh = lgh;
278 }
279 else {
280 free( mat );
281 mat = NULL;
282 }
283
284 return ROK;
285}
286
287PublicFuncZ R32 doImageIfftZ( ImageZ *pImgOut, CF64 *mat, U32 lgw, U32 lgh ) {
288 U32 x, y, w, h, i;
289 U32 re, im;
290 CF64 *ptrmat;
291 U08 *ptrimg;
292
293 if ( (NULL == pImgOut) || (NULL == mat) ) {
294 return RERR;
295 }
296 if ( NULL != (*pImgOut) ) {
297 destroyImageZ( *pImgOut );
298 *pImgOut = NULL;
299 }
300
301 w = PWR2_U32_Z( lgw );
302 h = PWR2_U32_Z( lgh );
303
304 *pImgOut = createImageZ( w, h, GRAY_NUM_Z );
305 if ( NULL == (*pImgOut) ) {
306 return RERR;
307 }
308
309 ptrimg = (*pImgOut)->pPalette;
310 for ( i = 0; i < GRAY_NUM_Z; ++i ) {
311 ptrimg[ IMAGEZ_OFFSET_BLUE_Z ] =
312 ptrimg[ IMAGEZ_OFFSET_GREEN_Z ] =
313 ptrimg[ IMAGEZ_OFFSET_RED_Z ] = (U08)(i);
314 ptrimg[ IMAGEZ_OFFSET_ALPHA_Z ] = 0;
315 ptrimg += IMAGEZ_COLOR_SIZE_Z;
316 }
317
318 doIfft2dZ( mat, lgw, lgh );
319
320 for ( y = 0; y < h; ++y ) {
321 ptrmat = mat + y * w;
322 ptrimg = (*pImgOut)->pPixel + y * (*pImgOut)->linePitch;
323 for ( x = 0; x < w; ++x ) {
324 MOV_U32_CF64_Z( re, im, ptrmat[ x ] );
325 *ptrimg++ = (U08)(re);
326 }
327 }
328
329 return ROK;
330}
331

posted on 2011-11-25 23:03 coreBugZJ 閱讀(11663) 評論(4)  編輯 收藏 引用 所屬分類: VideoImage 、Algorithm 、Mathematics 、課內作業