coreBugZJ

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

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

關(guān)于變量 X 的次數(shù)界為 n 多項(xiàng)式P(X)，

其系數(shù)表示法表示為

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

其點(diǎn)值表示法表示為

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

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

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

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

引理：對(duì)任何整數(shù) n>=0, k>=0, d>0, Wdn^dk = Wn^k 。

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

引理：如果 n>0 為偶數(shù)，n個(gè)n次單位復(fù)根的平方等于 n/2 個(gè) n/2 次單位復(fù)根。

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

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

由以上可以遞歸實(shí)現(xiàn) FFT 。

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

考察FFT遞歸的樹形結(jié)構(gòu)，可以將 A 中的元素按其在葉中出現(xiàn)的次序排序，然后從葉開始，一層層向根進(jìn)行迭代計(jì)算。

以上學(xué)習(xí)自《算法導(dǎo)論》，具體實(shí)現(xiàn)見代碼。

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

由其可分性，二維變換可以用二次一維變換實(shí)現(xiàn)，即先后在兩個(gè)“方向”上使用一維變換。

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

以下對(duì)于灰度圖來討論。

將圖像看做二維函數(shù)，圖像灰度值為函數(shù)在相應(yīng)XY處的函數(shù)值，對(duì)其進(jìn)行二維快速傅里葉變換，得到一個(gè)復(fù)數(shù)矩陣，將此矩陣水平循環(huán)移動(dòng)半寬，垂直循環(huán)移動(dòng)半高。新圖像的大小同此矩陣，而圖像中像素的灰度為復(fù)數(shù)矩陣中相應(yīng)位置的復(fù)數(shù)的模。新圖像還需進(jìn)行圖像增強(qiáng)。

逆變換，是在變換后得到的復(fù)數(shù)矩陣上進(jìn)行的二維傅里葉逆變換。

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

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/**//* 使用復(fù)數(shù) CF64 */
22 /**//* 由 base 開始，步長(zhǎng)為 step 的 (1<<lgn) 個(gè)元素按位翻轉(zhuǎn) */
23PublicFuncZ R32 bitReverseZ( CF64 *base, U32 step, U32 lgn );
24 /**//* 釋放變換當(dāng)中保存復(fù)數(shù)數(shù)據(jù)輸出所申請(qǐng)的內(nèi)存 */
25 /**//* 此函數(shù)只為避免申請(qǐng)與釋放內(nèi)存的方式不匹配 */
26PublicFuncZ R32 freeFftDataZ( CF64 **ppFftData );
27 /**//* 由 base 開始，步長(zhǎng)為 step 的 (1<<lgn) 個(gè)元素構(gòu)成的序列做 FFT */
28PublicFuncZ R32 doFftZ( CF64 *base, U32 step, U32 lgn );
29 /**//* 由 base 開始，步長(zhǎng)為 step 的 (1<<lgn) 個(gè)元素構(gòu)成的序列做 IFFT */
30PublicFuncZ R32 doIfftZ( CF64 *base, U32 step, U32 lgn );
31 /**//* 對(duì)寬為 (1<<lgw) 高為 (1<<lgh) 的矩陣做 2d FFT */
32 /**//* 點(diǎn)(x,y) 為 mat[ y * (1<<lgw) + x ] */
33 /**//* 先對(duì)每個(gè) (mat + y * (1<<lgw)) 做 FFT */
34PublicFuncZ R32 doFft2dZ( CF64 *mat, U32 lgw, U32 lgh );
35 /**//* 對(duì)寬為 (1<<lgw) 高為 (1<<lgh) 的矩陣做 2d IFFT */
36 /**//* 點(diǎn)的定位同 2d FFT */
37 /**//* 2d IFFT 的順序與 2d FFT 相反 */
38PublicFuncZ R32 doIfft2dZ( CF64 *mat, U32 lgw, U32 lgh );
39 /**//* 對(duì)圖像做 FFT，若圖像大小不合適，會(huì)自動(dòng)擴(kuò)充零，但保證不修改原圖 */
40 /**//* 若 NULL == pImgOut 則無(wú)圖像輸出，否則此參數(shù)輸出 FFT 后的圖像，會(huì)銷毀此參數(shù)原來的圖像 */
41 /**//* 若 NULL == ppFftOut 則無(wú)數(shù)據(jù)輸出，否則此參數(shù)輸出 FFT 后的復(fù)數(shù)數(shù)據(jù)，會(huì)釋放此參數(shù)原來的內(nèi)存 */
42 /**//* 若 NULL == ppFftOut 則忽略參數(shù) plgw, plgh，否則此參數(shù)輸出 FFT 后的復(fù)數(shù)數(shù)據(jù)的規(guī)模 */
43 /**//* 輸出時(shí)，(*ppFftOut) 中的數(shù)據(jù)格式同 2d FFT */
44PublicFuncZ R32 doImageFftZ( cImageZ src, ImageZ *pImgOut, CF64 **ppFftOut, U32 *plgw, U32 *plgh );
45 /**//* 對(duì) image FFT 后的數(shù)據(jù) mat 做 2d IFFT 得到原始圖像 */
46 /**//* 及原始圖像 image FFT 前的數(shù)據(jù) */
47 /**//* 會(huì)銷毀原來的 (*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/**//* 使用復(fù)數(shù) 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 閱讀(11604) 評(píng)論(4)  編輯 收藏 引用 所屬分類: VideoImage 、Algorithm 、Mathematics 、課內(nèi)作業(yè)