ffwt进行fft

来源：互联网发布：商桥2016软件下载编辑：程序博客网时间：2024/06/06 09:13

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

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

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

一、Windows安装FFTW

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

#include "fftw3.h"
#include <stdio.h>
#define N 8
int 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

[cpp] view plain copy

/* 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;
else
i0 = i - height / 2;
if (j < width / 2)
j0 = j + width / 2; // method 2
else
j0 = 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[百度文库]

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

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

#include <complex>#include <fftw3.h>#include <math.h>#include <iostream> #define N 10using 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;}

0 0