fftw进行fft

版权声明：本文为博主原创文章，未经博主允许不得转载。

FFTW首页： http://www.fftw.org/

据说FFTW是世界上最快的FFT。

下载地址： http://www.fftw.org/download.html

一、Windows安装FFTW

1. 从网址http://www.fftw.org/install/windows.html 上获得FFTW的windows dll预编译版本；
2. 解压缩文件，打开windows命令行窗口，就是那个cmd窗口。然后把当前目录转换到你解压缩文件的目录下。
3. 执行以下3个指令

        这里的def可以加路径  def:D:\fftw\libfftw3-3.def 这样。

               lib/machine:ix86/def:libfftw3-3.def

               lib/machine:ix86/def:libfftw3f-3.def

               lib/machine:ix86/def:libfftw3l-3.def

        这会在该目录下建三个相应的dll文件和lib文件

3-3对应double

3f-3对应float

3l-3d

fftw 3.3.4 32bit dll，lib下载： http://pan.baidu.com/s/1sjIAVKT  提取码：xmta

二、VS配置FFTW

• 链接器  |  输入  |  附加依赖项    加入libfftw3-3.lib，libfftw3f-3.lib，libfftw3l-3.lib 
• libfftw3-3.dll，libfftw3f-3.dll，libfftw3l-3.dll放在当前工程目录（或system32）或在vs属性中的Debugging  |  Environment 中添加路径（例：path=%path%;.\fftw\bin）
• 记得include“fftw3.h”

     测试demo

[cpp] view plain copy

1. #include "fftw3.h"  
2. #include <stdio.h>  
3. #define N 8  
4. int main()  
5. {  
6.     int i;  
7.     fftw_complex *din,*out;  
8.     fftw_plan p;  
9.     din  = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);  
10.     out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);  
11.     if((din==NULL)||(out==NULL))  
12.     {  
13.         printf("Error:insufficient available memory\n");  
14.     }  
15.     else  
16.     {  
17.         for(i=0; i<N; i++)/*测试数据*/  
18.         {  
19.             din[i][0] = i+1;  
20.             din[i][1] = 0;  
21.         }  
22.     }  
23.     p = fftw_plan_dft_1d(N, din, out, FFTW_FORWARD,FFTW_ESTIMATE);  
24.     fftw_execute(p); /* repeat as needed */  
25.     fftw_destroy_plan(p);  
26.     fftw_cleanup();  
27.     for(i=0;i<N;i++)/*OUTPUT*/  
28.     {  
29.         printf("%f,%fi\n",din[i][0],din[i][1]);  
30.     }  
31.     printf("\n");  
32.     for(i=0;i<N;i++)/*OUTPUT*/  
33.     {  
34.         printf("%f,%fi\n",out[i][0],out[i][1]);  
35.     }  
36.   
37.     if(din!=NULL) fftw_free(din);  
38.     if(out!=NULL) fftw_free(out);  
39.     getchar();  
40.     return 0;  
41. }  

结果输出如下：

三、对图像做FFT的demo

[cpp] view plain copy

1. /* load original image */  
2. IplImage *img_src = cvLoadImage("a.bmp", CV_LOAD_IMAGE_GRAYSCALE);  
3. if (img_src == 0)  
4. {  
5.     std::cout << "cannot load file" << std::endl;  
6.     return 0;  
7. }  
8.   
9. /* create new image for FFT & IFFT result */  
10. IplImage *img_fft = cvCreateImage(cvSize(img_src->width, img_src->height), IPL_DEPTH_8U, 1);  
11. IplImage *img_ifft = cvCreateImage(cvSize(img_src->width, img_src->height), IPL_DEPTH_8U, 1);  
12.   
13. /* get image properties */  
14. int width = img_src->width;  
15. int height = img_src->height;  
16. int step = img_src->widthStep;  
17. uchar *img_src_data = (uchar *)img_src->imageData;  
18. uchar *img_fft_data = (uchar *)img_fft->imageData;  
19. uchar *img_ifft_data = (uchar *)img_ifft->imageData;  
20.   
21. /* initialize arrays for fftw operations */  
22. fftw_complex *data_in = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * width * height);  
23. fftw_complex *fft = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * width * height);  
24. fftw_complex *ifft = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * width * height);  
25.   
26. /* create plans */  
27. fftw_plan plan_f = fftw_plan_dft_2d(height, width, data_in, fft, FFTW_FORWARD, FFTW_ESTIMATE);  
28. fftw_plan plan_b = fftw_plan_dft_2d(height, width, fft, ifft, FFTW_BACKWARD, FFTW_ESTIMATE);  
29.   
30. int i, j, k;  
31. /* load img_src's data to fftw input */  
32. for (i = 0, k = 0; i < height; ++i)  
33. {  
34.     for (j = 0; j < width; ++j)  
35.     {  
36.         // method 1: 输入数据乘以(-1)^（i+j），即可中心化  
37.         data_in[k][0] = /*pow(-1, i + j) * */(double)img_src_data[i * step + j];    
38.         data_in[k][1] = 0.0;  
39.         k++;  
40.     }  
41. }  
42.   
43. /* perform FFT */  
44. fftw_execute(plan_f);  
45.   
46. /* perform IFFT */  
47. fftw_execute(plan_b);  
48.   
49. /* normalize FFT result */  
50. double maxx = 0.0, minn = 10000000000.0;  
51. for (i = 0; i < width * height; ++i)  
52. {  
53.     fft[i][0] = log(sqrt(fft[i][0] * fft[i][0] + fft[i][1] * fft[i][1]));  
54.     maxx = fft[i][0] > maxx ? fft[i][0] : maxx;  
55.     minn = fft[i][0] < minn ? fft[i][0] : minn;  
56. }  
57.   
58. for (i = 0; i < width * height; ++i)  
59. {  
60.     fft[i][0] = 255.0 * (fft[i][0] - minn) / (maxx - minn);  
61. }  
62.   
63. /* copy FFT result to img_fft's data */  
64. int i0, j0;  
65. for (i = 0, k = 0; i < height; ++i)  
66. {  
67.     for (j = 0; j < width; ++j)  
68.     {  
69.         if (i < height / 2)  
70.             i0 = i + height / 2;  
71.         else  
72.             i0 = i - height / 2;  
73.         if (j < width / 2)  
74.             j0 = j + width / 2;   // method 2  
75.         else  
76.             j0 = j - width / 2;  
77.   
78.         img_fft_data[i * step + j] = (uchar)fft[/*k++*/i0 * width + j0][0];  
79.     }  
80. }  
81.   
82. /* normalize IFFT result */  
83. for (i = 0; i < width * height; ++i)  
84. {  
85.     ifft[i][0] /= width * height;  
86. }  
87.   
88. /* copy IFFT result to img_ifft's data */  
89. for (i = 0, k = 0; i < height; ++i)  
90. {  
91.     for (j = 0; j < width; ++j)  
92.     {  
93.         img_ifft_data[i * step + j] = (uchar)ifft[k++][0];  
94.     }  
95. }  
96.   
97. /* display images */  
98. cvNamedWindow("original_image", CV_WINDOW_AUTOSIZE);  
99. cvNamedWindow("FFT", CV_WINDOW_AUTOSIZE);  
100. cvNamedWindow("IFFT", CV_WINDOW_AUTOSIZE);  
101. cvShowImage("original_image", img_src);  
102. cvShowImage("FFT", img_fft);  
103. cvShowImage("IFFT", img_ifft);  
104.   
105. cvWaitKey(0);  
106.   
107. /* free memory */  
108. cvDestroyWindow("original_image");  
109. cvDestroyWindow("FFT");  
110. cvDestroyWindow("IFFT");  
111. cvReleaseImage(&img_src);  
112. cvReleaseImage(&img_fft);  
113. cvReleaseImage(&img_ifft);  
114. fftw_destroy_plan(plan_f);  
115. fftw_destroy_plan(plan_b);  
116. fftw_free(data_in);  
117. fftw_free(fft);  
118. fftw_free(ifft);  

参考资料：

[1]  FFTW_Documentation_CN[百度文库]

还有一段例子：
例如，如果需要用到双精度的实数FFT变换/FFTs，那么在编译的链接命令中需要按如下顺序加入
-ldrfftw -ldfftw参数

下面的是如何使用的一个例子

#include <complex>
#include <fftw3.h>
#include <math.h>
#include <iostream>
 
#define N 10
using namespace std;
int main(int argc, char * argv[]){
 
    fftw_complex in[N], out[N];
    fftw_plan p;
    p=fftw_plan_dft_1d(N,in,out,FFTW_FORWARD,FFTW_MEASURE);
    for(int i=0;i <N;i ++) {
        in[i][0]=i;
        in[i][1]=0.0;
    }
   
 
    fftw_execute(p);
 
    for(int i=0;i <N;i ++){
        cout<<out[i][0]<<" "<<out[i][1]<<endl;
    }
    complex<double> temp = 0.0;
    for(int k =0; k < N; k ++){
        double pi = 4*atan(1.0);
        temp += exp(complex<double>(0.0,-2.0*pi*3*k/N))*complex<double>(in[k][0],in[k][1]);
    }
    cout<<"out[3] is "<<temp<<endl;
 
    fftw_complex out1[N];
 
    fftw_plan p1;
 
    p1=fftw_plan_dft_1d(N,out1,in,FFTW_BACKWARD,FFTW_MEASURE);
 
    for(int i=0;i <N;i ++){
        out1[i][0]=out[i][0];
        out1[i][1]=out[i][1];
    }
 
    fftw_execute(p1);
    for(int i=0;i <N;i ++){
        cout<<in[i][0]<<" "<<in[i][1]<<endl;
    }
 
 
 
    fftw_destroy_plan(p);
    fftw_destroy_plan(p1);
    return 1;
}