任務(wù)介紹

時隔多年仍然逃不掉寫C的命運……因為這個任務(wù)周期不短還踩了好多坑，必須記錄一下了。
任務(wù)簡單要求就是使用C語言編寫一個GPU加速的快速傅里葉變換（FFT）
分為GPU加速的FFT代碼改寫、未使用GPU的FFT編寫、運算速度對比、運算結(jié)果測試（與matlab結(jié)果對比），只要按照我文章寫的順序做就行

環(huán)境所需相關(guān)軟件下載與安裝

Visual Studio 2010
要運行C語言代碼就要先下載Microsoft Visual Studio 編輯器，我的電腦是Win10系統(tǒng)，但為了與項目環(huán)境和大家的使用習(xí)慣保持一致，使用的是Visual Studio 2010版本，網(wǎng)上有安裝包的下載，這里提供一個我保存的中文版安裝包
鏈接：https://pan.baidu.com/s/12Eyw5Woh11gtP6hyxbpUmQ
提取碼：em9h
安裝過程還挺順利的，網(wǎng)上也有很多安裝教程，就是安裝時間挺久，我安裝在D盤，最后應(yīng)用圖標(biāo)在"D:\Visual Studio 2010\Common7\IDE\devenv.exe"，安裝在其他路徑的只要找到IDE文件夾就能找到。

CUDA7.5
然后需要下載CUDA 來調(diào)用GPU。首先保證自己的電腦是有顯卡GPU的（GPU和CUDA版本對應(yīng)關(guān)系網(wǎng)上可以找到，這里插一個寫挺好的文章鏈接CUDA和GPU版本對應(yīng)），我的筆記本是雙顯卡，獨立顯卡是NVIDIA GeForce GTX 1660Ti，最高能用CUDA11.7版本的

集成顯卡：為電腦基本顯卡，在正常情況下電腦會一直使用集成顯卡，以降低電腦運行負(fù)擔(dān)，減少散熱，提高筆記本使用壽命。
獨立顯卡：為電腦的高級顯卡，多用來運行大型游戲以及部分軟件。
要改為優(yōu)先使用獨立顯卡的話，就打開NIVDIA控制面板，在左面管理3D設(shè)置中全局設(shè)置->首選圖形處理器->下拉選擇高性能NIVDIA處理器，就設(shè)置成功了

CUDA和Visual Studio 2010是有版本兼容關(guān)系的，在安裝好Visual Studio 2010后安裝CUDA需要選擇對應(yīng)版本，網(wǎng)上找到Visual Studio 2013是安裝的CUDA7.5，但我想裝一個能兼容的最高級版，就從9.0試到7.5 （CUDA歷史版本下載鏈接），要是版本不兼容在安裝CUDA的時候會出現(xiàn)不兼容的提示，就沒辦法繼續(xù)安裝，那個圖我沒保存，但安裝的時候看到就能明白（安裝教程簡單，最好安裝在默認(rèn)路徑，網(wǎng)上可找到教程），最終確定兼容版本是CUDA7.5+Visual Studio 2010

環(huán)境配置
軟件安裝完后配置環(huán)境，這里我基本是按照這篇文章教程配置的
windows+VS2013+CUDA7.5配置
只有在這篇文章中的第3步我創(chuàng)建項目的方式不太一樣
我是打開VS2010創(chuàng)建一個空項目，在源文件中添加一個新建項，選擇C++文件，這里注意：不加后綴名生成的源文件是.cpp后綴名的，這是C++代碼文件后綴名，使用CUDA的C代碼要命名為.cu的格式，不使用CUDA的C代碼命名為.c格式
如果在這篇文章第3.f步的項類型里沒有找到CUDA C++的話，就右鍵 【工程】-【生成自定義】-勾選上CUDA 7.5就行，如下圖

最后最好拿他給的測試代碼試跑一下，能成功就算環(huán)境配置成功

C語言：不調(diào)用庫的GPU加速FFT代碼

開始網(wǎng)上找到這篇文章CUDA實現(xiàn)FFT并行計算
在我簡單調(diào)試更改后可以正常運行，很開心，我對代碼加了些注釋，這里把這個能正常運行的代碼貼上，分為兩個部分代碼。
先是需要調(diào)用的函數(shù)Complex.cu文件，我跟主文件放在一起

#define PI 3.1415926

class Complex 
{
public:
    double real;
    double imag;

    Complex() 
	{

    }

    // Wn 獲取n次單位復(fù)根中的主單位根
    __device__ static Complex W(int n) 
	{
        Complex res = Complex(cos(2.0 * PI / n), sin(2.0 * PI / n));
        return res;
    }

    // Wn^k 獲取n次單位復(fù)根中的第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;
    }
    
    // 實例化并返回一個復(fù)數(shù)（只能在Host調(diào)用）
    static Complex GetComplex(double real, double imag) 
	{
        Complex r;
        r.real = real;
        r.imag = imag;
        return r;
    }

    // 隨機返回一個復(fù)數(shù)
    static Complex GetRandomComplex() 
	{
        Complex r;
        r.real = (double)rand() / rand();
        r.imag = (double)rand() / rand();
        return r;
    }

    // 隨即返回一個實數(shù)
    static Complex GetRandomReal() 
	{
        Complex r;
        r.real = (double)rand() / rand();
        r.imag = 0;
        return r;
    }

 //   // 隨即返回一個純虛數(shù)
 //   static Complex GetRandomPureImag() 
	//{
 //       Complex r;
 //       r.real = 0;
 //       r.imag = (double)rand() / rand();
 //       return r;
 //   }

    // 構(gòu)造函數(shù)（只能在Device上調(diào)用）
    __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"  // 自定義的復(fù)數(shù)數(shù)據(jù)結(jié)構(gòu)
#include <iostream>
#include <string>
#include <stdlib.h>
#include <time.h>
#include <Windows.h>
#include<mmSystem.h>
#pragma comment(lib,"winmm.lib")
using namespace std;


// 根據(jù)數(shù)列長度n獲取二進制位數(shù)
int GetBits(int n) 
{
    int bits = 0;
    while (n >>= 1) 
	{
        bits++;
    }
    return bits;
}

// 在二進制位數(shù)為bits的前提下求數(shù)值i的二進制逆轉(zhuǎn)
__device__ int BinaryReverse(int i, int bits) 
{
    int r = 0;
    do 
	{
        r += i % 2 << --bits;
    } while (i /= 2);
    return r;
}

// 蝴蝶操作, 輸出結(jié)果直接覆蓋原存儲單元的數(shù)據(jù), factor是旋轉(zhuǎn)因子
__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();    //未定義標(biāo)識符
    }
    result[tid] = nums[BinaryReverse(tid, bits)];
}

// 打印數(shù)列
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));  // 設(shè)置當(dāng)前時刻為隨機數(shù)種子
    const int TPB = 1024;  // 每個Block的線程數(shù)，即blockDim.x，每個塊中的線程數(shù)量 <= 1024
    const int N = 4000000;  // 數(shù)列大小
    const int bits = GetBits(N);
    DWORD Start, End;

    // 隨機生成實數(shù)數(shù)列
    Complex *nums = (Complex*)malloc(sizeof(Complex) * N), *dNums, *dResult;  //申請一塊內(nèi)存，其大小是N個complex長度的總和。然后把這塊內(nèi)存的首地址強轉(zhuǎn)成complex *指針變量類型，賦給*nums
    for (int i = 0; i < N; ++i) 
	{   
		nums[i] = Complex::GetRandomReal();  //沒有在類的聲明里給出GetRandomReal的定義，那么在類外定義GetRandomReal時， 就要加冒號，表示這個GetRandomReal()函數(shù)是類Complex的成員函數(shù)
    }


    //printf("數(shù)列長度為: %d\n", N);
    /*printf("Before FFT: \n");
    printSequence(nums, N);*/
    
    // 保存開始時間
    Start = timeGetTime();
    //printf("開始時間為: %d\n", Start);

    // 分配device內(nèi)存，拷貝數(shù)據(jù)到device
    cudaMalloc((void**)&dNums, sizeof(Complex) * N);
    cudaMalloc((void**)&dResult, sizeof(Complex) * N);
    cudaMemcpy(dNums, nums, sizeof(Complex) * N, cudaMemcpyHostToDevice);//cudaMemcpy用于在主機（Host）和設(shè)備（Device）之間往返的傳遞數(shù)據(jù)
    //主機到設(shè)備：cudaMemcpy(d_A,h_A,nBytes,cudaMemcpyHostToDevice)
    

    // 調(diào)用kernel,在GPU進行的函數(shù)通常稱為核函數(shù),一般通過__global__修飾(在核函數(shù)里，都用雙下劃線來修飾)
    int threadPerBlock = TPB;
	//dim3 threadPerBlock = dim3(TPB);//threadPerBlock代表線程塊內(nèi)含有的線程數(shù)目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線程塊數(shù)目
	//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));
	}

    // 拷貝回結(jié)果
    cudaMemcpy(nums, dResult, sizeof(Complex) * N, cudaMemcpyDeviceToHost); 
    //設(shè)備到主機：cudaMemcpy(h_A,d_A,nBytes,cudaMemcpyDeviceToHost)

    // 計算用時
    End = timeGetTime();
	//printf("結(jié)束時間為: %d\n", End);
    /*printf("After FFT: \n");
    printSequence(nums, N);*/
    printf("變換用時: %dms\n", End - Start);
    

    // 釋放內(nèi)存
    free(nums);
    cudaFree(dNums);
    cudaFree(dResult);
}

關(guān)于CUDA加速應(yīng)用程序的教程這篇文章介紹的很詳細，看一遍差不多就串起來了
CUDA C/C++ 教程一：加速應(yīng)用程序
再貼一個看的時候記的一個草稿紙……

timeGetTime()計算時間精度比較高，我用了這個，最終計算4*10^6序列長度的數(shù)據(jù)，這個代碼運算時間竟然要七秒多……要知道不用GPU的FFT才160ms左右（這個代碼在下面會詳細介紹），反正這個運算時間讓我直接裂開，減沒了輸出內(nèi)容，將配置從Debug改為Release，選用獨立顯卡運算試了好多都不頂用，又一行又一行看代碼，也沒看出哪里有毛病。
后來看到說調(diào)用GPU運算數(shù)據(jù)量太少，CPU傳到GPU也要耗費時間，使用GPU要運算數(shù)據(jù)量大的才有優(yōu)勢，好嘛，我就加序列長度，沒加到10^8就爆內(nèi)存了，一看運算時間還是比不用GPU的FFT慢十倍多，我徹底萎掉了。

C語言：調(diào)用fftw庫的未使用GPU的FFT代碼

我是看到前輩調(diào)用了fftw的庫寫FFT簡直簡便了好多，這個庫分為32位的和64位的，網(wǎng)上可以下載，我還是貼一下我的下載鏈接fftw32位和64位庫百度網(wǎng)盤下載
鏈接：https://pan.baidu.com/s/1kILG2Y10KUGh7g1iIOZOqw
提取碼：bpgn
我的電腦是使用64位的庫，安裝教程參考這篇文章FFTW、Eigen庫在VisualStudio中的導(dǎo)入和使用
文章里有fftw的詳細導(dǎo)入教程，在打開開發(fā)人員命令提示符中，我的是下圖（因為之前裝過vs2015，但不影響，都能用）

在最后一步放的位置就是下圖

未調(diào)用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();
}

同樣數(shù)據(jù)量運算時間才一百多毫秒?。?！啊??！

C語言：調(diào)用cufft庫的GPU加速FFT

調(diào)用庫速度那么快，那使用CUDA有沒有庫？好家伙我一搜果然有，這個cufft就是我要找的庫，開始還找不到在哪里下載，找到個官網(wǎng)還只有l(wèi)inx版本的公開可以下載，到處去問人家博主，這個文章《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;  //創(chuàng)建句柄
	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)); //對句柄進行配置，主要是配置句柄對應(yīng)的信號長度，信號類型，在內(nèi)存中的存儲形式等信息。
	//第一個參數(shù)就是要配置的 cuFFT 句柄；
    //第二個參數(shù)為要進行 fft 的信號的長度；
    //第三個CUFFT_C2C為要執(zhí)行 fft 的信號輸入類型及輸出類型都為復(fù)數(shù)；
    //第四個參數(shù)BATCH表示要執(zhí)行 fft 的信號的個數(shù)，新版的已經(jīng)使用cufftPlanMany()來同時完成多個信號的 fft

	//Start = timeGetTime();
	//輸入數(shù)據(jù)


	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();

	//數(shù)據(jù)傳輸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);
	//第一個參數(shù)就是配置好的 cuFFT 句柄；
    //第二個參數(shù)為輸入信號的首地址；
    //第三個參數(shù)為輸出信號的首地址；
    //第四個參數(shù)CUFFT_FORWARD表示執(zhí)行的是 fft 正變換；CUFFT_INVERSE表示執(zhí)行 fft 逆變換。
    //!!!需要注意的是，執(zhí)行完逆 fft 之后，要對信號中的每個值乘以 1/N



	//數(shù)據(jù)傳輸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毫秒?。。。≌娴娜绻倬W(wǎng)所說的提升十倍的速度，具體為啥人家的快……我也不知道，我只是一個可憐弱小的C代碼小白
上面的代碼輸入的是我測試結(jié)果的非線性調(diào)頻信號(NLFM）二次函數(shù)形式（對于這個信號函數(shù)的matlab代碼如下），正經(jīng)代碼都在，只是注釋掉了，按需更改即可。

輸入信號NLFM的matlab代碼：

%% NLFM

T = 2e-6;         %發(fā)射脈寬
fs = 400e6;       %采樣頻率
n = T*fs;         %采樣點數(shù)
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運算結(jié)果對比

用gnuplot畫出C代碼的輸出數(shù)據(jù)，我是看到這篇文章這么做gnuplot畫圖
在文章的最后，我就也去安裝了個gnuplot來畫圖，安裝包和教程網(wǎng)上找吧，我不貼了，用的時候把代碼正確運行一遍后，就按他文章中下圖運行，不過不用運行第一句tcc就行
以NLFM信號為例（別問我為什么是這個，因為我最后測試的信號是這個）
輸出結(jié)果如下圖

matlab輸出結(jié)果是：

是不是幾乎一摸一樣??！啊哈哈哈哈，我成功啦~~

matlab測試信號和測試時的坑

我總共用matlab測試了三種信號，原來c代碼里用的等差數(shù)列看不出來，這三個matlab信號我貼在下面：

clear
%% 正弦信號

% fs=1000;%采樣率
% f1=200;%信號頻率
% T=1;%時寬1s
% n=round(T*fs);%采樣點個數(shù)
% t=linspace(0,T,n);%時域橫坐標(biāo)
% x = 3*sin(2*pi*f1*t);%形成三頻信號,注意第二個頻率信號幅度為2，直流幅度為3

%% LFM

% T = 2e-6;         %發(fā)射脈寬
% fs = 400e6;       %采樣頻率
% n = T*fs;         %采樣點數(shù)
% t = linspace(0,T,n);
% 
% f0 = 0.25*fs;%載頻[1/16fs-1/4fs]
% B = 0.25*fs; %帶寬[1/20fs-1/4fs]
% K = B/T;                                      %調(diào)頻率
% %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;         %發(fā)射脈寬
fs = 400e6;       %采樣頻率
n = T*fs;         %采樣點數(shù)
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);%頻域橫坐標(biāo)，注意奈奎斯特采樣定理，最大原信號最大頻率不超過采樣頻率的一半
figure(2)
plot(f,aa);%畫雙側(cè)頻譜幅度圖
xlabel("f/Hz")
ylabel("幅度")
grid on

現(xiàn)在還是那個NLFM信號，自己改前面注釋就行
首先我先測試了正弦信號，用前面調(diào)用cufft庫的C代碼，gnuplot畫出的結(jié)果圖，與maltab畫的圖對比。
這里要說明一下，C語言代碼運算結(jié)果沒有像matlab一樣用fftshift做了對稱變換，我matlab對應(yīng)代碼里已經(jīng)改了，開始的時候讓我奇怪的很，明明沒有計算錯誤就是峰值位置不對，原來是matlab用了fftshift……無語
還有在寫測試信號代碼時，我算是明白了C代碼的各種數(shù)值精度表示的區(qū)別，其中一個double類型除以double類型結(jié)果竟然都是大零蛋，我就在前面數(shù)值聲明里除好，后面就是相乘的形式，1要改成1.0，就結(jié)果正確了，我也不知道專業(yè)的C代碼大佬對于這種情況都是怎么處理的，我命名的變量名很垃圾湊合看吧，反正代碼能跑……

第二個測試信號是線性調(diào)頻信號（LFM），這個信號有用虛部表示，卡了我一天不知道C該怎么表示虛部，后來反應(yīng)過來用歐拉公式改成實部和虛部都是三角函數(shù)的形式
原來matlab里L(fēng)FM信號是這樣

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賦值的時候分別對應(yīng)實部和虛部就行，虛部不用加i，等于是實部和虛部分開運算了

……嗯……就寫到這里吧，該吃晚飯了，羅嗦了好多……文章來源地址http://www.zghlxwxcb.cn/news/detail-462739.html

到了這里，關(guān)于C語言使用CUDA中cufft函數(shù)做GPU加速FFT運算，與調(diào)用fftw函數(shù)的FFT做運算速度對比的文章就介紹完了。如果您還想了解更多內(nèi)容，請在右上角搜索TOY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！