一些FFT库的使用

记录一些FFT库的使用方法

1.kiss_fft

#include "kiss_fftr.h"
#include "kiss_fft.h"
#include "kiss_fftr.h"

#include<stdint.h>

#define WinLen 9
#define N_FFT 9

#include<iostream>
#include<fstream>
using namespace std;

int main()
{
	/* fft */
	kiss_fft_cpx cpx_in[WinLen];
	kiss_fft_cpx cx_out[N_FFT];
	kiss_fft_cfg cfg = kiss_fft_alloc(N_FFT, 0, NULL, NULL);

	for (int i = 0; i < WinLen; i++)
	{
		cpx_in[i].r = i+1;
		cpx_in[i].i = 0;
	}

	kiss_fft(cfg, cpx_in, cx_out);

	for (int i = 0; i < N_FFT; i++)
	{
		cout << cx_out[i].r << "+j"<<cx_out[i].i<<endl;
	}

	free(cfg);

	return 0;

	cout << "+j" << endl;

	/* fftr */
	kiss_fft_scalar cx_in[WinLen];
	kiss_fftr_cfg cfg_r = kiss_fftr_alloc(N_FFT, 0, NULL, NULL);

	for (int i = 0; i < WinLen; i++)
	{
		cx_in[i] = i + 1;
	}

	kiss_fftr(cfg_r, cx_in, cx_out);

	for (int i = 0; i < N_FFT; i++)
	{
		cout << cx_out[i].r << "+j" << cx_out[i].i << endl;
	}
	free(cfg_r);
   
   	/* inverse FFT */
	kiss_fftr_cfg icfg = kiss_fftr_alloc(N_FFT, 1, 0, 0);
	kiss_fft_cpx freq_data[N_FFT / 2 + 1];
	kiss_fft_scalar time_data[N_FFT];

	for (int i = 0; i < N_FFT / 2 + 1; i++)
	{
		freq_data[i].r = cx_out[i].r;
		freq_data[i].i = cx_out[i].i;
	}

	kiss_fftri(icfg, freq_data, time_data);
	for (int i = 0; i < N_FFT; i++)
	{
		cout << "time_data[i]=" << time_data[i]/ N_FFT << endl;
	}
	free(icfg);



}

2. eigen

eigen的fft模块并不常用

	FFT<float> fft;

	Matrix < complex<float>, N_FFT, 1> A;
	MatrixXcf B(N_FFT, 1);  

	for (int i = 0; i < N_FFT; i++)
	{
		A(i) = i+1;
	}

	for (int k = 0; k<A.cols(); ++k)
		fft.fwd(B.col(k),A.col(k));
	cout << "B = " << B << endl;

	for (int k = 0; k<A.cols(); ++k)
		fft.inv(A.col(k), B.col(k));
	cout << "A = " << A << endl;

3.GSL

GSL是一个高效的C数值计算库，文档很详细，看官方使用手册就可以使用
一个简单的实数基-2 FFT例子如下

#include <gsl/gsl_fft_real.h>
int test_radix2_real_fft()
{
    int i; double data[128];

    for (i = 0; i < 128; i++)
      {
         data[i] = 0.0;
      }

    for (i = 1; i <= 10; i++)
      {
         data[i] = 1.0;
      }


    gsl_fft_real_radix2_transform (data, 1, 128);


    return 0;
}

文档中称作为in-place FFT，也就是输出的实部和虚部都是存在原先的数组中，上面的例子中的结果就还是在data数组中，实部跟虚部的存储排列方式如下

complex[0].real    =    data[0]
complex[0].imag    =    0
complex[1].real    =    data[1]
complex[1].imag    =    data[n-1]
...............         ................
complex[k].real    =    data[k]
complex[k].imag    =    data[n-k]
...............         ................
complex[n/2].real  =    data[n/2]
complex[n/2].imag  =    0
...............         ................
complex[k'].real   =    data[k]        k' = n - k
complex[k'].imag   =   -data[n-k]
...............         ................
complex[n-1].real  =    data[1]
complex[n-1].imag  =   -data[n-1]

复数的基-2 FFT使用方式如下：

#include <gsl/gsl_fft_complex.h>

#define REAL(z,i) ((z)[2*(i)])
#define IMAG(z,i) ((z)[2*(i)+1])
int test_complex_fft()
{
    int i; double data[2*128];

    for (i = 0; i < 128; i++)
      {
         REAL(data,i) = 0.0; IMAG(data,i) = 0.0;
      }

    REAL(data,0) = 1.0;

    for (i = 1; i <= 10; i++)
      {
         REAL(data,i) = REAL(data,128-i) = 1.0;
      }

    for (i = 0; i < 128; i++)
      {
        printf ("%d %e %e\n", i,
                REAL(data,i), IMAG(data,i));
      }
    printf ("\n\n");

    gsl_fft_complex_radix2_forward (data, 1, 128);

    for (i = 0; i < 128; i++)
      {
        printf ("%d %e %e\n", i,
                REAL(data,i)/sqrt(128),
                IMAG(data,i)/sqrt(128));
      }

    return 0;
}

关于一些其它的使用方式，如混合基等，官方文档中都有详细的介绍，这里就不赘述了