一、FFT介紹
傅里葉變換是數字信號處理領域一個很重要的數學變換，它用來實現將信號從時域到頻域的變換，在物理學、數論、組合數學、信號處理、概率、統計、密碼學、聲學、光學等領域有廣泛的應用。離散傅里葉變換(Discrete Fourier Transform，DFT)是連續傅里葉變換在離散系統中的表示形式，由於DFT的計算量很大，因此在很長一段時間內其應用受到了很大的限制。20世紀60年代（1965年）由Cooley和Tukey提出了快速傅里葉變換(Fast Fourier Transform，FFT)算法，它是DFT的快速算法，使得離散傅里葉變換和卷積這類難度很大的計算工作的複雜度從N2量級降到了Nlog2N量級，大大提高了DFT的運算速度，從而使DFT在實際應用中得到了廣泛的應用。
傳統上，GPU只負責圖形渲染，而大部分的處理都交給了CPU。自二十世紀九十年代開始，GPU的發展迅速。由於GPU具有強大的並行計算能力，加之其可編程能力的不斷提高，GPU也用於通用計算，為科學計算的應用提供了新的選擇。
2007年6月，NVIDIA公司推出了CUDA (Compute Unified Device Architecture)，CUDA 不需要藉助圖形學API，而是採用了類C語言進行開發。同時，CUDA採用了統一處理架構，降低了編程的難度，同時，NVIDIA GPU引入了片內共享存儲器，提高了效率。這兩項改進使CUDA架構更加適合進行GPU通用計算。
二、快速傅里葉變換（FFT）
快速傅里葉變換（英語：Fast Fourier Transform, FFT），是快速計算序列的離散傅里葉變換（DFT）或其逆變換的方法[1]。傅里葉分析將信號從原始域（通常是時間或空間）轉換到頻域的表示或者逆過來轉換。
FFT會通過把DFT矩陣分解為稀疏（大多為零）因子之積來快速計算此類變換。因此，它能夠將計算DFT的複雜度從只用DFT定義計算需要的 O（n2），降低到O（nlog n），其中n為數據大小。
快速傅里葉變換廣泛的應用於工程、科學和數學領域。這裏的基本思想在1965年才得到普及，但早在1805年就已推導出來。1994年美國數學家吉爾伯特·斯特朗把FFT描述為“我們一生中最重要的數值算法”，它還被IEEE科學與工程計算期刊列入20世紀十大算法。
三、FFT的CPU實現
一維FFT基2算法的實現
我們使用按頻率抽取的方法實現了一維FFT基2算法。算法的關鍵代碼如下：
(1) 聲明雙精度複數的結構體：
struct Complex
{
double re; //複數的實部
double im; //複數的虛部
};(2) 通過冪數獲得快速傅里葉變換的長度，並初始化：
count = 1 << power; //power為冪數，count為快速傅里葉變換的長度
a = new Complex[count]; //a為初始化的數組，用來存放時域複數的值
b = new Complex[count]; //b為變換後存放結果的數組
memcpy(a, t, sizeof(Complex)*count); //初始化，將時域數據存放在a中，t為時域數據(3) 計算旋轉因子：
w = new Complex[count / 2];
for (i = 0; i<count / 2; i++)
{
angle = -i*pi * 2 / count;
w[i].re = cos(angle);
w[i].im = sin(angle);
}(4) 採用頻率分解法進行蝶形運算：
for (k = 0; k<power; k++)
{
for (j = 0; j<1 << k; j++)
{
bfsize = 1 << (power - k);
for (i = 0; i<bfsize / 2; i++)
{
p = j*bfsize;
b[i + p] = Add(a[i + p], a[i + p + bfsize / 2]); //Add()函數實現複數的加法
b[i + p + bfsize / 2] = Mul(Sub(a[i + p], a[i + p + bfsize / 2]), w[i*(1 << k)]);//Mul()函數實現複數的乘法，Sub()函數實現複數的減法
}
}
x = a;
a = b;
b = x;
}(5) 蝶形運算全部完成後，對結果進行整序，從而得到正確的輸出順序：
for (j = 0; j<count; j++)
{
p = 0;
for (i = 0; i<power; i++)
{
if (j&(1 << i))
p += 1 << (power - i - 1);
}
f[j] = a[p];
}二維FFT基2算法的實現
通過兩次調用一維快速傅里葉變換即可實現二維快速傅里葉變換。關鍵代碼如下：
(1) 計算進行傅里葉變換的寬度、高度，以及垂直方向上、水平方向上的迭代次數：
while (w * 2 <= width) //計算進行傅立葉變換的寬度（2的整數次方）
{
w *= 2; //w為變換的寬度
wp++; //wp為垂直方向上的迭代次數
}
while (h * 2 <= height)//計算進行傅立葉變換的高度（2的整數次方）
{
h *= 2; //h為變換的高度
hp++; //hp為水平方向上的迭代次數
}(2) 兩次調用一維快速傅里葉變換：
for (j = 0; j<h; j++) //在垂直方向上進行快速傅立葉變換
{
FFT1D(&t[w*j], &f[w*j], wp);
}
for (j = 0; j<h; j++) //轉換變換結果
{
for (i = 0; i<w; i++)
{
t[j + h*i] = f[i + w*j];
}
}
for (j = 0; j<w; j++) //水平方向進行快速傅立葉變換
{
FFT1D(&t[j*h], &f[j*h], hp);
}四、FFT的GPU實現
對一維或多維信號進行離散傅里葉變換的FFT變換自身具有可“分治”實現的特點，因此能高效地在GPU平台上實現。CUDA提供了封裝好的CUFFT庫，它提供了與CPU上的FFTW庫相似的接口，能夠讓使用者輕易地挖掘GPU的強大浮點處理能力，又不用自己去實現專門的FFT內核函數。使用者通過調用CUFFT庫的API函數，即可完成FFT變換。
常見的FFT庫在功能上有很多不同。有些庫採用了基2變換，只能處理長度為2的指數的FFT，而其他的一些可以處理任意長度的FFT。CUFFT4.0提供了以下功能：
(1) 對實數或複數進行一維、二維或三維離散傅里葉變換；
(2) 可以同時並行處理一批一維離散傅里葉變換；
(3) 對任意維的離散傅里葉變換，單精度最大長度可達到6400萬，雙精度最大長度可達到12800萬（實際長度受限於GPU存儲器的可用大小）；
(4) 對實數或複數進行的FFT，結果輸出位置可以和輸入位置相同（原地變換），也可以不同；
(5) 可在兼容硬件（GT200以及之後的GPU）上運行雙精度的變換；
(6)支持流執行：數據傳輸過程中可以同時執行計算過程。
一維FFT算法的CUDA實現
關鍵代碼如下：
#include <cufft.h> //CUFFT文件頭
#define NX 1024
#define BATCH 1
cufftDoubleComplex *data; //顯存數據指針
//在顯存中分配空間
cudaMalloc((void**)&data, sizeof(cufftDoubleComplex)*NX*BATCH);
//創建CUFFT句柄
cufftHandle plan;
cufftPlan1d(&plan, NX, CUFFT_Z2Z, BATCH);
//執行CUFFT
cufftExecZ2Z(plan, data, data, CUFFT_FORWARD); //快速傅里葉正變換
//銷燬句柄，並釋放空間
cufftDestroy(plan);
cudaFree(data);二維FFT算法的CUDA實現
二維FFT算法實現中，同一維FFT不同的是：
(1) 輸入參數：沒有BATCH，增加了NY。NX為行數，NY為列數；
(2) FFT的正變換的輸入位置和輸出位置不同。
關鍵代碼如下：
#include <cufft.h> //CUFFT文件頭
#define NX 1024
#define NY 1024
cufftDoubleComplex *idata, *odata; //顯存數據指針
//在顯存中分配空間
cudaMalloc((void**)&idata, sizeof(cufftDoubleComplex)*NX*NY);
cudaMalloc((void**)&odata, sizeof(cufftDoubleComplex)*NX*NY);
//創建CUFFT句柄
cufftHandle plan;
cufftPlan2d(&plan, NX, NY, CUFFT_Z2Z);
//執行CUFFT
cufftExecZ2Z(plan, idata, odata, CUFFT_FORWARD); //快速傅里葉正變換
//銷燬句柄，並釋放空間
cufftDestroy(plan);
cudaFree(idata);
cudaFree(odata);三維FFT算法的CUDA實現
#define NX 64
#define NY 64
#define NZ 128
cufftHandle plan;
cufftComplex *data1, *data2;
cudaMalloc((void**)&data1, sizeof(cufftComplex)*NX*NY*NZ);
cudaMalloc((void**)&data2, sizeof(cufftComplex)*NX*NY*NZ);
/* Create a 3D FFT plan. */
cufftPlan3d(&plan, NX, NY, NZ, CUFFT_C2C);
/* Transform the first signal in place. */
cufftExecC2C(plan, data1, data1, CUFFT_FORWARD);
/* Transform the second signal using the same plan. */
cufftExecC2C(plan, data2, data2, CUFFT_FORWARD);
/* Destroy the cuFFT plan. */
cufftDestroy(plan);
cudaFree(data1); cudaFree(data2);五、實例
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
// Include CUDA runtime and CUFFT
#include <cuda_runtime.h>
#include <cufft.h>
// Helper functions for CUDA
#include "device_launch_parameters.h"
#define pi 3.1415926535
#define LENGTH 100 //signal sampling points
int main()
{
// data gen
float Data[LENGTH] = { 1,2,3,4 };
float fs = 1000000.000;//sampling frequency
float f0 = 200000.00;// signal frequency
for (int i = 0; i < LENGTH; i++)
{
Data[i] = 1.35*cos(2 * pi*f0*i / fs);//signal gen,
}
cufftComplex *CompData = (cufftComplex*)malloc(LENGTH * sizeof(cufftComplex));//allocate memory for the data in host
int i;
for (i = 0; i < LENGTH; i++)
{
CompData[i].x = Data[i];
CompData[i].y = 0;
}
cufftComplex *d_fftData;
cudaMalloc((void**)&d_fftData, LENGTH * sizeof(cufftComplex));// allocate memory for the data in device
cudaMemcpy(d_fftData, CompData, LENGTH * sizeof(cufftComplex), cudaMemcpyHostToDevice);// copy data from host to device
cufftHandle plan;// cuda library function handle
cufftPlan1d(&plan, LENGTH, CUFFT_C2C, 1);//declaration
cufftExecC2C(plan, (cufftComplex*)d_fftData, (cufftComplex*)d_fftData, CUFFT_FORWARD);//execute
cudaDeviceSynchronize();//wait to be done
cudaMemcpy(CompData, d_fftData, LENGTH * sizeof(cufftComplex), cudaMemcpyDeviceToHost);// copy the result from device to host
for (i = 0; i < LENGTH / 2; i++)
{
printf("i=%d\tf= %6.1fHz\tRealAmp=%3.1f\t", i, fs*i / LENGTH, CompData[i].x*2.0 / LENGTH);
printf("ImagAmp=+%3.1fi", CompData[i].y*2.0 / LENGTH);
printf("\n");
}
cufftDestroy(plan);
free(CompData);
cudaFree(d_fftData);
}
