離散傅里葉變換的實現與在圖像中的應用

c++實現二維離散傅里葉變換：

基於快速傅里葉蝶形算法

  1 #include <stdio.h>
  2 #include <iostream>
  3 #include <complex>
  4 #include <bitset>
  5 #include <fstream>
  6 #include <algorithm>
  7 #include <iterator>
  8 #include <math.h>
  9 #include "cv.h"
 10 #include "highgui.h"
 11 #include "CImg.h"
 12 #define pi 3.1415926535
 13 using std::iostream;
 14 using std::bitset;
 15 using std::complex;
 16 using namespace std;
 17 using namespace cimg_library;
 18 
 19 int rev(int k,int n)
 20 {
 21      bitset<32> tmp(k);
 22      bitset<32> dist;
 23      for(int m = 0;m<n;m++)
 24      {
 25         if(tmp.test(m))
 26         {
 27             dist.set(n-1-m);
 28         }
 29      }
 30      int revk = (int)dist.to_ulong();
 31      return revk;
 32 }
 33 //重載形式
 34 complex<double>* FFT(const complex<double> *srcimg,int n)
 35 {
 36     double flag = log10((double)n)/log10(2.0);
 37     int N = n;
 38     if(flag - (int)flag != 0){   //判斷是否爲2的指數次
 39         cout <<"the length of srcimg is wrong"<< endl;
 40         /*填充成2的指數項*/
 41         cout <<"need padding"<<endl;
 42         N = pow(2,(double)((int)flag+1));
 43         flag = log10((double)N)/log10(2.0);
 44     }
 45     /*改變成奇偶順序*/
 46     complex<double> *arr= new complex<double>[N];
 47     int sub;
 48     for(int k =0;k<N;k++)
 49     {
 50         sub =rev(k,(int)flag); 
 51         if(sub <= n-1){
 52             arr[k] = *(srcimg + sub);
 53         }else{
 54             complex<double> t = complex<double>(0,0);
 55             arr[k] = t;
 56         }
 57     }
 58 //    cout<<"------------after padding and retrival-----------------"<<endl;
 59 //    for(size_t y=0;y<N;y++)
 60 //    {
 61 //        cout << arr[y].real()<<"  "<<arr[y].imag()<<endl;
 62 //    }
 63 
 64     /*基於迭代的蝶形快速傅立葉變換，自底向上*/
 65      for(int s =1;s <= flag;s++)
 66     {
 67         int m = pow(2,(double)s);
 68         complex<double> wm = complex<double>(cos(2*pi/m),sin(2*pi/m));//wm1：從W21開始,週期變換
 69         for(int p = 0;p<N;p+=m)
 70         {
 71             complex<double> w(1,0);
 72             for(int j = 0;j<=m/2-1;j++)
 73             {
 74                 complex<double> t = w * arr[p+j+m/2];
 75                 complex<double> u = arr[p+j];
 76                 arr[p+j] = u+t;
 77                 arr[p+j+m/2] = u-t;
 78                 w = w*wm;
 79             }
 80         }
 81     }
 82     return arr;
 83 }
 84 /***********一維快速傅立葉變換********************
 85 *srcimg : 原始一維序列                          *
 86 *n     ：一維序列的長度
 87 *************************************************/
 88 complex<double>* FFT(const double *srcimg,int n)
 89 {
 90     double flag = log10((double)n)/log10(2.0);
 91     int N = n;
 92     if(flag - (int)flag != 0){   //判斷是否爲2的指數次
 93 //        cout <<"the length of srcimg is wrong"<< endl;
 94         /*填充成2的指數項*/
 95 //        cout <<"need padding"<<endl;
 96         N = pow(2,(double)((int)flag+1));
 97         flag = log10((double)N)/log10(2.0);
 98     }
 99     /*改變成奇偶順序*/
100     complex<double> *arr= new complex<double>[N];
101     int sub;
102     for(int k =0;k<N;k++)
103     {
104         sub =rev(k,(int)flag); 
105         if(sub <= n-1){
106             arr[k] = complex<double>(*(srcimg + sub),0);
107         }else{
108             complex<double> t = complex<double>(0,0);
109             arr[k] = t;
110         }
111     }
112 //    cout<<"------------after padding and retrival-----------------"<<endl;
113 //    for(size_t y=0;y<N;y++)
114 //    {
115 //        cout << arr[y].real()<<"  "<<arr[y].imag()<<endl;
116 //    }
117 
118     /*基於迭代的蝶形快速傅立葉變換，自底向上*/
119      for(int s =1;s <= flag;s++)
120     {
121         int m = pow(2,(double)s);
122         complex<double> wm = complex<double>(cos(2*pi/m),sin(2*pi/m));//wm1：從W21開始,週期變換
123         for(int p = 0;p<N;p+=m)
124         {
125             complex<double> w(1,0);
126             for(int j = 0;j<=m/2-1;j++)
127             {
128                 complex<double> t = w * arr[p+j+m/2];
129                 complex<double> u = arr[p+j];
130                 arr[p+j] = u+t;
131                 arr[p+j+m/2] = u-t;
132                 w = w*wm;
133             }
134         }
135     }
136     return arr;
137 }
138 int countPadding(int n);
139 /*****************一維快速傅立葉逆變換********************/
140 /*fftimg:原始一維傅立葉序列
141   n     : 序列長度
142 *******************************************************/
143 complex<double>* IFFT(const complex<double> *fftimg,int n)
144 {
145     n = countPadding(n);
146     double flag = log10((double)n)/log10(2.0);
147     int N = n;
148     if(flag - (int)flag != 0){   //判斷是否爲2的指數次
149         cout <<"the length of srcimg is wrong"<< endl;
150         /*填充成2的指數項*/
151         cout <<"need padding"<<endl;
152         N = pow(2,(double)((int)flag+1));
153         flag = log10((double)N)/log10(2.0);
154     }
155     /*改變成奇偶順序*/
156     complex<double> * spit = new complex<double>[N];
157     int sub=0;
158     for(int k =0;k<N;k++)
159     {
160         sub =rev(k,(int)flag); 
161         if(sub < n){
162             spit[k] = complex<double>(*(fftimg + sub));
163         }else{
164             spit[k] = complex<double>(0,0);
165         }
166     }
167 
168     for(int s =1;s <= flag;s++)
169     {
170         int m = pow(2,(double)s);
171         complex<double> wm = complex<double>(cos(-2*pi/m),sin(-2*pi/m));//wm1：從W2(-1)開始
172         for(int p = 0;p<N;p+=m)
173         {
174             complex<double> w(1,0);
175             for(int j = 0;j<=m/2-1;j++)
176             {
177                 complex<double> t = w * spit[p+j+m/2];
178                 complex<double> u = spit[p+j];
179                 spit[p+j] = u+t;
180                 spit[p+j+m/2] = u-t;
181                 w = w*wm;
182             }
183         }
184     }
185 
186     for(size_t p =0;p<n;p++)
187     {
188         spit[p] = spit[p]/complex<double>(N,0);
189     }
190     return spit;
191 }
192 /*******使用共軛的快速傅立葉逆變換**************/
193 complex<double>* IFFT2(const complex<double> *fftimg,int n)
194 {
195     n = countPadding(n);
196     complex<double> *gfftimg = new complex<double>[n];
197     for(size_t i = 0;i<n;i++){
198         gfftimg[i] = complex<double>(fftimg[i].real(),-fftimg[i].imag());
199     }
200     complex<double> *ifft = FFT(gfftimg,n); 
201     for(size_t j = 0;j<n;j++)
202     {
203         ifft[j] = ifft[j]/complex<double>(n,0);
204     }
205     delete gfftimg;
206     return ifft;
207 }
208 
209 /*************二維快速傅立葉變換**************************
210 *srcimg:用一維表示的二維原始序列
211 *width ：寬度
212 height：高度
213 ********************************************************/
214 complex<double>* twoDFFT(const double *srcimg,int width,int height)
215 {
216     int w = countPadding(width);
217     int h = countPadding(height);
218     int pixes = w * h;
219     complex<double> *hdirection = new complex<double>[w];
220     complex<double> *vdirection = new complex<double>[h];
221     complex<double> *fourier = new complex<double>[pixes];
222     /********移位居中************/
223     /*二維水平方向*/
224     for(size_t i = 0;i<h;i++){
225         for(size_t j = 0;j<w;j++){
226             if(i>=height || j >=width){
227                 hdirection[j] = complex<double>(0,0);
228             }else{
229                 hdirection[j] = complex<double>(srcimg[i*width + j],0);
230             }
231     //        cout << hdirection[j] << " ";
232         }
233     //    cout <<""<<endl;
234         complex<double> *hfourier = FFT(hdirection,w);
235         for(size_t m = 0;m<w;m++){
236             fourier[i*w+m] = hfourier[m];
237         }
238         delete hfourier;
239     }
240     /*二維垂直方向*/
241     for(size_t ii = 0;ii<w;ii++){
242         for(size_t jj = 0;jj<h;jj++){
243             vdirection[jj] = fourier[jj*w + ii];
244         }
245         complex<double> *vfourier = FFT(vdirection,h);
246         for(size_t mm = 0;mm < h;mm++){
247             fourier[mm*w +ii] = vfourier[mm];
248         }
249         delete vfourier;
250     }
251     delete hdirection;
252     delete vdirection;
253     return fourier;
254 
255 }
256 /**************二維快速傅立葉逆變換*************************
257 *fourier : 一維表示的二維傅立葉變換序列
258 *width:寬
259 *height:高
260 ********************************************/
261 
262 complex<double>* twoDIFFT(const complex<double> *fourier,int width,int height)
263 {
264     width = countPadding(width);
265     height = countPadding(height);
266     int fpoints = width * height;
267     complex<double> *hdirection = new complex<double>[width];
268     complex<double> *vdirection = new complex<double>[height];
269     complex<double> *ifourier = new complex<double>[fpoints];
270 
271     for(size_t ii = 0;ii<height;ii++)
272     {
273         for(size_t jj = 0;jj<width;jj++){
274             hdirection[jj] = fourier[ii*width+jj];
275         }
276         complex<double> *hifour = IFFT(hdirection,width);//臨時變量
277         for(size_t mm = 0;mm<width;mm++){
278             ifourier[ii*width+mm] = hifour[mm];
279         }
280         delete hifour;
281     }
282     for(size_t i = 0;i<width;i++){
283         for(size_t j = 0;j<height;j++){
284             vdirection[j] = ifourier[j*width+i];
285         }
286         complex<double> *vifour = IFFT(vdirection,height);
287         for(size_t m = 0;m<height;m++){
288             ifourier[m*width+i] = vifour[m];
289         }
290         delete vifour;
291     }    
292     delete hdirection;
293     delete vdirection;
294     return ifourier;
295 }
296 /******計算填充數*********************/
297 int countPadding(int n)
298 {
299     double lg = log10((double)n)/log10(2.0);
300     if((lg - (int)lg) == 0){
301         return n;
302     }
303     int N = pow(2.0,((int)lg+1));
304     return N;
305 }
306 /*測試一維填充形式傅立葉變換*/
307 /*
308 int main()
309 {
310     ifstream infile("D:\\ftest.txt");
311     if(!infile.is_open())
312     {
313         cout<<"file read error"<<endl;
314         return -1;
315     }
316     istream_iterator<double> df(infile);
317     istream_iterator<double> eof;
318     int N = *df++;
319     double *src = new double[N];
320     size_t i = 0;
321     while(df != eof)
322     {
323         src[i] = *df++;
324         i++;
325     }
326 /*    for(size_t a = 0;a < N;a++)
327     {
328         cout << src[a]<<endl;
329     }
330 */    
331 //    complex<double> *fourier = new complex<double>[N];
332 //    complex<double> *spit = new complex<double>[N];
333 /*    complex<double> *fourier = FFT(src,N);
334     cout <<"--------fourier sieres----------"<<endl;
335     for(size_t x = 0;x <16;x++)
336     {
337         cout << (*(fourier +x)).real()<<"  "<<(*(fourier+x)).imag() << endl;
338     }
339     cout <<"----------------inverse fourier -----------------------"<<endl;
340     complex<double> *spit = IFFT(fourier,N);
341     for(size_t y = 0;y <16;y++)
342     {
343         cout << spit[y].real()<<"  "<<spit[y].imag() << endl;
344     }
345     delete src;
346 //    delete fourier;
347 //    delete spit;
348     return 1;
349 }*/
350 /*測試二維零填充傅立葉正負變換*/
351 
352 /*
353 int main()
354 {
355     ifstream infile("D:\\ftest2.txt");
356     if(!infile.is_open())
357     {
358         cout<<"file read error"<<endl;
359         return -1;
360     }
361     istream_iterator<double> df(infile);
362     istream_iterator<double> eof;
363     int N = *df++;
364     double *src = new double[N];
365     size_t i = 0;
366     while(df != eof)
367     {
368         src[i] = *df++;
369         i++;
370     }
371     complex<double> *fourier = twoDFFT(src,3,6);
372     cout <<"--------fourier sieres----------"<<endl;
373     for(size_t x = 0;x <32;x++)
374     {
375         cout << (*(fourier +x)).real()<<"  "<<(*(fourier+x)).imag() << endl;
376     }
377     cout <<"----------------inverse fourier -----------------------"<<endl;
378     complex<double> *spit = twoDIFFT(fourier,3,6);
379     for(size_t y = 0;y <32;y++)
380     {
381         cout << spit[y].real()<<"  "<<spit[y].imag() << endl;
382     }
383     delete src;
384 //    delete fourier;
385 //    delete spit;
386     return 1;
387 }*/
388 
389 /*測試一維傅里葉圖像變換*//*
390 int main(int argc,char ** argv)
391 {
392     IplImage *img;
393     img = cvLoadImage("F:\\Fig0222(a)(face).tif",0);
394     int dim = img->imageSize;
395     //從圖像矩陣複製出單維數組；
396     double * src = new double[dim];
397     size_t j =0;
398     for(int y =0;y<img->height;y++)
399     {
400         uchar * ptr = (uchar*)(img->imageData + y * img->widthStep);
401         for(int x =0;x<img->width;x++){
402             src[j] = (double)ptr[x]/256;
403             j++;
404         }
405     }
406     //一維數組作傅立葉變換
407     complex<double>* fourier = FFT(src,dim);
408     //計算填充後傅氏數組大小
409     double lg = log10((double)dim)/log10(2.0);
410     int n = pow(2,(double)((int)lg+1));
411     //傅立葉逆變換
412     complex<double> *ifourier = IFFT(fourier,n);
413     double ipix = 0;
414     size_t po = 0;
415     //重填充圖像
416     for(int y =0;y<img->height;y++)
417     {
418         uchar * ptr = (uchar*)(img->imageData + y * img->widthStep);
419         for(int x =0;x<img->width;x++){
420             ipix = ifourier[po].real();
421             ptr[x] = (uchar)(ipix * 256);
422             po++;
423         }
424     }
425     cvNamedWindow("fourier_t",CV_WINDOW_AUTOSIZE);
426     cvShowImage("fourier_t",img);
427     cvWaitKey(0);
428     cvReleaseImage(&img);
429     cvDestroyWindow("fourier_t");
430 
431     return 1;
432 /*    ifstream infile("F:\\Fig0222(a)(face).tif",ios::binary);
433     vector<double> collection;
434     istream_iterator<byte> in_iter(infile);
435     istream_iterator<byte> eof;
436     double temp = 0;
437     unsigned char c = 0;
438     while(in_iter != eof)
439     {
440         c = *in_iter++;
441         temp = (double)c/256;
442         collection.push_back(temp);
443     }
444     int amount = collection.size();
445     double * src = new double[amount];  
446     size_t j = 0;
447     for(vector<double>::iterator iter = collection.begin();iter != collection.end();++iter)
448     {
449         src[j] = *iter;
450         j++;
451     }*/
452 /*    delete src;
453     delete fourier;
454     delete ifourier;
455 } 
456 */
457 
458 int main(int argc,char **argv)
459 {
460     IplImage *img;
461     if((img = cvLoadImage("F:\\Fig0222(a)(face).tif",0))!=0){
462         int dim = img->imageSize;
463     //從圖像矩陣複製出單維數組；
464         double * src = new double[dim];
465         size_t j =0;
466         for(int y =0;y<img->height;y++)
467         {
468             uchar * ptr = (uchar*)(img->imageData + y * img->widthStep);
469             for(int x =0;x<img->width;x++){
470                 src[j] = (double)ptr[x]*pow(-1,(double)(x+y))/256;
471                 j++;
472             }
473         }
474         int w = img->width;
475         int h = img->height;
476         int pwid = countPadding(w);
477         int phei = countPadding(h);
478         complex<double> *twodfourier = twoDFFT(src,w,h);
479         double *vals = new double[pwid*phei];
480         char * pixval = new char[pwid*phei];
481         CvMat freq;
482         double max = 0;
483         double min = 255;
484         for(size_t p = 0;p<pwid*phei;p++){
485             vals[p]=log(1+sqrt(pow(twodfourier[p].real(),2.0)+pow(twodfourier[p].imag(),2.0)));//對數級的幅度鋪
486             if(vals[p] > max){
487                 max = vals[p];
488             }
489             if(vals[p] < min){
490                 min = vals[p];
491             }
492         }
493         cout<<min<< " "<<max<<endl;
494         for(size_t s = 0;s<pwid*phei;s++){
495             pixval[s] = (char)((vals[s]-min)/(max-min)*255);
496         }
497 /*        for(size_t q = 0;q < pwid*phei;q++){
498             cout <<vals[q]<<" "<<endl;
499         }*/
500         cvInitMatHeader(&freq,pwid,phei,CV_8UC1,pixval);
501         IplImage *ipl = cvCreateImage(cvGetSize(&freq),8,1);
502         cvGetImage(&freq,ipl);
503         complex<double> *twodifourier = twoDIFFT(twodfourier,w,h);
504         double ipix = 0;
505         size_t po = 0;
506         for(int y =0;y<img->height;y++)
507         {
508             uchar * ptr = (uchar*)(img->imageData + y * img->widthStep);
509             for(int x =0;x<img->width;x++){
510                 ipix = twodifourier[po*pwid+x].real();
511                 ptr[x] = (uchar)(ipix * 256)*pow(-1,(double)(x+y));
512             }
513             po++;
514         }
515         cvNamedWindow("frequence domain",CV_WINDOW_AUTOSIZE);
516         cvShowImage("frequence domain",ipl);
517         cvNamedWindow("fourier_t",CV_WINDOW_AUTOSIZE);
518         cvShowImage("fourier_t",img);
519         cvWaitKey(0);
520         cvReleaseImage(&ipl);
521         cvReleaseImage(&img);
522         cvDestroyWindow("fourier_t");
523         cvDestroyWindow("frequence domain");
524         delete vals;
525         delete pixval;
526         delete src;
527         delete twodfourier;
528         delete twodifourier;
529         return 1;
530     }
531     return 0;
532 }
533 
534 /*
535 int main()
536 {
537     CImg<unsigned char> srcimg("F:\\Fig0222(a)(face).tif");
538     //CImgDisplay main_disp(srcimg,"lina");
539     srcimg.display();
540     return 1;
541 
542 }
543 */

 輸出：

原圖：

頻譜圖：