FFTW的使用

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

#include "fftw3.h"#include <stdio.h>#define N 8int main(){int i;fftw_complex *din,*out;fftw_plan p;din  = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);if((din==NULL)||(out==NULL)){printf("Error:insufficient available memory\n");}else{for(i=0; i<N; i++)/*测试数据*/{din[i][0] = i+1;din[i][1] = 0;}}p = fftw_plan_dft_1d(N, din, out, FFTW_FORWARD,FFTW_ESTIMATE);fftw_execute(p); /* repeat as needed */fftw_destroy_plan(p);fftw_cleanup();for(i=0;i<N;i++)/*OUTPUT*/{printf("%f,%fi\n",din[i][0],din[i][1]);}printf("\n");for(i=0;i<N;i++)/*OUTPUT*/{printf("%f,%fi\n",out[i][0],out[i][1]);}if(din!=NULL) fftw_free(din);if(out!=NULL) fftw_free(out);getchar();return 0;}

结果输出如下：

三、对图像做FFT的demo

/* load original image */IplImage *img_src = cvLoadImage("a.bmp", CV_LOAD_IMAGE_GRAYSCALE);if (img_src == 0){std::cout << "cannot load file" << std::endl;return 0;}/* create new image for FFT & IFFT result */IplImage *img_fft = cvCreateImage(cvSize(img_src->width, img_src->height), IPL_DEPTH_8U, 1);IplImage *img_ifft = cvCreateImage(cvSize(img_src->width, img_src->height), IPL_DEPTH_8U, 1);/* get image properties */int width = img_src->width;int height = img_src->height;int step = img_src->widthStep;uchar *img_src_data = (uchar *)img_src->imageData;uchar *img_fft_data = (uchar *)img_fft->imageData;uchar *img_ifft_data = (uchar *)img_ifft->imageData;/* initialize arrays for fftw operations */fftw_complex *data_in = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * width * height);fftw_complex *fft = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * width * height);fftw_complex *ifft = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * width * height);/* create plans */fftw_plan plan_f = fftw_plan_dft_2d(height, width, data_in, fft, FFTW_FORWARD, FFTW_ESTIMATE);fftw_plan plan_b = fftw_plan_dft_2d(height, width, fft, ifft, FFTW_BACKWARD, FFTW_ESTIMATE);int i, j, k;/* load img_src's data to fftw input */for (i = 0, k = 0; i < height; ++i){for (j = 0; j < width; ++j){// method 1: 输入数据乘以(-1)^（i+j），即可中心化data_in[k][0] = /*pow(-1, i + j) * */(double)img_src_data[i * step + j];  data_in[k][1] = 0.0;k++;}}/* perform FFT */fftw_execute(plan_f);/* perform IFFT */fftw_execute(plan_b);/* normalize FFT result */double maxx = 0.0, minn = 10000000000.0;for (i = 0; i < width * height; ++i){fft[i][0] = log(sqrt(fft[i][0] * fft[i][0] + fft[i][1] * fft[i][1]));maxx = fft[i][0] > maxx ? fft[i][0] : maxx;minn = fft[i][0] < minn ? fft[i][0] : minn;}for (i = 0; i < width * height; ++i){fft[i][0] = 255.0 * (fft[i][0] - minn) / (maxx - minn);}/* copy FFT result to img_fft's data */int i0, j0;for (i = 0, k = 0; i < height; ++i){for (j = 0; j < width; ++j){if (i < height / 2)i0 = i + height / 2;elsei0 = i - height / 2;if (j < width / 2)j0 = j + width / 2;   // method 2elsej0 = j - width / 2;img_fft_data[i * step + j] = (uchar)fft[/*k++*/i0 * width + j0][0];}}/* normalize IFFT result */for (i = 0; i < width * height; ++i){ifft[i][0] /= width * height;}/* copy IFFT result to img_ifft's data */for (i = 0, k = 0; i < height; ++i){for (j = 0; j < width; ++j){img_ifft_data[i * step + j] = (uchar)ifft[k++][0];}}/* display images */cvNamedWindow("original_image", CV_WINDOW_AUTOSIZE);cvNamedWindow("FFT", CV_WINDOW_AUTOSIZE);cvNamedWindow("IFFT", CV_WINDOW_AUTOSIZE);cvShowImage("original_image", img_src);cvShowImage("FFT", img_fft);cvShowImage("IFFT", img_ifft);cvWaitKey(0);/* free memory */cvDestroyWindow("original_image");cvDestroyWindow("FFT");cvDestroyWindow("IFFT");cvReleaseImage(&img_src);cvReleaseImage(&img_fft);cvReleaseImage(&img_ifft);fftw_destroy_plan(plan_f);fftw_destroy_plan(plan_b);fftw_free(data_in);fftw_free(fft);fftw_free(ifft);

参考资料：

[1] FFTW_Documentation_CN[百度文库]