vc++实现的傅立叶变换，参考算法导论中的快速傅立叶算法，使用openCV做图像的读入与显示，加入了高斯函数低通滤波；
1 /************************************************************ 2 *二维离散傅立叶变换和滤波函数（高斯低通）的实现 3 *用于计算图像离散序列的傅立叶频谱图，并用于频域图像处理 4 *算法：多项式快速傅立叶算法（蝶形） 理论基础：算法导论,图像处理 5 *时间复杂度：一维O(NlogN)，二维O(N^2) 6 *工具：openCV读取图像与展示效果 7 *版本：测试 8 *@auto:Bruce mu 9 ************************************************************/ 10 #include <stdio.h> 11 #include <iostream> 12 #include <complex> 13 #include <bitset> 14 #include <fstream> 15 #include <algorithm> 16 #include <iterator> 17 #include <math.h> 18 #include "cv.h" 19 #include "highgui.h" 20 #include "CImg.h" 21 #define pi 3.1415926535 22 #define DELTA 5.0 23 using std::iostream; 24 using std::bitset; 25 using std::complex; 26 using namespace std; 27 using namespace cimg_library; 28 29 /*******为自底向上的迭代重排序列计算位置************/ 30 int rev(int k,int n) 31 { 32 bitset<32> tmp(k); 33 bitset<32> dist; 34 for(int m = 0;m<n;m++) 35 { 36 if(tmp.test(m)) 37 { 38 dist.set(n-1-m); 39 } 40 } 41 int revk = (int)dist.to_ulong(); 42 return revk; 43 } 44 //重载形式 45 complex<double>* FFT(const complex<double> *srcimg,int n) 46 { 47 double flag = log10((double)n)/log10(2.0); 48 int N = n; 49 if(flag - (int)flag != 0){ //判断是否为2的指数次 50 cout <<"the length of srcimg is wrong"<< endl; 51 /*填充成2的指数项*/ 52 cout <<"need padding"<<endl; 53 N = pow(2,(double)((int)flag+1)); 54 flag = log10((double)N)/log10(2.0); 55 } 56 /*改变成奇偶顺序*/ 57 complex<double> *arr= new complex<double>[N]; 58 int sub; 59 for(int k =0;k<N;k++) 60 { 61 sub =rev(k,(int)flag); 62 if(sub <= n-1){ 63 arr[k] = *(srcimg + sub); 64 }else{ 65 complex<double> t = complex<double>(0,0); 66 arr[k] = t; 67 } 68 } 69 for(int s =1;s <= flag;s++) 70 { 71 int m = pow(2,(double)s); 72 complex<double> wm = complex<double>(cos(2*pi/m),sin(2*pi/m));//wm1：从W21开始,周期变换 73 for(int p = 0;p<N;p+=m) 74 { 75 complex<double> w(1,0); 76 for(int j = 0;j<=m/2-1;j++) 77 { 78 complex<double> t = w * arr[p+j+m/2]; 79 complex<double> u = arr[p+j]; 80 arr[p+j] = u+t; 81 arr[p+j+m/2] = u-t; 82 w = w*wm; 83 } 84 } 85 } 86 return arr; 87 } 88 /***********一维快速傅立叶变换******************** 89 *srcimg : 原始一维序列 * 90 *n ：一维序列的长度 91 *************************************************/ 92 complex<double>* FFT(const double *srcimg,int n) 93 { 94 double flag = log10((double)n)/log10(2.0); 95 int N = n; 96 if(flag - (int)flag != 0){ //判断是否为2的指数次 97 // cout <<"the length of srcimg is wrong"<< endl; 98 /*填充成2的指数项*/ 99 // cout <<"need padding"<<endl; 100 N = pow(2,(double)((int)flag+1)); 101 flag = log10((double)N)/log10(2.0); 102 } 103 /*改变成奇偶顺序*/ 104 complex<double> *arr= new complex<double>[N]; 105 int sub; 106 for(int k =0;k<N;k++) 107 { 108 sub =rev(k,(int)flag); 109 if(sub <= n-1){ 110 arr[k] = complex<double>(*(srcimg + sub),0); 111 }else{ 112 complex<double> t = complex<double>(0,0); 113 arr[k] = t; 114 } 115 } 116 // cout<<"------------after padding and retrival-----------------"<<endl; 117 // for(size_t y=0;y<N;y++) 118 // { 119 // cout << arr[y].real()<<" "<<arr[y].imag()<<endl; 120 // } 121 /*基于迭代的蝶形快速傅立叶变换，自底向上*/ 122 for(int s =1;s <= flag;s++) 123 { 124 int m = pow(2,(double)s); 125 complex<double> wm = complex<double>(cos(2*pi/m),sin(2*pi/m));//wm1：从W21开始,周期变换 126 for(int p = 0;p<N;p+=m) 127 { 128 complex<double> w(1,0); 129 for(int j = 0;j<=m/2-1;j++) 130 { 131 complex<double> t = w * arr[p+j+m/2]; 132 complex<double> u = arr[p+j]; 133 arr[p+j] = u+t; 134 arr[p+j+m/2] = u-t; 135 w = w*wm; 136 } 137 } 138 } 139 return arr; 140 } 141 142 int countPadding(int n); 143 /*****************一维快速傅立叶逆变换********************/ 144 /*fftimg:原始一维傅立叶序列 145 n : 序列长度 146 *******************************************************/ 147 complex<double>* IFFT(const complex<double> *fftimg,int n) 148 { 149 n = countPadding(n); 150 double flag = log10((double)n)/log10(2.0); 151 int N = n; 152 if(flag - (int)flag != 0){ //判断是否为2的指数次 153 cout <<"the length of srcimg is wrong"<< endl; 154 /*填充成2的指数项*/ 155 cout <<"need padding"<<endl; 156 N = pow(2,(double)((int)flag+1)); 157 flag = log10((double)N)/log10(2.0); 158 } 159 /*改变成奇偶顺序*/ 160 complex<double> * spit = new complex<double>[N]; 161 int sub=0; 162 for(int k =0;k<N;k++) 163 { 164 sub =rev(k,(int)flag); 165 if(sub < n){ 166 spit[k] = complex<double>(*(fftimg + sub)); 167 }else{ 168 spit[k] = complex<double>(0,0); 169 } 170 } 171 172 for(int s =1;s <= flag;s++) 173 { 174 int m = pow(2,(double)s); 175 complex<double> wm = complex<double>(cos(-2*pi/m),sin(-2*pi/m));//wm1：从W2(-1)开始 176 for(int p = 0;p<N;p+=m) 177 { 178 complex<double> w(1,0); 179 for(int j = 0;j<=m/2-1;j++) 180 { 181 complex<double> t = w * spit[p+j+m/2]; 182 complex<double> u = spit[p+j]; 183 spit[p+j] = u+t; 184 spit[p+j+m/2] = u-t; 185 w = w*wm; 186 } 187 } 188 } 189 190 for(size_t p =0;p<n;p++) 191 { 192 spit[p] = spit[p]/complex<double>(N,0); 193 } 194 return spit; 195 } 196 197 /*******使用共轭的快速傅立叶逆变换**************/ 198 complex<double>* IFFT2(const complex<double> *fftimg,int n) 199 { 200 n = countPadding(n); 201 complex<double> *gfftimg = new complex<double>[n]; 202 for(size_t i = 0;i<n;i++){ 203 gfftimg[i] = complex<double>(fftimg[i].real(),-fftimg[i].imag()); 204 } 205 complex<double> *ifft = FFT(gfftimg,n); 206 for(size_t j = 0;j<n;j++) 207 { 208 ifft[j] = ifft[j]/complex<double>(n,0); 209 } 210 delete gfftimg; 211 return ifft; 212 } 213 214 /*************二维快速傅立叶变换************************** 215 *srcimg: 用一维表示的二维原始序列 216 *width ：宽度 217 *height：高度 218 ********************************************************/ 219 complex<double>* twoDFFT(const double *srcimg,int width,int height) 220 { 221 int w = countPadding(width); 222 int h = countPadding(height); 223 int pixes = w * h; 224 complex<double> *hdirection = new complex<double>[w]; 225 complex<double> *vdirection = new complex<double>[h]; 226 complex<double> *fourier = new complex<double>[pixes]; 227 /*二维水平方向*/ 228 for(size_t i = 0;i<h;i++){ 229 for(size_t j = 0;j<w;j++){ 230 if(i>=height || j >=width){ 231 hdirection[j] = complex<double>(0,0); 232 }else{ 233 hdirection[j] = complex<double>(srcimg[i*width + j],0); 234 } 235 // cout << hdirection[j] << " "; 236 } 237 // cout <<""<<endl; 238 complex<double> *hfourier = FFT(hdirection,w); 239 for(size_t m = 0;m<w;m++){ 240 fourier[i*w+m] = hfourier[m]; 241 } 242 delete hfourier; 243 } 244 /*二维垂直方向*/ 245 for(size_t ii = 0;ii<w;ii++){ 246 for(size_t jj = 0;jj<h;jj++){ 247 vdirection[jj] = fourier[jj*w + ii]; 248 } 249 complex<double> *vfourier = FFT(vdirection,h); 250 for(size_t mm = 0;mm < h;mm++){ 251 fourier[mm*w +ii] = vfourier[mm]; 252 } 253 delete vfourier; 254 } 255 delete hdirection; 256 delete vdirection; 257 return fourier; 258 259 } 260 /**************二维快速傅立叶逆变换************************* 261 *fourier : 一维表示的二维傅立叶变换序列 * 262 *width :宽 * 263 *height :高 * 264 ***********************************************************/ 265 complex<double>* twoDIFFT(const complex<double> *fourier,int width,int height) 266 { 267 width = countPadding(width); 268 height = countPadding(height); 269 int fpoints = width * height; 270 complex<double> *hdirection = new complex<double>[width]; 271 complex<double> *vdirection = new complex<double>[height]; 272 complex<double> *ifourier = new complex<double>[fpoints]; 273 274 for(size_t ii = 0;ii<height;ii++) 275 { 276 for(size_t jj = 0;jj<width;jj++){ 277 hdirection[jj] = fourier[ii*width+jj]; 278 } 279 complex<double> *hifour = IFFT(hdirection,width);//临时变量 280 for(size_t mm = 0;mm<width;mm++){ 281 ifourier[ii*width+mm] = hifour[mm]; 282 } 283 delete hifour; 284 } 285 for(size_t i = 0;i<width;i++){ 286 for(size_t j = 0;j<height;j++){ 287 vdirection[j] = ifourier[j*width+i]; 288 } 289 complex<double> *vifour = IFFT(vdirection,height); 290 for(size_t m = 0;m<height;m++){ 291 ifourier[m*width+i] = vifour[m]; 292 } 293 delete vifour; 294 } 295 delete hdirection; 296 delete vdirection; 297 return ifourier; 298 } 299 /******************计算填充数******************************************** 300 *蝶形傅立叶变换算法只计算2的指数次的离散序列，对于非2的指数次的序列，使用零填充* 301 ***********************************************************************/ 302 inline int countPadding(int n) 303 { 304 double lg = log10((double)n)/log10(2.0); 305 if((lg - (int)lg) == 0){ 306 return n; 307 } 308 int N = pow(2.0,((int)lg+1)); 309 return N; 310 } 311 312 /*****高斯低通滤波函数************************* 313 *src： 输入频谱 314 *width： 输入频谱宽度 315 *height：输入频谱高度 316 *D： 高斯函数方差，即滤波阈值 317 *只对于移位居中后的傅立叶频谱进行高斯低通滤波 318 */ 319 void guass(complex<double> *src,int width,int height,double D) 320 { 321 //find centre point 322 int orgx = floor(width/2.0); 323 int orgy = floor(height/2.0); 324 double distence = 0; 325 for(size_t i = 0;i<height;i++) 326 { 327 for(size_t j = 0;j<width;j++) 328 { 329 distence = sqrt(pow(abs((int)(i-orgy)),2.0)+pow(abs((int)(j-orgx)),2.0)); 330 src[i*width+j] = src[i*width+j] * exp(-distence/(2*pow(D,2.0))); 331 332 } 333 //cout<<i<<" "<<distence<<" "<<endl; 334 335 } 336 // cout<<orgx<<" "<<orgy<<endl; 337 } 338 /************复数傅立叶频谱数组转换成256级灰度数组*****************/ 339 void toGray(complex<double> *twodfourier,int pwid,int phei,char* pixval) 340 { 341 double *vals = new double[pwid*phei]; 342 double max = 0; 343 double min = 255; 344 for(size_t p = 0;p<pwid*phei;p++){ 345 vals[p]=log(1+sqrt(pow(twodfourier[p].real(),2.0)+pow(twodfourier[p].imag(),2.0)));//对数级的幅度铺 346 if(vals[p] > max){ 347 max = vals[p]; 348 } 349 if(vals[p] < min){ 350 min = vals[p]; 351 } 352 } 353 cout<<min<< " "<<max<<endl; 354 cout<<pwid<<" "<<phei<<endl; 355 for(size_t s = 0;s<pwid*phei;s++){ 356 pixval[s] = (char)((vals[s]-min)/(max-min)*255); 357 } 358 delete vals; 359 } 360 /******************opencv 读取图像与展示效果***********************/ 361 int main(int argc,char **argv) 362 { 363 IplImage *img; 364 if((img = cvLoadImage("F:\\Fig0222(a)(face).tif",0))!=0){ 365 int dim = img->imageSize; 366 //从图像矩阵复制出单维数组，并做频谱居中操作，对应改写接口，返回单维数组； 367 double * src = new double[dim]; 368 size_t j =0; 369 for(int y =0;y<img->height;y++) 370 { 371 uchar * ptr = (uchar*)(img->imageData + y * img->widthStep); 372 for(int x =0;x<img->width;x++){ 373 /***输入函数乘以（-1）的（x+y）次方，频谱将居中*/ 374 src[j] = (double)ptr[x]*pow(-1,(double)(x+y))/256; 375 j++; 376 } 377 } 378 int w = img->width; 379 int h = img->height; 380 int pwid = countPadding(w); 381 int phei = countPadding(h); 382 383 complex<double> *twodfourier = twoDFFT(src,w,h); 384 char * pixval = new char[pwid*phei]; 385 CvMat freq; 386 toGray(twodfourier,pwid,phei,pixval);//复数频谱转256灰度 387 cvInitMatHeader(&freq,pwid,phei,CV_8UC1,pixval); 388 IplImage *ipl = cvCreateImage(cvGetSize(&freq),8,1); 389 cvGetImage(&freq,ipl);//获取频谱图像 390 391 guass(twodfourier,pwid,phei,DELTA);//方差DELTA的高斯平滑（高斯低通滤波） 392 CvMat gaufre; 393 char *pixvals = new char[pwid*phei]; 394 toGray(twodfourier,pwid,phei,pixvals); 395 cvInitMatHeader(&gaufre,pwid,phei,CV_8UC1,pixvals); 396 IplImage *gausf = cvCreateImage(cvGetSize(&gaufre),8,1); 397 cvGetImage(&gaufre,gausf);//高斯平滑后的频谱图像 398 399 complex<double> *twodifourier = twoDIFFT(twodfourier,w,h); 400 double ipix = 0; 401 size_t po = 0; 402 for(int y =0;y<img->height;y++) 403 { 404 uchar * ptr = (uchar*)(img->imageData + y * img->widthStep); 405 for(int x =0;x<img->width;x++){ 406 ipix = twodifourier[po*pwid+x].real(); 407 ptr[x] = (uchar)(ipix * 256)*pow(-1,(double)(x+y)); 408 } 409 po++; 410 } 411 cvNamedWindow("frequence_domain",CV_WINDOW_AUTOSIZE); 412 cvShowImage("frequence_domain",ipl); 413 cvNamedWindow("gauss_fre",CV_WINDOW_AUTOSIZE); 414 cvShowImage("gauss_fre",gausf); 415 cvNamedWindow("fourier_t",CV_WINDOW_AUTOSIZE); 416 cvShowImage("fourier_t",img); 417 cvWaitKey(0); 418 cvReleaseImage(&ipl); 419 cvReleaseImage(&gausf); 420 cvReleaseImage(&img); 421 cvDestroyWindow("fourier_t"); 422 cvDestroyWindow("frequence_domain"); 423 cvDestroyWindow("gauss_fre"); 424 delete pixval; 425 delete pixvals; 426 delete src; 427 delete twodfourier; 428 delete twodifourier; 429 return 1; 430 } 431 return 0; 432 }
显示输出：
原图
 A |
 B |
|
C |
D |
A为原图，B为使用均值为0，方差为5.0的高斯函数做低通滤波后的效果图。C为原图的傅立叶频谱图（已居中显示），D为使用方差5.0的高斯滤波后的频谱图；
