coreBugZJ

此 blog 已弃。

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

Feedback

# re: 数字图像处理上机之四:灰度图 快速傅里叶变换 ( FFT IFFT 一维 二维 )[未登录] 2012-02-28 12:51 U

你好!那个编译错误显示找不到TypeZ.h是为什么??新手上路。。这两天急用。。望赐教,谢谢~~!!  回复  更多评论   

# re: 数字图像处理上机之四:灰度图 快速傅里叶变换 ( FFT IFFT 一维 二维 ) 2012-02-29 12:15 coreBugZJ

@U
数字图像处理上机是一个系列的,你在其它文件里找找。  回复  更多评论   

# re: 数字图像处理上机之四:灰度图 快速傅里叶变换 ( FFT IFFT 一维 二维 ) 2012-03-16 10:21 C小加

radon变换你会么、?  回复  更多评论   

# re: 数字图像处理上机之四:灰度图 快速傅里叶变换 ( FFT IFFT 一维 二维 ) 2012-03-16 19:44 coreBugZJ

@C小加
不会呀,求指教。。。  回复  更多评论   



只有注册用户登录后才能发表评论。
网站导航: 博客园   IT新闻   BlogJava   知识库   博问   管理