任务介绍
时隔多年仍然逃不掉写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.1415926
class 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 <iostream>
#include <string>
#include <stdlib.h>
#include <time.h>
#include <Windows.h>
#include<mmSystem.h>
#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，每个块中的线程数量 <= 1024
const int N = 4000000; // 数列大小
const int bits = GetBits(N);
DWORD Start, End;
// 随机生成实数数列
Complex *nums = (Complex*)malloc(sizeof(Complex) * N), *dNums, *dResult; //申请一块内存，其大小是N个complex长度的总和。然后把这块内存的首地址强转成complex *指针变量类型，赋给*nums
for (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内存，拷贝数据到device
cudaMalloc((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 <stdio.h>
#include <stdlib.h>
#include <windows.h>
#include <MMSystem.h>
#pragma comment(lib,"winmm.lib")
#define N 4000000
void 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 <string>
#include <stdlib.h>
#include <time.h>
#include <Windows.h>
#include <mmSystem.h>
#pragma comment(lib,"winmm.lib")
using namespace std;
//#include "Complex.cu"
#include<iostream>
#include <math.h>
//#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->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);
//第一个参数就是配置好的 cuFFT 句柄；
//第二个参数为输入信号的首地址；
//第三个参数为输出信号的首地址；
//第四个参数CUFFT_FORWARD表示执行的是 fft 正变换；CUFFT_INVERSE表示执行 fft 逆变换。
//!!!需要注意的是，执行完逆 fft 之后，要对信号中的每个值乘以 1/N
//数据传输gpu->cpu
CHECK(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代码：
%% NLFM
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 = 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信号
%% NLFM
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 = 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 on
X = 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，等于是实部和虚部分开运算了文章来源：https://www.toymoban.com/news/detail-462739.html
……嗯……就写到这里吧，该吃晚饭了，罗嗦了好多……文章来源地址https://www.toymoban.com/news/detail-462739.html
到了这里，关于C语言使用CUDA中cufft函数做GPU加速FFT运算，与调用fftw函数的FFT做运算速度对比的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！
