如何使用开源FFTW库进行FFT和逆FFT运算

1. 首先下载FFTW3库，在Windows平台上使用Developer Command Prompt工具生成静态库.lib文件。
2. 创建VC测试项目，将库中提供的libfftw3-3.dll和fftw3.h，以及生成的libfftw3-3.lib放在测试项目中，配置好include path和link path/item。

// FFT_CPP_Test.cpp : This file contains the 'main' function. Program execution begins and ends there.
//

#include <iostream>
#include "fftw3.h"

void PrintFFT(fftw_complex *val, size_t size) {
   std::cout << "===============" << std::endl;
   for (int i = 0; i < size; ++i) {
      std::cout << val[i][0] << '\t' << val[i][1] << std::endl;
   }
}

int main()
{
   fftw_complex *in, *fftout, *ifftout;
   fftw_plan forwardPlan, backwordPlan;

   int n = 5;
   in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * n);
   fftout = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * n);
   ifftout = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * n);

   for (int i = 0; i < n; ++i) {
      in[i][0] = i;
      in[i][1] = 0;
   }
   std::cout << "Origin: " << std::endl;
   PrintFFT(in, n);

   forwardPlan = fftw_plan_dft_1d(5, in, fftout, FFTW_FORWARD, FFTW_ESTIMATE);
   backwordPlan = fftw_plan_dft_1d(5, fftout, ifftout, FFTW_BACKWARD, FFTW_ESTIMATE);

   fftw_execute(forwardPlan); // FFT
   // normalization
   for (int i = 0; i < n; ++i) {
      fftout[i][0] /= n;
      fftout[i][1] /= n;
   }
   std::cout << "After FFT: " << std::endl;
   PrintFFT(fftout, n);

   fftw_execute(backwordPlan); // Inverse FFT
   std::cout << "Inverse FFT: " << std::endl;
   PrintFFT(ifftout, n);

   fftw_destroy_plan(forwardPlan);
   fftw_destroy_plan(backwordPlan);
   fftw_free(in);
   fftw_free(fftout);
   fftw_free(ifftout);

   return 0;
}

测试结果：

Origin:
===============
0       0
1       0
2       0
3       0
4       0
After FFT:
===============
2       0
-0.5    0.688191
-0.5    0.16246
-0.5    -0.16246
-0.5    -0.688191
Inverse FFT:
===============
0       0
1       0
2       0
3       0
4       0

reference：https://stackoverflow.com/questions/39675436/how-to-get-fftw-working-on-windows-for-dummies

数据经过FFT后，在计算频谱之前会进行shift，FFT shift代码如下：

void swap(complex *v1, complex *v2)
{
    complex tmp = *v1;
    *v1 = *v2;
    *v2 = tmp;
}

void fftshift(complex *data, int count)
{
    int k = 0;
    int c = (int) floor((float)count/2);
    // For odd and for even numbers of element use different algorithm
    if (count % 2 == 0)
    {
        for (k = 0; k < c; k++)
            swap(&data[k], &data[k+c]);
    }
    else
    {
        complex tmp = data[0];
        for (k = 0; k < c; k++)
        {
            data[k] = data[c + k + 1];
            data[c + k + 1] = data[k + 1];
        }
        data[c] = tmp;
    }
}

void ifftshift(complex *data, int count)
{
    int k = 0;
    int c = (int) floor((float)count/2);
    if (count % 2 == 0)
    {
        for (k = 0; k < c; k++)
            swap(&data[k], &data[k+c]);
    }
    else
    {
        complex tmp = data[count - 1];
        for (k = c-1; k >= 0; k--)
        {
            data[c + k + 1] = data[k];
            data[k] = data[c + k];
        }
        data[c] = tmp;
    }
}

Reference: https://stackoverflow.com/questions/5915125/fftshift-ifftshift-c-c-source-code

通常会使用多线程的方式，同时对多路数据进行FFT计算，这时会出现fftw_dft_plan read violation的问题，为了避免冲突，可以对plan加锁：

#include <mutex>

std::mutex m;

m.lock(); // fftw plan read voilation when multi thread begin
forwardPlan = fftw_plan_dft_1d(n, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
m.unlock();
fftw_execute(forwardPlan);

Reference: https://www.cnblogs.com/zizbee/p/13520823.html