fftw的使用

1、下载编译

官网：http://www.fftw.org/index.html
官方文档：http://www.fftw.org/fftw3_doc/

数据类型： FFTW 有三个版本的数据类型:double、float 和 long double，使用方法如下:

1. 链接对应的库(比如 libfftw3-3、libfftw3f-3、或 ibfftw3l-3)；
2. 包含同样的头文件 fftw3.h；将所有以小写"fftw_“开头的名字替换为"fftwf_”(float 版本)或"fftwl_"(long double 版本)。比如将 fftw_complex 替换为 fftwf_complex,将 fftw_execute替换为 fftwf_execute 等。
3. 所有以大写"FFTW_"开头的名字不变
4. 将函数参数中的 double 替换为 float 或 long double
5. 最后,虽然 long double 是 C99 的标准,但你的编译器可能根本不支持该类型,或它并不能提供比 double 更高的精度。
6. fftw_malloc 考虑了数据对齐,以便使用 SIMD 指令加速,所以最好不要用 C 函数malloc 替代,而且不要将 fftw_malloc、fftw_free 和 malloc、free、 delete 等混用。
   尽量使用 fftw_malloc 分配空间,而不是使用的静态数组,因为静态数组是在栈上分配的,而栈空间较小;还因为这种方式没有考虑数据对齐,不便应用SIMD 指令。

编译构造：

1. 默认FFTW编译生成double类型，加入参数“–enable-single”或“–enable-float”编译单精度(float)，加入参数“–enable-long-double”支持长双进度类型
2. 另外,在ARM平台上可增加“–enable-neon”,使能ARM NEON流媒体加速核心

1.1 Linux下编译

fftw 编译安装说明

./configure --prefix=/home/xxx/usr/fftw --enable-mpi --enable-openmp --enable-threads --enable-shared MPICC=mpicc CC=gcc F77=gfortran
make
make install

可通过./configure --help命令查看具体配置：

–prefix：设定安装目录
–enable-mpi：是否编译mpi版的fftw库
–enable-openmp：是否使用OpenMP指令进行并行
–enable-threads：是否编译FFTW SMP线程库
–enable-shared：是否编译动态库
MPICC：指定C语言的MPI编译器
CC：指定C语言的编译器
F77：指定Fortran 77的编译器

2、FFT基础知识

2.1 概念

• FFT分辨率可以表示为：f_s/N_fft
  频率分辨率的物理量就是：观测信号的时间窗长度， 时间窗越长（N大）， 对应频率分辨率就越高。
  提高频率分辨率的两种方法：增加fft点数N（计算量大），降低采样率（时间分辨率降低了）。

• 输出取N/2个频点，N点fft的输出值也是N个（对称的），故只取一半。

• 时域上的补0相当于频域上的插值，通过补0增加的fft点数无法提高fft精度

• 幅度、相位
  

• 示例：fs =1024/10 N=1024 △f=0.1Hz
  
  

2.2 对1024个点的信号，做4次256点FFT和1次1024点FFT的区别

有信号如下:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

如果拆分成4个4点FFT,或者1个16点,其中的数学关系不能凭直觉发现,但是,如果把这段信号拆分成:

1 2 3 4 | 0 0 0 0 | 0  0  0  0 | 0  0  0  0
0 0 0 0 | 5 6 7 8 | 0  0  0  0 | 0  0  0  0
0 0 0 0 | 0 0 0 0 | 9 10 11 12 | 0  0  0  0
0 0 0 0 | 0 0 0 0 | 0  0  0  0 |13 14 15 16

然后做4次16点的FFT,那么这4次16点的FFT频谱叠加,最后得到的叠加后的频谱就=一次16点的FFT.
那么,问题就变成了: 1个带12个0的16点的FFT跟1个4点的FFT之间是什么关系?
1 2 3 4 0 0 0 0 0 0 0 0 0 0 0 0和1 2 3 4,他们的频谱,就相当于1234的4点频谱做插值成16个点.
至于0 0 0 0 5 6 7 8 0 0 0 0 0 0 0 0和5 6 7 8 之间的关系,你知道根据时移定理,对频谱乘一个e^-i2pikm相当于时域上延迟k个点:0 0 0 0 5 6 7 8 0 0 0 0 0 0 0 0可以变形为5 6 7 8 0 0 0 0 0 0 0 0 0 0 0 0
也就是说,4点4FFT与1点的16FFT的关系,大致相当于:

FFT(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16)=插值(FFT(1,2,3,4))+插值(FFT(5,6,7,8))*时移+插值(FFT(9,10,11,12))*时移^2+插值(FFT(13,14,15,16))*时移^3

3、FFTW基础

全球最快的傅里叶变换算法（FFTW）

3.1 基本函数

plan= fftw_plan_dft_1d(_N,_in,_out,FFTW_BACKWARD,FFTW_MEASURE);//配置计算计划
倒数第二个参数sign， FFTW_FORWARD：正变换 ；FFTW_BACKWARD：逆变换。
最后的一个参数flags，FFTW_MEASURE：表示先计算一些FFT并测量所用的时间，以便为大小为n的变换寻找最优的计算方法。FFTW_MEASURE模式下in和out数组中的值会被覆盖，所以该方式的plan应该在用户初始化输入数据in之前完成。

//c-c
fftw_plan fftw_plan_dft_1d(int n, fftw_complex *in, fftw_complex *out, int sign, unsigned flags);
fftw_plan fftw_plan_dft_2d(int n0, int n1, fftw_complex *in, fftw_complex *out, int sign, unsigned flags);
fftw_plan fftw_plan_dft_3d(int n0, int n1, int n2, fftw_complex *in, fftw_complex *out, int sign, unsigned flags);
fftw_plan fftw_plan_dft(int rank, const int *n, fftw_complex *in, fftw_complex *out, int sign, unsigned flags);
//r-c
fftw_plan fftw_plan_dft_r2c_1d(int n, double *in, fftw_complex *out, unsigned flags);
fftw_plan fftw_plan_dft_r2c_2d(int n0, int n1, double *in, fftw_complex *out, unsigned flags);
fftw_plan fftw_plan_dft_r2c_3d(int n0, int n1, int n2, double *in, fftw_complex *out, unsigned flags);
fftw_plan fftw_plan_dft_r2c(int rank, const int *n, double *in, fftw_complex *out, unsigned flags);
//r-r
fftw_plan fftw_plan_r2r_1d(int n, double *in, double *out, fftw_r2r_kind kind, unsigned flags);
fftw_plan fftw_plan_r2r_2d(int n0, int n1, double *in, double *out, fftw_r2r_kind kind0, fftw_r2r_kind kind1, unsigned flags);
fftw_plan fftw_plan_r2r_3d(int n0, int n1, int n2, double *in, double *out, fftw_r2r_kind kind0, fftw_r2r_kind kind1, fftw_r2r_kind kind2, unsigned flags);
fftw_plan fftw_plan_r2r(int rank, const int *n, double *in, double *out, const fftw_r2r_kind *kind, unsigned flags);

void fftw_execute(const fftw_plan p);
void fftw_execute_dft(const fftw_plan p, fftw_complex *in, fftw_complex *out);     
void fftw_execute_split_dft(const fftw_plan p, double *ri, double *ii, double *ro, double *io);     
void fftw_execute_dft_r2c(const fftw_plan p,double *in, fftw_complex *out);     
void fftw_execute_split_dft_r2c(const fftw_plan p, double *in, double *ro, double *io);
void fftw_execute_dft_c2r(const fftw_plan p, fftw_complex *in, double *out);
void fftw_execute_split_dft_c2r(const fftw_plan p, double *ri, double *ii, double *out);
void fftw_execute_r2r(const fftw_plan p, double *in, double *out);

参数sign可以为正变换FFTW_FORWARD(-1)，也可以为逆变换FFTW_BACKWORD(+1)，实际上就是变换公式中指数项的符号。需注意FFTW的逆变换没有除以N，即数据正变换再反变换后是原始数据的N倍。
单序列正反变换的归一化处理时：除以N/2

3.2数据类型

FFTW有三个版本的数据类型：double、float、long double，
使用方法如下：

• 链接对应的库（比如libfftw3-3、libfftw3f-3、或ibfftw3l-3） 包含同样的头文件fftw3.h
  单精度 前缀 “-fftwf” 编译选项 “-lfftw3f”
  双精度 前缀 “-fftwl” 编译选项 “-lfftw3l”
• 将所有以小写"fftw_“开头的名字替换为"fftwf_”（float版本）或"fftwl_"（long double版本）。
  比如将fftw_complex替换为fftwf_complex，
  将fftw_execute替换为fftwf_execute等。
• 所有以大写"FFTW_"开头的名字不变，将函数参数中的double替换为float或long double

最后，虽然long double是C99的标准，但你的编译器可能根本不支持该类型，或它并不能提供比double更高的精度。

3.3 计算公式

官方计算公式介绍

3、示例

Qt中用fftw对音频数据变换
全球最快的傅里叶变换算法（FFTW）

#include <fftw3.h>
     ...
{
    fftw_complex *in, *out;
    fftw_plan p;
    ...
    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);    
	...
    p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);	//变换类型
    fftw_execute(p); 				// 执行变换
    ...
    fftw_destroy_plan(p);
    fftw_free(in); 
    fftw_free(out);
}

int main()
{
	int Nfft = 1024, Ns = 600;
	double fc1 = 0e6;
	double fc2 = 1.2e6;
	double Fs = 3.5e6;
	fftw_complex *fftIn = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*Nfft);
	fftw_complex *fftOut = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*Nfft);
	fftw_plan plan = fftw_plan_dft_1d(Nfft, fftIn, fftOut, FFTW_FORWARD, FFTW_MEASURE);
	for (int i = 0; i < Ns; i++)
	{
		fftIn[i][0] = cos(2 * PI*fc1*i / Fs)+ cos(2 * PI*fc2*i / Fs); // real part
		fftIn[i][1] = sin(2 * PI*fc1*i / Fs)+ sin(2 * PI*fc2*i / Fs); // imag part
	}
	fftw_execute(plan);
	fftw_destroy_plan(plan);
	double *fftres = new double[Nfft];
	VectorXd freqpoint = VectorXd::LinSpaced(Nfft, -Nfft / 2, Nfft / 2 - 1);
	freqpoint = freqpoint*Fs / Nfft/1e6;
	for (int i = 0; i < Nfft; i++)
	{
		if (i<Nfft/2)
			fftres[i + Nfft / 2] = 10 * log10(sqrt(fftOut[i][0] * fftOut[i][0] + fftOut[i][1] * fftOut[i][1]));
		else
			fftres[i - Nfft / 2] = 10 * log10(sqrt(fftOut[i][0] * fftOut[i][0] + fftOut[i][1] * fftOut[i][1]));
	}
}

音频变换的例子：

	int N = filesize/2;     			//计算数据个数
	short *in_buf;						//自己存输入数据
	//为fft输入计算分配空间
	double * in = (double*)fftw_malloc(sizeof(double) * N);
	for(int i=0; i<N; i++)
	{
	   in[i] = in_buf[i]; 			 	//将pcm文件中的数据复制到fft的输入
	}
	
	//为fft输出算分配空间
	fftw_complex * out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);
	
	//进行fft变换，fftw_plan_dft_c2r_1d函数进行反变换
	fftw_plan p = FFTW3_H::fftw_plan_dft_r2c_1d(N, in, out, FFTW_ESTIMATE);
	fftw_execute(p);
	
	double dx3 = (double)SAMPLE_RATE / N;
	
	polt_1->xAxis->setRange(0, SAMPLE_RATE/2, Qt::AlignLeft);
	
	//根据FFT计算的复数计算振幅谱
	for( int i=0; i<N/2; i++ )
	{
	   double val = sqrt(out[i][0] * out[i][0] + out[i][1] * out[i][1]);
	   val = val / (N / 2);
	   polt_1->graph(0)->addData( dx3 * i, val );
	
	   if( val > val_max )
	   {
	       val_max = val;
	   }
	
	   double db = log(val);
	   //qDebug("frequency = %f, amplitude = %f, db = %f", dx3 * i, val / (N / 2), db);
	}
	
	polt_1->yAxis->setRange(val_max*0.6, val_max*1.2, Qt::AlignBottom);
	polt_1->replot();
	
	//根据FFT计算的复数计算相位谱
	polt_2->xAxis->setRange(0, SAMPLE_RATE/2, Qt::AlignLeft);
	polt_2->yAxis->setRange(0, 10, Qt::AlignBaseline);	
	for( int i=0; i<N/2; i++ )
	{
	   double val = atan2(out[i][1], out[i][0]);
	   polt_2->graph(0)->addData( dx3 * i, val );
	}
	polt_2->replot();
	
	//fftw销毁
	fftw_destroy_plan(p);	
	fftw_free(in);
	fftw_free(out);

N点的FFT，以 “采样率/N” 为频率间隔，

4、注意事项

C++ FFTW3的fftw_plan_dft_r2c_2d()和Eigen::MatrixXd的傅里叶正变换