1. 引言

傅里叶变换是数字信号处理和数字图像处理中的常用算法。FFTW是一个非常优秀的面向C/C++语言的开源一维/多维离散傅里叶变换工具包，由MIT 的 Matteo Frigo and Steven G. Johnson 开发，关于其下载地址详见 http://www.fftw.org/

目前的最新版本是3.3.8；Windows用户可以直接下载3.3.5的官方编译版本。

但是严格来说，FFTW的使用是面向C设计的，浏览了一下国内很多博客页对其用法介绍都是基于单次使用的，但在实际的使用中我们可能需要处理很多同样参数的FFT计算，譬如，对麦克风接收的声音信号实时计算频谱信号。本文开始从一个不怎么好的函数调用，进一步分析FFTW的工作细节，最后给出一个C++中可复用的例子抛砖引玉，如果读者能有更好的实现方式，欢迎留言讨论。

2.从一个不太友好的调用函数开头

首先声明，这个函数是绝对没有问题且可用的。以下是一个朋友刚看了一些博客教程以后写的一个计算信号的幅频响应的（PSD）的函数，如果对傅里叶变换和FFTW一无所知，在看了一些教程后，是非常有可能写出这样的函数的。

以3.3版本为例，首先需要把fft编译出来的dll文件和lib文件添加到工作目录

#include "fftw3.h"

#pragma comment(lib, "libfftw3-3.lib") // double版本

备注：这里我们用的是double版本；FFTW包中对应的有三个版本，如图所示

根据使用需要，在 comment 语句的时候，libfftw3-3是double 版本，libfftw3f-3是float版本，libfftw3l-3是long double版本；

调用函数如下：

void function MyPSD(double in[], double out[], int Nfft)
{
    // Part I
    double *fftIn =  (double*)fftw_malloc(sizeof(double)* Nfft);		    //fft 输入缓存
	fftw_complex *fftOut =(double*)fftw_malloc(sizeof(double)* Nfft);	    //fft 输出缓存

    fftw_plan fftPlan = fftw_plan_dft_r2c_1d(Nfft, fftIn, fftOut, FFTW_ESTIMATE); //一维正向傅里叶变换
    
    //Part II;
    for(auto i = 0;i<Nfft;i++)
        fftIn[i] = in[i];
    
    fftw_execute(fftPlan);

    for(auto i = 0;i<(Nfft/2+1);i++)
        Out[i] = (fftOut[i][0]*fftOut[i][0]+fftOut[i][1]*fftOut[i][1])/Nfft;
    
    // Part III
    fftw_free(fftIn);
    fftw_free(fftOut); 
    fftw_destroy_plan(fftPlan);        // 销毁fftw初始化的内容

    return;
}

函数在功能上一点问题都没有，但是，在实际使用中速度是比较慢的，原因将在下一节中详细展开。

3.逐行分析FFTW的工作流程结构

我在上一段函数的注释中将用注释分了三段，从功能来说；

• Part I 准本工作
• Part II 传入数据，执行计算，拷贝数据
• Part III 内存释放

其实，FFTW的工作原理和STL的设计有些类似，将工作分解为：输入数据容器，输出数据容器，算法执行器构成。为了快速计算需要，FFTW针对确定点数的FFT运算，实时上会预生成一些参量和申请特定的空间。

double *fftIn =  (double*)fftw_malloc(sizeof(double)* Nfft);		    //fft 输入缓存
fftw_complex *fftOut =(double*)fftw_malloc(sizeof(double)* Nfft);	    //fft 输出缓存

这里的fftin 和 fftw_complex 定义的 fftOut就是输入和输出缓存；double 类型用 molloc 初始化和上述语句是没有区别的；

初始化FFT计算的参数，关联内存

fftw_plan fftPlan = fftw_plan_dft_r2c_1d(Nfft, fftIn, fftOut, FFTW_ESTIMATE); 

这一步很重要，我们通过一个叫fftw_plan 类型的东西，将输入容器，输出容器，FFT的初始化参数

• fftw_plan_dft_r2c_1d: 一维正向傅里叶变换
  • Nfft:    FFT 点数
  • fftIn:    输入缓存地址
  • fftOut: 输出缓存地址
  • FFTW_ESTIMATE： FFT的计算精度，宏定义；用FFTW_MEASURE 精度会更高一些，相应速度更慢

拷贝数据

for(auto i = 0;i<Nfft;i++)
        fftIn[i] = in[i];
    

执行FFT计算

    fftw_execute(fftPlan);

  

计算并输出功率谱

  for(auto i = 0;i<(Nfft/2+1);i++)
        Out[i] = (fftOut[i][0]*fftOut[i][0]+fftOut[i][1]*fftOut[i][1])/Nfft;

销毁内存

    fftw_free(fftIn);
    fftw_free(fftOut); 

事实上，在使用过程中，准备（初始化）与收尾（销毁内存）对于 重复执行很多相同参数的傅里叶变换来说，是多余的。初始化过程实际上占用了这个函数的绝大部分时间。

4. 一个关于C++中复用的FFT计算器类的实现

基于3 中的分析，我们很简单的可以知道，对于重复使用的过程在每次计算中保留，其他的东西只做一次最好。

以下是一个简单的改进（备注：并没有考虑各类安全问题，仅做功能说明）

我们将初始化放到构造函数，内存清理放到析构函数，PSD函数就仅包含注入缓存，计算和提取输出

class MyFFT
{
private:
    int Nfft;
    double *_fftIn;
    fftw_complex *_fftOut;
    fftw_plan _fftPlan;
public:
    MyFFT()
    {
        Nfft   = 0;
        _fftIn  = nullptr;
        _fftOut = nullptr;
    }
    MyFFT(int N)
    {
        Nfft = N;
        _fftIn =  (double*)fftw_malloc(sizeof(double)* Nfft);
        _fftOut =(double*)fftw_malloc(sizeof(double)* Nfft);
        _fftPlan = fftw_plan_dft_r2c_1d(Nfft, fftIn, fftOut, FFTW_ESTIMATE);
    }
    ~MyFFT
    {
        if(fftIn != nullptr)
            fftw_free(_fftIn);
        if(fftOut != nullptr)
            fftw_free(_fftOut);
        fftw_destroy_plan(_fftPlan);
    }
    void PSD(double in[], double out[])
    {
        for(auto i = 0;i<Nfft;i++)
            _fftIn[i] = in[i];
    
        fftw_execute(_fftPlan);

        for(auto i = 0;i<(Nfft/2+1);i++)
            Out[i] = (_fftOut[i][0]*_fftOut[i][0]+_fftOut[i][1]*_fftOut[i][1])/Nfft;
    }
    

}

调用的时候，例如我们要执行N点傅里叶变换求PSD，可以首先初始化一个实例或者通过指针生成以供不同对象调用；

MyFFT *FFT_ALg = new MyFFT(N);

然后在每次计算的时候，使用

FFT_ALg->PSD(in ,out);

就可以了，在MyFFT 中，还可以进一步添加幅频响应，相频相应等若干函数。

对于不同点数的FFT需求，需要分别初始化不同点数的实例就可以了。

5. 总结

总而言之，上述面对C++的调用FFTW的使用方案：将初始化，计算与后勤处理隔离开

在连续计算相同点数的FFT应用中，是一种比较省时的用法。