coreBugZJ

数字图像处理上机之四：灰度图 快速傅里叶变换 （ 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 阅读(11589) 评论(4)  编辑 收藏 引用 所属分类: VideoImage 、Algorithm 、Mathematics 、课内作业