fftw3/gsl/kissfft/OouraFFT库中傅里叶变换/反傅里叶变换函数和Matlab中的fft/ifft的对应关系

先分析一维度的

一、fftw_plan_dft_1d

正变换：

 fftw_complex *in = fftw_malloc ( sizeof ( fftw_complex ) * n );

 fftw_complex *out = fftw_malloc ( sizeof ( fftw_complex ) * n );

plan_forward = fftw_plan_dft_1d ( n, in, out, FFTW_FORWARD, FFTW_ESTIMATE );

fftw_execute ( plan_forward );

n就是总体个数，不是啥一半还是一半加一那种，in和out都是fftw_complex这种类型

反变换：

 fftw_complex *in2 = fftw_malloc ( sizeof ( fftw_complex ) * n );
plan_backward = fftw_plan_dft_1d ( n, out, in2, FFTW_BACKWARD, FFTW_ESTIMATE );

fftw_execute ( plan_backward );

这时候的in2并不是原来的in，要除以n才可以完全回去和最初的align

plan_backward 的结果需要除以n才可以和 matlab中的 ifft()对应上，而fft不用管，和matlab中的fft是一模一样的

使用完之后，释放内存：

  fftw_destroy_plan ( plan_forward );
  fftw_destroy_plan ( plan_backward );

  fftw_free ( in );
  fftw_free ( in2 );

  fftw_free ( out );

matlab code: fft(a) 和 ifft(a)

二、fftw_plan_dft_r2c_1d 和 fftw_plan_dft_c2r_1d

double *in = fftw_malloc ( sizeof ( double ) * n );

nc = ( n / 2 ) + 1;

fftw_complex* out = fftw_malloc ( sizeof ( fftw_complex ) * nc );

plan_forward = fftw_plan_dft_r2c_1d ( n, in, out, FFTW_ESTIMATE );

fftw_execute ( plan_forward );

c code

 0     55.000000      0.000000
 1     -5.000000     15.388418
 2     -5.000000      6.881910
 3     -5.000000      3.632713

 4     -5.000000      1.624598

5     -5.000000      0.000000

matlab code

fft([1:10]) =
55
-5 +      15.3884176858763i
-5 +      6.88190960235587i
-5 +       3.6327126400268i
-5 +      1.62459848116453i
-5
-5 -      1.62459848116453i
-5 -       3.6327126400268i
-5 -      6.88190960235587i
-5 -      15.3884176858763i

反变换

in2 = fftw_malloc ( sizeof ( double ) * n );

out = fftw_malloc ( sizeof ( fftw_complex ) * nc );

plan_backward = fftw_plan_dft_c2r_1d ( n, out, in2, FFTW_ESTIMATE );

fftw_execute ( plan_backward );

要完全回去还得除以n  

in2[i] / ( double ) ( n )

Gsl 库中有下面几种

这几种是基2的

gsl_fft_real_radix2_transform
gsl_fft_real_float_radix2_transform
gsl_fft_complex_radix2_forward
gsl_fft_complex_float_radix2_forward

这几种是任意尺寸的FFT变换
gsl_fft_real_transform
gsl_fft_real_float_transform
gsl_fft_complex_forward
gsl_fft_complex_float_forward

重点介绍 gsl_fft_complex_radix2_forward （double版本））gsl_fft_complex_float_radix2_forward （float版本）

GSL里面正变换、反变换和 逆变换的差别 具体内容见( 一般都只有正反变换，哪有逆变换，或者把逆变换和反变换默认一个东西）

GSL 学习笔记（快速傅立叶变换）_liyuanbhu的博客-CSDN博客

gsl_fft_complex_radix2_forward  FFT的长度是N的话，那么数组要造2N这么长，

所谓的造数据就是gls库中的 fft的输入输出都是double或者float的数组，它不需要啥gsl_complex 或者gsl_complex_float 这种玩意

下面原型函数中的 gsl_complex_packed_array 和  gsl_complex_packed_array  对应如下：

/* 2N consecutive built-in types as N complex numbers */
typedef double *       gsl_complex_packed_array ;
typedef float *        gsl_complex_packed_array_float  ;

它的数据排列是

data0_real_part  data0_image_part data1_real_part  data2_image_part ....

你需要加的库文件

对于double：gsl_fft_complex.h

对于float：gsl_fft_complex_float.h

原型函数:

GSL_FUN int gsl_fft_complex_radix2_forward (gsl_complex_packed_array data,
                                    const size_t stride,
                                    const size_t n);

GSL_FUN int gsl_fft_complex_float_radix2_forward (gsl_complex_packed_array_float data,
                                          const size_t stride,
                                          const size_t n);

执行gsl_fft_complex_radix2_forward 之后，结果还保存在 data里面（数据会被冲掉！）

下面是原型：

int gsl_fft_complex_radix2_forward (gsl_complex_packed_array data, size t stride, size t n)
int gsl_fft_complex_radix2_backward (gsl_complex_packed_array data, size t stride, size t n)
int gsl_fft_complex_radix2_inverse (gsl_complex_packed_array data, size t stride, size t n)

下面是使用方法：

gsl_fft_complex_radix2_forward(data, 1, 128);
gsl_fft_complex_radix2_backward(data, 1, 128);
gsl_fft_complex_radix2_inverse(data, 1, 128);

backward是不除以N, 和matlab ifft 对比，需要matlab里面的ifft结果乘以 N

forward~fft 正变换

inverse~ifft 逆变换

他们和matlab的 fft ifft 等效

FFTW3中的二维变换的例子：【二维变换就是按行优先码成一维向量】

typedef struct _COMPLEX
{
    double real;
    double img;
}COMPLEX;

正变换

void FFT2(double **input, COMPLEX **output, int height, int width)
{
    fftw_plan planR;
    fftw_complex *inR, *outR;
    inR = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * height * width);
    outR = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * height * width);
    for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
        {
            inR[i * width + j][0] = input[i][j];
            inR[i * width + j][1] = 0;
        }
    planR = fftw_plan_dft_2d(height, width, inR, outR, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(planR);
    for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
        {
            output[i][j]._COMPLEX::real = outR[i * width + j][0];
            output[i][j]._COMPLEX::img = outR[i * width + j][1];
        }
    fftw_destroy_plan(planR);
    fftw_free(inR);
    fftw_free(outR);
}

反变换

void IFFT2(COMPLEX **input, double **output, int height, int width)
{
    fftw_plan planR;
    fftw_complex *inR, *outR;
    inR = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * height * width);
    outR = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * height * width);
    for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
        {
            inR[i * width + j][0] = input[i][j]._COMPLEX::real;
            inR[i * width + j][1] = input[i][j]._COMPLEX::img;
        }
    planR = fftw_plan_dft_2d(height, width, inR, outR, FFTW_BACKWARD, FFTW_ESTIMATE);
    fftw_execute(planR);
    for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
        {
            output[i][j] = outR[i * width + j][0] / (height*width);
        }
    fftw_destroy_plan(planR);
    fftw_free(inR);
    fftw_free(outR);
}

gsl 貌似没有2D的fft和ifft？ （待确认。。。）

下面是一个使用fftw来搞定gsl 2D fft的一个例子

#include <fftw3.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>

int main()
{
    //Create 2 GSL Complex matrices
    gsl_matrix_complex* mc = gsl_matrix_complex_alloc(3, 3);
    gsl_matrix_complex* mc2 = gsl_matrix_complex_alloc(3, 3);
    int i, j;

    /* FFTW requires data stored in row major order. Real and Imaginary parts of complex number stored as double[2]
    Fortunately, GSL implements gsl_matrix_complex in the required format. */
    fftw_complex* in = (fftw_complex*)mc->data;
    fftw_complex* out = (fftw_complex*)mc2->data;

    fftw_plan fp = fftw_plan_dft_2d(3, 3, in, in, FFTW_FORWARD, FFTW_MEASURE);
    fftw_plan fp2 = fftw_plan_dft_2d(3, 3, in, in, FFTW_BACKWARD, FFTW_MEASURE);

    for (i = 0; i < 3; i++)
    {
        for (j = 0; j < 3; j++)
        {
            gsl_matrix_complex_set(mc, i, j, gsl_complex_rect(i, j));
        }
    }

    //Calculate FFT
    fftw_execute(fp);

    for (i = 0; i < 3; i++)
    {
        for (j = 0; j < 3; j++)
        {
            printf("%f %f\n", GSL_REAL(gsl_matrix_complex_get(mc, i, j)), GSL_IMAG(gsl_matrix_complex_get(mc, i, j)));
        }
    }

    printf("\n");

    gsl_matrix_complex_scale(mc, gsl_complex_rect((1 / 9.0), 0));
    fftw_execute(fp2);

    for (i = 0; i < 3; i++)
    {
        for (j = 0; j < 3; j++)
        {
            printf("%f %f\n", GSL_REAL(gsl_matrix_complex_get(mc, i, j)), GSL_IMAG(gsl_matrix_complex_get(mc, i, j)));
        }
    }
    return 0;

}
 

再加上kissfft的一些用法吧

typedef struct {
    kiss_fft_scalar r;
    kiss_fft_scalar i;
}kiss_fft_cpx;
自定义类一个结构体，包含实部和虚部

kiss_fft_scalar 就是float或者double

kiss_fftr_cfg kiss_fftr_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem)

inverse_fft 0 正变换

inverse_fft 1 反变换，后面 mem lenmem 都搞NULL

正变换，由 实数 变换到 复数

void kiss_fftr(kiss_fftr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata)

反变换，由 复数 变换到 实数

void kiss_fftri(kiss_fftr_cfg st,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata)

数组分配情况

kiss_fft_scalar cx_in[FFT_LEN];

kiss_fft_cpx freq_data[FFT_LEN / 2 + 1];

OouraFFT库 中的

cdft (complex DFT），数据按实部虚部挨个顺序给如

rdft  （real DFT），只是实数fft

虚部差个负号