目录
- 任务介绍
- 环境所需相关软件下载与安装
- C语言：不调用库的GPU加速FFT代码
- C语言：调用fftw库的未使用GPU的FFT代码
- C语言：调用cufft库的GPU加速FFT
- gnuplot安装画图，maltab编写的FFT运算结果对比
- matlab测试信号和测试时的坑
任务介绍
时隔多年仍然逃不掉写C的命运……因为这个任务周期不短还踩了好多坑，必须记录一下了。
任务简单要求就是使用C语言编写一个GPU加速的快速傅里叶变换（FFT）
分为GPU加速的FFT代码改写、未使用GPU的FFT编写、运算速度对比、运算结果测试（与matlab结果对比），只要按照我文章写的顺序做就行
环境所需相关软件下载与安装
Visual Studio 2010
要运行C语言代码就要先下载Microsoft Visual Studio 编辑器，我的电脑是Win10系统，但为了与项目环境和大家的使用习惯保持一致，使用的是Visual Studio 2010版本，网上有安装包的下载，这里提供一个我保存的中文版安装包
链接：https://pan.baidu.com/s/12Eyw5Woh11gtP6hyxbpUmQ
提取码：em9h
安装过程还挺顺利的，网上也有很多安装教程，就是安装时间挺久，我安装在D盘，最后应用图标在”D:\Visual Studio 2010\Common7\IDE\devenv.exe”，安装在其他路径的只要找到IDE文件夹就能找到。
CUDA7.5
然后需要下载CUDA 来调用GPU。首先保证自己的电脑是有显卡GPU的（GPU和CUDA版本对应关系网上可以找到，这里插一个写挺好的文章链接CUDA和GPU版本对应），我的笔记本是双显卡，独立显卡是NVIDIA GeForce GTX 1660Ti，最高能用CUDA11.7版本的
集成显卡：为电脑基本显卡，在正常情况下电脑会一直使用集成显卡，以降低电脑运行负担，减少散热，提高笔记本使用寿命。
独立显卡：为电脑的高级显卡，多用来运行大型游戏以及部分软件。
要改为优先使用独立显卡的话，就打开NIVDIA控制面板，在左面管理3D设置中全局设置->首选图形处理器->下拉选择高性能NIVDIA处理器，就设置成功了
CUDA和Visual Studio 2010是有版本兼容关系的，在安装好Visual Studio 2010后安装CUDA需要选择对应版本，网上找到Visual Studio 2013是安装的CUDA7.5，但我想装一个能兼容的最高级版，就从9.0试到7.5 （CUDA历史版本下载链接），要是版本不兼容在安装CUDA的时候会出现不兼容的提示，就没办法继续安装，那个图我没保存，但安装的时候看到就能明白（安装教程简单，最好安装在默认路径，网上可找到教程），最终确定兼容版本是CUDA7.5+Visual Studio 2010
环境配置
软件安装完后配置环境，这里我基本是按照这篇文章教程配置的
windows+VS2013+CUDA7.5配置
只有在这篇文章中的第3步我创建项目的方式不太一样
我是打开VS2010创建一个空项目，在源文件中添加一个新建项，选择C++文件，这里注意：不加后缀名生成的源文件是.cpp后缀名的，这是C++代码文件后缀名，使用CUDA的C代码要命名为.cu的格式，不使用CUDA的C代码命名为.c格式
如果在这篇文章第3.f步的项类型里没有找到CUDA C++的话，就右键 【工程】-【生成自定义】-勾选上CUDA 7.5就行，如下图
最后最好拿他给的测试代码试跑一下，能成功就算环境配置成功
C语言：不调用库的GPU加速FFT代码
开始网上找到这篇文章CUDA实现FFT并行计算
在我简单调试更改后可以正常运行，很开心，我对代码加了些注释，这里把这个能正常运行的代码贴上，分为两个部分代码。
先是需要调用的函数Complex.cu文件，我跟主文件放在一起
#define PI 3.1415926class Complex {public:double real;double imag;Complex() {}// Wn 获取n次单位复根中的主单位根__device__ static Complex W(int n) {Complex res = Complex(cos(2.0 * PI / n), sin(2.0 * PI / n));return res;}// Wn^k 获取n次单位复根中的第k个__device__ static Complex W(int n, int k) {Complex res = Complex(cos(2.0 * PI * k / n), sin(2.0 * PI * k / n));return res;}// 实例化并返回一个复数（只能在Host调用）static Complex GetComplex(double real, double imag) {Complex r;r.real = real;r.imag = imag;return r;}// 随机返回一个复数static Complex GetRandomComplex() {Complex r;r.real = (double)rand() / rand();r.imag = (double)rand() / rand();return r;}// 随即返回一个实数static Complex GetRandomReal() {Complex r;r.real = (double)rand() / rand();r.imag = 0;return r;} // // 随即返回一个纯虚数 // static Complex GetRandomPureImag() //{ // Complex r; // r.real = 0; // r.imag = (double)rand() / rand(); // return r; // }// 构造函数（只能在Device上调用）__device__ Complex(double real, double imag) {this->real = real;this->imag = imag;}// 运算符重载__device__ Complex operator+(const Complex &other) {Complex res(this->real + other.real, this->imag + other.imag);return res;}__device__ Complex operator-(const Complex &other) {Complex res(this->real - other.real, this->imag - other.imag);return res;}__device__ Complex operator*(const Complex &other) {Complex res(this->real * other.real - this->imag * other.imag, this->imag * other.real + this->real * other.imag);return res;}};然后是主文件
// 一维FFT#include "cuda_runtime.h"#include "device_launch_parameters.h"#include "Complex.cu"// 自定义的复数数据结构#include #include #include #include #include #include#pragma comment(lib,"winmm.lib")using namespace std;// 根据数列长度n获取二进制位数int GetBits(int n) {int bits = 0;while (n >>= 1) {bits++;}return bits;}// 在二进制位数为bits的前提下求数值i的二进制逆转__device__ int BinaryReverse(int i, int bits) {int r = 0;do {r += i % 2 << --bits;} while (i /= 2);return r;}// 蝴蝶操作, 输出结果直接覆盖原存储单元的数据, factor是旋转因子__device__ void Bufferfly(Complex *a, Complex *b, Complex factor) {Complex a1 = (*a) + factor * (*b);Complex b1 = (*a) - factor * (*b);*a = a1;*b = b1; }__global__ void FFT(Complex nums[], Complex result[], int n, int bits) {int tid = threadIdx.x + blockDim.x * blockIdx.x;//threadIdx,线程ID的索引//blockDim:线程块的维度//blockIdx,线程块ID的索引if (tid >= n) return;for (int i = 2; i < 2 * n; i *= 2){if (tid % i == 0) {int k = i;if (n - tid < k) k = n - tid;for (int j = 0; j < k / 2; ++j) {Bufferfly(&nums[BinaryReverse(tid + j, bits)], &nums[BinaryReverse(tid + j + k / 2, bits)], Complex::W(k, j));}}__syncthreads();//未定义标识符}result[tid] = nums[BinaryReverse(tid, bits)];}// 打印数列void printSequence(Complex nums[], const int N) {printf("[");for (int i = 0; i < N; ++i) {double real = nums[i].real, imag = nums[i].imag;if (imag == 0) printf("%.16f", real);else {if (imag > 0) printf("%.16f+%.16fi", real, imag);else printf("%.16f%.16fi", real, imag);}if (i != N - 1) printf(", ");}printf("]\n");}int main() {//srand(time(0));// 设置当前时刻为随机数种子const int TPB = 1024;// 每个Block的线程数，即blockDim.x，每个块中的线程数量 <= 1024const int N = 4000000;// 数列大小const int bits = GetBits(N);DWORD Start, End;// 随机生成实数数列Complex *nums = (Complex*)malloc(sizeof(Complex) * N), *dNums, *dResult;//申请一块内存，其大小是N个complex长度的总和。然后把这块内存的首地址强转成complex *指针变量类型，赋给*numsfor (int i = 0; i < N; ++i) { nums[i] = Complex::GetRandomReal();//没有在类的声明里给出GetRandomReal的定义，那么在类外定义GetRandomReal时， 就要加冒号，表示这个GetRandomReal()函数是类Complex的成员函数}//printf("数列长度为: %d\n", N);/*printf("Before FFT: \n");printSequence(nums, N);*/// 保存开始时间Start = timeGetTime();//printf("开始时间为: %d\n", Start);// 分配device内存，拷贝数据到devicecudaMalloc((void**)&dNums, sizeof(Complex) * N);cudaMalloc((void**)&dResult, sizeof(Complex) * N);cudaMemcpy(dNums, nums, sizeof(Complex) * N, cudaMemcpyHostToDevice);//cudaMemcpy用于在主机（Host）和设备（Device）之间往返的传递数据//主机到设备：cudaMemcpy(d_A,h_A,nBytes,cudaMemcpyHostToDevice)// 调用kernel,在GPU进行的函数通常称为核函数,一般通过__global__修饰(在核函数里，都用双下划线来修饰)int threadPerBlock = TPB;//dim3 threadPerBlock = dim3(TPB);//threadPerBlock代表线程块内含有的线程数目thread//printf("%d\n%d\n%d\n",threadPerBlock.x, threadPerBlock.y, threadPerBlock.z);/*printf("%d\n", threadPerBlock);*/ int blockNum = (N + threadPerBlock - 1) / threadPerBlock;//dim3 blockNum = dim3((N + threadPerBlock.x - 1) / threadPerBlock.x); //blockNum代表block线程块数目//printf("%d\n", blockNum); //printf("%d\n%d\n%d\n",blockNum.x, blockNum.y, blockNum.z);FFT <<< blockNum, threadPerBlock >>>(dNums, dResult, N, bits);cudaError_t err;err = cudaGetLastError(); // `cudaGetLastError` 会捕获上面代码中的最近的一个错误if (err != cudaSuccess) {printf("Error: %s\n", cudaGetErrorString(err));}// 拷贝回结果cudaMemcpy(nums, dResult, sizeof(Complex) * N, cudaMemcpyDeviceToHost); //设备到主机：cudaMemcpy(h_A,d_A,nBytes,cudaMemcpyDeviceToHost)// 计算用时End = timeGetTime();//printf("结束时间为: %d\n", End);/*printf("After FFT: \n");printSequence(nums, N);*/printf("变换用时: %dms\n", End - Start);// 释放内存free(nums);cudaFree(dNums);cudaFree(dResult);}关于CUDA加速应用程序的教程这篇文章介绍的很详细，看一遍差不多就串起来了
CUDA C/C++ 教程一：加速应用程序
再贴一个看的时候记的一个草稿纸……
timeGetTime()计算时间精度比较高，我用了这个，最终计算4*10^6序列长度的数据，这个代码运算时间竟然要七秒多……要知道不用GPU的FFT才160ms左右（这个代码在下面会详细介绍），反正这个运算时间让我直接裂开，减没了输出内容，将配置从Debug改为Release，选用独立显卡运算试了好多都不顶用，又一行又一行看代码，也没看出哪里有毛病。
后来看到说调用GPU运算数据量太少，CPU传到GPU也要耗费时间，使用GPU要运算数据量大的才有优势，好嘛，我就加序列长度，没加到10^8就爆内存了，一看运算时间还是比不用GPU的FFT慢十倍多，我彻底萎掉了。
C语言：调用fftw库的未使用GPU的FFT代码
我是看到前辈调用了fftw的库写FFT简直简便了好多，这个库分为32位的和64位的，网上可以下载，我还是贴一下我的下载链接fftw32位和64位库百度网盘下载
链接：https://pan.baidu.com/s/1kILG2Y10KUGh7g1iIOZOqw
提取码：bpgn
我的电脑是使用64位的库，安装教程参考这篇文章FFTW、Eigen库在VisualStudio中的导入和使用
文章里有fftw的详细导入教程，在打开开发人员命令提示符中，我的是下图（因为之前装过vs2015，但不影响，都能用）
在最后一步放的位置就是下图
未调用GPU的FFT代码如下：
#define _CRT_SECURE_NO_WARNINGS 1#include "fftw3.h"#include #include #include #include #pragma comment(lib,"winmm.lib")#define N 4000000void main(){int i;fftw_complex *din, *out;DWORD start, end;fftw_plan p;din = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*N);out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*N);if(din == NULL || out == NULL){printf("Error:insufficient avaliable memory\n");}else{for(i=0; i<N; i++){din[i][0] = i+1;din[i][1] = 0;}}start = timeGetTime();p = fftw_plan_dft_1d(N, din, out, FFTW_FORWARD, FFTW_ESTIMATE);fftw_execute(p);fftw_destroy_plan(p);fftw_cleanup();end = timeGetTime();printf("所用时间为%dms", end - start);//for(i=0; i<N; i++)//{//printf("%f,%fi\n", din[i][0], din[i][1]);//}//printf("\n");//for(i=0; i<N; i++)//{//printf("%f,%fi\n", out[i][0], out[i][1]);//}if(din != NULL) fftw_free(din);if(out != NULL) fftw_free(out);getchar();}同样数据量运算时间才一百多毫秒！！！啊！！
C语言：调用cufft库的GPU加速FFT
调用库速度那么快，那使用CUDA有没有库？好家伙我一搜果然有，这个cufft就是我要找的库，开始还找不到在哪里下载，找到个官网还只有linx版本的公开可以下载，到处去问人家博主，这个文章《CUFFT使用启蒙》一下子让我意识到这个cufft库就在CUDA里啊啊啊啊！！，我C盘一搜，果然有cufft.h这个文件，果断用上，代码主要是这个文章《FFTW cuFFT的使用记录》的，我改过的代码如下：
#include "cuda_runtime.h"#include "device_launch_parameters.h"#include "cufft.h"#include #include #include #include #include #pragma comment(lib,"winmm.lib")using namespace std;//#include "Complex.cu"#include#include //#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 = 3.141592653; //void test_FFTW();int main() {const int NX = 800;const int fs = 4000000;double T = 2e-6;const int BATCH = 1;//DWORD Start, End;int bei = fs/NX;double timebei = T/NX;double aaa = 1.0/NX;double Lk = 1/T;//printf("%f\n", aaa);//printf("%.20lf",timebei);cufftHandle plan;//创建句柄cufftComplex *data;cufftComplex *data_cpu;double t[NX]; 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)); //对句柄进行配置，主要是配置句柄对应的信号长度，信号类型，在内存中的存储形式等信息。//第一个参数就是要配置的 cuFFT 句柄；//第二个参数为要进行 fft 的信号的长度；//第三个CUFFT_C2C为要执行 fft 的信号输入类型及输出类型都为复数；//第四个参数BATCH表示要执行 fft 的信号的个数，新版的已经使用cufftPlanMany()来同时完成多个信号的 fft//Start = timeGetTime();//输入数据for (int i = 0; i < NX; ++i) {t[i]=i*timebei;//printf("%.15lf\n", t[i]);data_cpu[i].x = cos(PI*8*0.25*fs/(3*T*T)*pow(t[i]-T/2.0, 3)+2*PI*0.25*fs*t[i]);data_cpu[i].y = sin(PI*8*0.25*fs/(3*T*T)*pow(t[i]-T/2.0, 3)+2*PI*0.25*fs*t[i]);}//for (int i = 0; i < NX; ++i) //{////printf("%f%f\n", data_cpu[i].x, data_cpu[i].y);//double mo = sqrt(data_cpu[i].x * data_cpu[i].x + data_cpu[i].y * data_cpu[i].y);//double time = i*timebei;//printf("%.15f%.20f\n", time, mo);//}//Start = timeGetTime();//数据传输cpu->gpuCHECK(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);//第一个参数就是配置好的 cuFFT 句柄；//第二个参数为输入信号的首地址；//第三个参数为输出信号的首地址；//第四个参数CUFFT_FORWARD表示执行的是 fft 正变换；CUFFT_INVERSE表示执行 fft 逆变换。//!!!需要注意的是，执行完逆 fft 之后，要对信号中的每个值乘以 1/N//数据传输gpu->cpuCHECK(cudaMemcpy(data_cpu, data, sizeof(cufftComplex) * NX * BATCH, cudaMemcpyDeviceToHost));CHECK(cudaDeviceSynchronize());for (int i = 0; i < NX; ++i) {//printf("%f%f\n", data_cpu[i].x, data_cpu[i].y);double mo = sqrt(data_cpu[i].x * data_cpu[i].x + data_cpu[i].y * data_cpu[i].y);//double Fre = i*F;printf("%d%.15f\n", i*bei, mo*aaa);}//End = timeGetTime();//printf("变换用时: %dms\n", End - Start);cufftDestroy(plan);//释放 GPU 资源cudaFree(data);//printf("CUFFT_FORWARD:\n");system("pause");return 0;}这个4*10^6序列长度运算时间只需要10毫秒！！！！真的如官网所说的提升十倍的速度，具体为啥人家的快……我也不知道，我只是一个可怜弱小的C代码小白
上面的代码输入的是我测试结果的非线性调频信号(NLFM）二次函数形式（对于这个信号函数的matlab代码如下），正经代码都在，只是注释掉了，按需更改即可。
输入信号NLFM的matlab代码：
%% NLFMT = 2e-6; %发射脉宽fs = 400e6; %采样频率n = T*fs; %采样点数t = linspace(0,T,n);f0 = 0.25*fs;%载频[1/16fs-1/4fs]B = 0.25*fs; %带宽[1/20fs-1/4fs]K = 8*B/(3*T^2);s = exp(1i*pi*K*(t-T/2).^3+1i*2*pi*f0*t);%% 绘图figure(1);plot(t,s);%画时域图xlabel("t/s")gnuplot安装画图，maltab编写的FFT运算结果对比
用gnuplot画出C代码的输出数据，我是看到这篇文章这么做gnuplot画图
在文章的最后，我就也去安装了个gnuplot来画图，安装包和教程网上找吧，我不贴了，用的时候把代码正确运行一遍后，就按他文章中下图运行，不过不用运行第一句tcc就行
以NLFM信号为例（别问我为什么是这个，因为我最后测试的信号是这个）
输出结果如下图
matlab输出结果是：
是不是几乎一摸一样！！啊哈哈哈哈，我成功啦~~
matlab测试信号和测试时的坑
我总共用matlab测试了三种信号，原来c代码里用的等差数列看不出来，这三个matlab信号我贴在下面：
clear%% 正弦信号% fs=1000;%采样率% f1=200;%信号频率% T=1;%时宽1s% n=round(T*fs);%采样点个数% t=linspace(0,T,n);%时域横坐标% x = 3*sin(2*pi*f1*t);%形成三频信号,注意第二个频率信号幅度为2，直流幅度为3%% LFM% T = 2e-6; %发射脉宽% fs = 400e6; %采样频率% n = T*fs; %采样点数% t = linspace(0,T,n);% % f0 = 0.25*fs;%载频[1/16fs-1/4fs]% B = 0.25*fs; %带宽[1/20fs-1/4fs]% K = B/T;%调频率% %s = cos(2*pi*f0*t+2*pi*1/2*K*t.^2)+1i*sin(2*pi*f0*t+2*pi*1/2*K*t.^2);% s = exp(1i*2*pi*f0*t+1i*2*pi*1/2*K*t.^2); %LFM信号%% NLFMT = 2e-6; %发射脉宽fs = 400e6; %采样频率n = T*fs; %采样点数t = linspace(0,T,n);f0 = 0.25*fs;%载频[1/16fs-1/4fs]B = 0.25*fs; %带宽[1/20fs-1/4fs]K = 8*B/(3*T^2);s = exp(1i*pi*K*(t-T/2).^3+1i*2*pi*f0*t);%% 绘图figure(1);plot(t,s);%画时域图xlabel("t/s")grid onX = fft(s./(n));%X = fftshift(fft(s./n)); %用fft得出离散傅里叶变换aa = abs(X);%f=linspace(0,T,n);f=linspace(0,fs,n);%频域横坐标，注意奈奎斯特采样定理，最大原信号最大频率不超过采样频率的一半figure(2)plot(f,aa);%画双侧频谱幅度图xlabel("f/Hz")ylabel("幅度")grid on现在还是那个NLFM信号，自己改前面注释就行
首先我先测试了正弦信号，用前面调用cufft库的C代码，gnuplot画出的结果图，与maltab画的图对比。
这里要说明一下，C语言代码运算结果没有像matlab一样用fftshift做了对称变换，我matlab对应代码里已经改了，开始的时候让我奇怪的很，明明没有计算错误就是峰值位置不对，原来是matlab用了fftshift……无语
还有在写测试信号代码时，我算是明白了C代码的各种数值精度表示的区别，其中一个double类型除以double类型结果竟然都是大零蛋，我就在前面数值声明里除好，后面就是相乘的形式，1要改成1.0，就结果正确了，我也不知道专业的C代码大佬对于这种情况都是怎么处理的，我命名的变量名很垃圾凑合看吧，反正代码能跑……
第二个测试信号是线性调频信号（LFM），这个信号有用虚部表示，卡了我一天不知道C该怎么表示虚部，后来反应过来用欧拉公式改成实部和虚部都是三角函数的形式
原来matlab里LFM信号是这样
s = exp(1i*2*pi*f0*t+1i*2*pi*1/2*K*t.^2); %LFM信号拆成实部和虚部的形式就是这样
s = cos(2*pi*f0*t+2*pi*1/2*K*t.^2)+1i*sin(2*pi*f0*t+2*pi*1/2*K*t.^2);在C里面就是data_cpu[i].x 和data_cpu[i].y赋值的时候分别对应实部和虚部就行，虚部不用加i，等于是实部和虚部分开运算了
……嗯……就写到这里吧，该吃晚饭了，罗嗦了好多……