FFTW3

在做ITK的配准时发现了document的一个bug，itkcurvatureRegistrationFilter无法显示,在Nabble里看了一下，原来06年就有人提出这个问题了，不过在10年11-24这个问题才解决。解决的方法是：在加上预编译#ifndef USE_FFTWF #define USE_FFTWF
#ifndef USE_FFTWD  #define USE_FFTWD。但是如果你的电脑里没有fftw3（世界上最快的傅立叶变化）那么这个filter还是无法进行，下面的是转的一个兄弟的贴子，他详细的介绍了方法，以及测试结果。

 

Fastest Fourier Transform in the West

         FFTW is a comprehensive collection of fast C routines for computing the discrete Fourier transform (DFT) and various special cases thereof.

        据说是世界上最快的FFT。   

        机子里的FFTW库下了很长时间了，总也没有去搞。唉，有很多东西就是这样，千方百计搞过来，搞到手了就晾在那里了。记得《黄生借书说》里面说得对啊“书非借不能读也”。

         好了，感慨完毕。归入正题。

Windows下FFTW库的安装

1、 从网址http://www.fftw.org/install/windows.html上获得FFTW的windows dll预编译版本；

2、 解压缩文件，打开windows命令行窗口，就是那个cmd窗口啦。然后把当前目录转换到你解压缩文件的目录下。

3、 执行以下3个指令

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

这会在该目录下建三个相应的dll文件和lib文件。注意第三个.def文件中的“3l-3”中的是字母L的小写，不是数字一。因为这个问题，我搞了半个小时，呵呵。

4、 将libfftw3l-3.dll, libfftw3f-3.dll, libfftw3-3.dll 文件复制到文件夹system32中。这一步是为了你以后都不用在你的可执行文件所在的文件夹中带上这3个拖油瓶，因为系统直接会去system32中找。

5、 在 VC 中指定 libfftw3l-3.lib, libfftw3f-3.lib, libfftw3-3.lib 这3个lib文件及 fftw3.h 文件所在的目录。也就是在vc++的tools->options的 Directories选项中的Include Files和Library Files中把这两个目录加上，使得以后VC编译的时候知道该到哪个目录中去找。

6、 最后一步就是，在你新建工程的时候，记得#include“fftw3.h”,然后把你要用的lib写到Project->setting->link->General里面的Object/library modules里面去。

VS2008中则, alt + F7，-> input -> additional dependences 中加入libfftw3l-3.lib   libfftw3f-3.lib libfftw3-3.lib 。注意：三者之间是空格。

7、 下面，你就可以放心大胆地去使用fftw的库编程了。为了熟悉FFTW的调用方式和数据结构，你还可以从http://www.fftw.org/#documentation下载一份manual好好钻研钻研。

可能出现的错误：

1.      LNK1181：cannot open input file“…”: 出现这个错误的原因有: (1)你文件名打错了; (2) 你的当前目录不对，当前目录应该是你解压后的文件目录。

2.      fatal error LNK1104: cannot open file 'libfftw3l-3.obj'或者源文件编译时报找不到lib，那是因为你第5步或第6步没做好，建议重做一下第5步和第6步。

 

使用FFTW编写测试程序

上面的搞好后，就写一个小的测试代码试一下效果。我就抄了一个网上的代码：

//Code
#include "fftw3.h"
int main()
{
    fftw_complex *in, *out;
    fftw_plan p;
    int N= 8;
    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    for( int i=0; i < N; i++)
           {
        in[i][0] = 1.0;
        in[i][1] = 0.0;
        printf("%6.2f ",in[i][0]);
        }
    printf("/n");
    p=fftw_plan_dft_1d(N,in,out, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(p); /* repeat as needed*/
    for(int j = 0;j < N;j++)
          {
        printf("%6.2f ",out[j][0]);
          }
    printf("/n");
    fftw_destroy_plan(p);
    fftw_free(in);
    fftw_free(out);
    return 0;
}

程序给了一个直流的时域数据，应该出来一个只有直流分量的DFT数据。
请看大屏幕：

数据正确，说明我们的调用流程没有问题，就可以使用了。

参考文献：

1、    http://www.fftw.org/install/windows.html

2、    FFTW mannual

3、    http://blog.chinson.idv.tw/ 好像不太能登上，请用网页快照，呵呵。

 

例子二：

 

今天测试了FFTW-Windows下的DLL。本部分为第一部分，介绍了基本的用法，如何动态分配数组，释放数组，执行FFT等。在第二部分中，将介绍如何将PLAN写入文件。

要用lib命令，生成.lib，如lib /def:libfftw3-3.def。代码如下：

// Testing code for FFTW using VS7.
//
// Written by Weijia Sun, November 24, 2008.
// Copyright reserved.

#include "stdafx.h"
#include <complex>
#include "fftw3.h"

int _tmain(int argc, _TCHAR* argv[])
{
int N=16000; // length of FFT
fftw_complex *in,*out, *in2; // array for I/O

cout << "Size of fftw_complex: " << sizeof(fftw_complex) << endl;
in = (fftw_complex*) fftw_malloc(N * sizeof(fftw_complex)); // for Inputting
in2 = (fftw_complex*) fftw_malloc(N * sizeof(fftw_complex)); // for output
out = (fftw_complex*) fftw_malloc(N * sizeof(fftw_complex)); // for comparison

// Initializing ...
for(int i = 0; i < N; i++){
   in[i][0] = i;
   in[i][1] = 0;
}

// making a plan and execute it, then destroy it. Forward FFT
fftw_plan pf = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(pf);
fftw_destroy_plan(pf);

// making a plan and execute it, then destroy it. Backward FFT
fftw_plan pb = fftw_plan_dft_1d(N, out, in2, FFTW_BACKWARD, FFTW_ESTIMATE);
fftw_execute(pb);
fftw_destroy_plan(pb);

// output the fore 20 results
for(i = 0; i < 20; i++){
   cout << in[i][0] << " " << in2[i][0]/N << endl;
}

// free space
fftw_free(in);
fftw_free(in2);
fftw_free(out);

return 0;
}

例子三：

使用FFTW3做二维DFT的示例代码

#include <fftw3.h>
#include <stdio.h>

#define N 3
#define ELEM(r,c) (r*N+c)

int showresult(fftw_complex* in, fftw_complex* out)
{
    int i, j;
     printf("In:/n");
    for (i=0; i<N; i++) {
        for (j=0; j<N; j++) {
             printf("%lf/t", in[ELEM(i, j)][0]);
         }
         printf("/n");
     }
     printf("Out:/n");
    for (i=0; i<N; i++) {
        for (j=0; j<N; j++) {
             printf("%lf/t", out[ELEM(i, j)][0]);
         }
         printf("/n");
     }
    return 1;
}

int main()
{
     fftw_complex *in, *out;
     fftw_plan p;
    int i, j;

    // 分配存储空间
    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N * N);
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N * N);

    // 设置变换计划
     p = fftw_plan_dft_2d(N, N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);

    // 设置测试数据
    for (i=0; i<N; i++) {
        for (j=0; j<N; j++) {
            in[ELEM(i, j)][0] = 1;
            in[ELEM(i, j)][1] = 0;
         }
     }
    in[ELEM(1, 2)][0] = 5;
    in[ELEM(2, 2)][0] = 3;

    // 执行变换
     fftw_execute(p); /**//* repeat as needed */

    // 显示测试结果
     showresult(in, out);

    // 释放内存
     fftw_destroy_plan(p);
     fftw_free(in);
     fftw_free(out);

    return 1;
}

makefile：

all: fft_test

fft_test: fft_test.c
     gcc -o fft_test fft_test.c -lfftw3

clean:
     rm -f fft_test

结果：

In:
1.000000    1.000000    1.000000    
1.000000    1.000000    5.000000    
1.000000    1.000000    3.000000    
Out:
15.000000     -3.000000     -3.000000    
-3.000000    3.000000    0.000000    
-3.000000    0.000000    3.000000