FFTW cuFFT的使用记录

好久没写点什么东西了，今天饶有兴趣，总结一下FFTW、cuFFT的调用方法。

一些知识点的回顾

弄懂了FT、DTFT、DFT三者之间的关系

傅里叶变换（Fourier Transform，FT），表示能够将一定条件的某个函数表示为三角函数或者它们的积分的线性组合。从连续域到连续域。

离散时间傅里叶变换（Discrete-time Fourier Transform，DTFT），表示将以离散时间nT作为变量的函数变换到连续频域，频谱会周期延拓。

离散傅里叶变换（Discrete Fourier Transform，DFT），将时域信号的采样变换为在DTFT频域的采样。形式上，两端都是有限长，实际上，两组序列都应该被认为是离散周期信号的主值序列。即使对有限长的离散信号做DFT，也应该将其看做经过周期延拓成为周期信号再做变换。

说到底，还是为了让计算机可以处理傅里叶变换，因此需要在时域和频域都是离散信号。

弄懂了如何从离散傅里叶级数到离散傅里叶变换，其中包含公式推导

将要处理的离散非周期信号看做是一个离散周期信号的主值序列，然后对其做DFS，取结果的主值序列（DFS的结果也是周期的）。

循环卷积的计算步骤：

1. 有限长序列构造周期序列
2. 计算周期卷积
3. 周期卷积取主值

循环卷积需要两个序列的长度一样，也就是周期一样。

使用循环卷积来计算线性卷积：两个序列长度分别为N,M，那么取L>max(N,M)作为两个序列的周期，在两个序列后面补零，然后计算周期卷积，取主值序列就是循环卷积。

从DFT到FFT：

1. 分奇偶；
2. DFT是DFS的主值序列，本身实质上是隐含周期，因此可以快算计算出另一半结果；
3. 于是就将算法复杂度转为了对数复杂度；

之前网上找的一个FFT实现，现在看，更加透彻了

1. 先对输入信号的位置做一个反转，此处用到了位运算和动态规划；
2. 一排一排的计算，文中的h代表当前处理长度的一半；

FFTW实现

FFTW是一个开源的fft库，首先需要去官网下载压缩包windows版本，然后使用VS Command Prompt工具生成lib，首先cd到压缩包的目录下，然后执行命令lib /machine:x64 /def:libfftw3l-3.def，注意，此处有一个大坑：linux下的可以直接cd到指定目录，而在win下，竟然需要先到D盘，然后cd。

基本流程

1. 创建句柄plan，可以将其看做一个对象，里面需要包含输入in，输出out的地址；
2. 初始化数据
3. 执行
4. 销毁

#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); /* repeat as needed */
...
fftw_destroy_plan(p);
fftw_free(in); fftw_free(out);
}

简单测试

做一个8点的FFT，使用实数形式和复数形式：

#include<iostream>

#include"fftw3.h"
#pragma comment(lib, "libfftw3-3.lib") // double版本
// #pragma comment(lib, "libfftw3f-3.lib")// float版本
// #pragma comment(lib, "libfftw3l-3.lib")// long double版本


const double PI = acos(-1); //很巧妙的一种定义pi的方式


int main() {
	int len = 8;

	// 如果要使用float版本,需先引用float版本的lib库,然后在fftw后面加上f后缀即可.
	double *in = NULL;
	fftw_complex *out = NULL;// fftwf_complex --> 即为float版本
	fftw_complex *in2 = NULL;
	fftw_complex *out2 = NULL;// fftwf_complex --> 即为float版本
	fftw_plan p, p2;

	//分配内存空间
	in = (double *)fftw_malloc(sizeof(double) * len);
	out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * len);
	in2 = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * len);
	out2 = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * len);


	// 创建句柄
	p = fftw_plan_dft_r2c_1d(len, in, out, FFTW_ESTIMATE);
	p2 = fftw_plan_dft_1d(len, in2, out2, FFTW_FORWARD, FFTW_ESTIMATE);

	// 输入
	printf("in:\n");
	for (int i = 0; i < len; i++)
	{
		//in[i] = sin(2 * PI * dx*i) + sin(4 * PI * dx*i);
		in[i] = i + 1;
		printf("%.2f ", in[i]);
	}
	printf("\n\n");

	printf("in2:\n");
	for (int i = 0; i < len; i++)
	{
		in2[i][0] = i + 1;
		in2[i][1] = 0;
		printf("%.2f %.2f\n", in2[i][0], in2[i][1]);
	}
	printf("\n\n");

	//执行
	fftw_execute(p);
	fftw_execute(p2);

	// 输出
	printf("out1:\n");
	for (int i = 0; i < len; i++)
	{
		printf("%.5f ,%.5f  \n", out[i][0], out[i][1]);
	}
	printf("\n");

	printf("out2:\n");
	for (int i = 0; i < len; i++)
	{
		printf("%.5f ,%.5f  \n", out2[i][0], out2[i][1]);
	}
	printf("\n");

	// 释放资源
	fftw_destroy_plan(p);
	fftw_free(in);
	fftw_free(out);

	fftw_destroy_plan(p2);
	fftw_free(in2);
	fftw_free(out2);

	//system("pause");
	return 0;

}

对于一个数据文件，需要对其做FFT，应该是一段一段的做，得到的其实是短时傅里叶变换，单纯的傅里叶变换频谱不能反映出时间信息，短时傅里叶变换可以体现出这一点；

cuFFT实现

cuFFT的思想和FFTW基本一致，唯一有点区别的是，cuFFT需要将数据从cpu拷贝到gpu，然后计算，最后从gpu拷贝回cpu;

1. 创建句柄，分配内存
2. 初始化，内存拷贝cpu->gpu
3. 执行
4. 内存拷贝gpu->cpu
5. 销毁

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "cufft.h"


#include<iostream>


#include"fftw3.h"
#pragma comment(lib, "libfftw3-3.lib") // double版本
// #pragma comment(lib, "libfftw3f-3.lib")// float版本
// #pragma comment(lib, "libfftw3l-3.lib")// long double版本



#define CHECK(call)\
{\
	if ((call) != cudaSuccess)\
			{\
		printf("Error: %s:%d, ", __FILE__, __LINE__);\
		printf("code:%d, reason: %s\n", (call), cudaGetErrorString(cudaGetLastError()));\
		exit(1);\
			}\
}


const double PI = acos(-1); //很巧妙的一种定义pi的方式

void test_FFTW();

int main() {

	const int NX = 8;
	const int BATCH = 1;

	cufftHandle plan;
	cufftComplex *data;
	cufftComplex *data_cpu;

	data_cpu = (cufftComplex *)malloc(sizeof(cufftComplex) * NX * BATCH);
	if (data_cpu == NULL) return -1;

	CHECK(cudaMalloc((void**)&data, sizeof(cufftComplex) * NX * BATCH));
	
	CHECK(cufftPlan1d(&plan, NX, CUFFT_C2C, BATCH)); 

	//输入数据
	for (int i = 0; i < NX; ++i) {
		data_cpu[i].x = i + 1;
		data_cpu[i].y = 0;
	}

	//数据传输cpu->gpu
	CHECK(cudaMemcpy(data, data_cpu, sizeof(cufftComplex) * NX * BATCH, cudaMemcpyHostToDevice));
	CHECK(cudaDeviceSynchronize());
	
	CHECK(cufftExecC2C(plan, data, data, CUFFT_FORWARD)); 
	//CHECK(cufftExecC2C(plan, data, data, CUFFT_INVERSE) != CUFFT_SUCCESS);

	//数据传输gpu->cpu
	CHECK(cudaMemcpy(data_cpu, data, sizeof(cufftComplex) * NX * BATCH, cudaMemcpyDeviceToHost));
	CHECK(cudaDeviceSynchronize());
	

	cufftDestroy(plan);
	cudaFree(data);

	printf("CUFFT_FORWARD:\n");
	for (int i = 0; i < NX; ++i) {
		printf("%f , %f\n", data_cpu[i].x, data_cpu[i].y);
	}

	system("pause");
	return 0;

}