FFTW3學習筆記3：FFTW 和 CUFFT 的使用對比

一、流程

1.使用cufftHandle創建句柄

2.使用cufftPlan1d(),cufftPlan3d(),cufftPlan3d(),cufftPlanMany()對句柄進行配置，主要是配置句柄對應的信號長度，信號類型，在內存中的存儲形式等信息。 

cufftPlan1d()：針對單個 1 維信號

cufftPlan2d()：針對單個 2 維信號

cufftPlan3d()：針對單個 3 維信號

cufftPlanMany()：針對多個信號同時進行 fft

3.使用cufftExec()函數執行 fft

4.使用cufftDestroy()函數釋放 GPU 資源

 

二、單個 1 維信號的 fft

假設要執行 fft 的信號data_dev的長度為N，並且已經傳輸到 GPU 顯存中，data_dev數據的類型為cufftComplex，可以用一下方式產生主機段的data_dev。

cufftComplex *data_Host = (cufftComplex*)malloc(NX*BATCH * sizeof(cufftComplex)); // 主機端數據頭指針 // 初始數據
for (int i = 0; i < NX; i++) { data_Host[i].x = float((rand() * rand()) % NX) / NX; data_Host[i].y = float((rand() * rand()) % NX) / NX; }

然后用cudaMemcpy()將主機端的data_host拷貝到設備端的data_dev，即可用下述方法執行 fft ：

cufftHandle plan; // 創建cuFFT句柄
cufftPlan1d(&plan, N, CUFFT_C2C, BATCH); cufftExecC2C(plan, data_dev, data_dev, CUFFT_FORWARD); // 執行 cuFFT，正變換

cufftPlan1d()：

• 第一個參數就是要配置的 cuFFT 句柄；
• 第二個參數為要進行 fft 的信號的長度；
• 第三個CUFFT_C2C為要執行 fft 的信號輸入類型及輸出類型都為復數；CUFFT_C2R表示輸入復數，輸出實數；CUFFT_R2C表示輸入實數，輸出復數；CUFFT_R2R表示輸入實數，輸出實數；
• 第四個參數BATCH表示要執行 fft 的信號的個數，新版的已經使用cufftPlanMany()來同時完成多個信號的 fft。

cufftExecC2C()：

• 第一個參數就是配置好的 cuFFT 句柄；
• 第二個參數為輸入信號的首地址；
• 第三個參數為輸出信號的首地址；
• 第四個參數CUFFT_FORWARD表示執行的是 fft 正變換；CUFFT_INVERSE表示執行 fft 逆變換。

需要注意的是，執行完逆 fft 之后，要對信號中的每個值乘以 1/N

三、代碼實現

#include <iostream>
#include <time.h>
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <cufft.h>

#define NX 3335 // 有效數據個數
#define N 5335 // 補0之后的數據長度
#define BATCH 1
#define BLOCK_SIZE 1024
using std::cout;
using std::endl;


/**
* 功能：判斷兩個 cufftComplex 數組的是否相等
* 輸入：idataA 輸入數組A的頭指針
* 輸入：idataB 輸出數組B的頭指針
* 輸入：size 數組的元素個數
* 返回：true | false
*/
bool IsEqual(cufftComplex *idataA, cufftComplex *idataB, const int size)
{
    for (int i = 0; i < size; i++)
    {
        if (abs(idataA[i].x - idataB[i].x) > 0.000001 || abs(idataA[i].y - idataB[i].y) > 0.000001)
            return false;
    }

    return true;
}



/**
* 功能：實現 cufftComplex 數組的尺度縮放，也就是乘以一個數
* 輸入：idata 輸入數組的頭指針
* 輸出：odata 輸出數組的頭指針
* 輸入：size 數組的元素個數
* 輸入：scale 縮放尺度
*/
static __global__ void cufftComplexScale(cufftComplex *idata, cufftComplex *odata, const int size, float scale)
{
    const int threadID = blockIdx.x * blockDim.x + threadIdx.x;

    if (threadID < size)
    {
        odata[threadID].x = idata[threadID].x * scale;
        odata[threadID].y = idata[threadID].y * scale;
    }
}

int main()
{
    cufftComplex *data_dev; // 設備端數據頭指針
    cufftComplex *data_Host = (cufftComplex*)malloc(NX*BATCH * sizeof(cufftComplex)); // 主機端數據頭指針
    cufftComplex *resultFFT = (cufftComplex*)malloc(N*BATCH * sizeof(cufftComplex)); // 正變換的結果
    cufftComplex *resultIFFT = (cufftComplex*)malloc(NX*BATCH * sizeof(cufftComplex)); // 先正變換后逆變換的結果

    // 初始數據
    for (int i = 0; i < NX; i++)
    {
        data_Host[i].x = float((rand() * rand()) % NX) / NX;
        data_Host[i].y = float((rand() * rand()) % NX) / NX;
    }


    dim3 dimBlock(BLOCK_SIZE); // 線程塊
    dim3 dimGrid((NX + BLOCK_SIZE - 1) / dimBlock.x); // 線程格

    cufftHandle plan; // 創建cuFFT句柄
    cufftPlan1d(&plan, N, CUFFT_C2C, BATCH);

    // 計時
    clock_t start, stop;
    double duration;
    start = clock();

    cudaMalloc((void**)&data_dev, sizeof(cufftComplex)*N*BATCH); // 開辟設備內存
    cudaMemset(data_dev, 0, sizeof(cufftComplex)*N*BATCH); // 初始為0
    cudaMemcpy(data_dev, data_Host, NX * sizeof(cufftComplex), cudaMemcpyHostToDevice); // 從主機內存拷貝到設備內存

    cufftExecC2C(plan, data_dev, data_dev, CUFFT_FORWARD); // 執行 cuFFT，正變換
    cudaMemcpy(resultFFT, data_dev, N * sizeof(cufftComplex), cudaMemcpyDeviceToHost); // 從設備內存拷貝到主機內存

    cufftExecC2C(plan, data_dev, data_dev, CUFFT_INVERSE); // 執行 cuFFT，逆變換
    cufftComplexScale << <dimGrid, dimBlock >> > (data_dev, data_dev, N, 1.0f / N); // 乘以系數
    cudaMemcpy(resultIFFT, data_dev, NX * sizeof(cufftComplex), cudaMemcpyDeviceToHost); // 從設備內存拷貝到主機內存

    stop = clock();
    duration = (double)(stop - start) * 1000 / CLOCKS_PER_SEC;
    cout << "時間為 " << duration << " ms" << endl;

    cufftDestroy(plan); // 銷毀句柄
    cudaFree(data_dev); // 釋放空間

    cout << IsEqual(data_Host, resultIFFT, NX) << endl;

    return 0;
}

 

四、用fftw和cufft實現傅里葉變換

1.創建C++的文件命名為fftw.cpp，配置fftw環境（環境配置移步：這里），復制以下代碼

#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "fftw3.h"
#include <windows.h>
#include <Eigen/Dense>
#include <iostream>  
#include <opencv2/core/eigen.hpp>
#include <opencv2/opencv.hpp>
#include <iostream>

using namespace cv;
using namespace std;
using namespace Eigen;

#define COLS 3
#define ROWS 3

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

extern "C"    void iteration_mat1();

/**********************************主函數****************************************/
int main()
{

    fftw_complex*result_temp_din, *result_temp_out;
    fftw_plan p;

    result_temp_din = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*COLS*ROWS);
    result_temp_out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*COLS*ROWS);
    cout << "fftw" << endl;

    for (size_t j = 0; j < ROWS; j++)
    {
        for (size_t i = 0; i < COLS; i++)
        {

            result_temp_din[i + j*COLS][0] = (i+1)*(j+1);
            cout << result_temp_din[i + j*COLS][0] << " ";
            result_temp_din[i + j*COLS][1] = 0;
        }
    }

    //forward fft  
    p = fftw_plan_dft_2d(ROWS, COLS, result_temp_din, result_temp_out, FFTW_FORWARD, FFTW_ESTIMATE);
    
    fftw_execute(p);
    cout << endl;
    for (size_t j = 0; j < ROWS; j++)
    {
        for (size_t i = 0; i < COLS; i++)
        {
            cout << result_temp_out[i + j*COLS][0] << " ";//實部  
            cout << result_temp_out[i + j*COLS][1] << endl;//虛部  
        }
    }

    cout << "cuda" << endl;
    iteration_mat1();
    system("pause");
    return 0;
}

 

2.創建cuda文件命名為cufft.cu，配置環境（環境配置移步：這里），復制以下代碼

注： cufftPlan2d(&p, ROWS, COLS, CUFFT_C2C); 看清楚rows和cols，千萬別出錯！

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <cufft.h>
#include <stdio.h>
#include <opencv2/opencv.hpp>
#include <iostream>

using namespace std;
using namespace cv;

#define COLS 3
#define ROWS 3

extern "C"    void iteration_mat1()
{
    cufftComplex *result_temp_din = (cufftComplex*)malloc(COLS*ROWS * sizeof(cufftComplex));
    cufftHandle p;
    //輸入賦值數據
    for (size_t j = 0; j < ROWS; j++)
    {
        for (size_t i = 0; i < COLS; i++)
        {
            result_temp_din[i + j*COLS].x = (i + 1)*(j + 1);
            cout << result_temp_din[i + j*COLS].x << " ";
            result_temp_din[i + j*COLS].y = 0;
        }
    }
    cout << endl;

    size_t pitch;

    cufftComplex *t_result_temp_din;
    cudaMallocPitch((void**)&t_result_temp_din, &pitch, COLS * sizeof(cufftComplex), ROWS);
     
    cufftComplex *t_result_temp_out;
    cudaMallocPitch((void**)&t_result_temp_out, &pitch, COLS * sizeof(cufftComplex), ROWS);

    //將值輔到Device
    //cudaMemcpy2D(t_result_temp_din, pitch, result_temp_din, COLS * sizeof(cufftComplex), COLS * sizeof(cufftComplex), ROWS, cudaMemcpyHostToDevice);
    cudaMemcpy(t_result_temp_din,result_temp_din,  ROWS * sizeof(cufftComplex)* COLS, cudaMemcpyHostToDevice);

    //forward fft  制定變換規則
    cufftPlan2d(&p, ROWS, COLS, CUFFT_C2C);

    //執行變換
    cufftExecC2C(p, (cufftComplex*)t_result_temp_din, (cufftComplex*)t_result_temp_out, CUFFT_FORWARD);

    //將值輔到host
    cudaMemcpy(result_temp_din,  t_result_temp_out, ROWS * sizeof(cufftComplex)* COLS, cudaMemcpyDeviceToHost);
    //cudaMemcpy2D(result_temp_din, pitch, t_result_temp_out, COLS * sizeof(cufftComplex), sizeof(cufftComplex)* ROWS, COLS, cudaMemcpyDeviceToHost);


    //提取實部和虛部
    for (size_t j = 0; j < ROWS; j++)
    {
        for (size_t i = 0; i < COLS; i++)
        {
            cout << result_temp_din[i + j*COLS].x << " ";//實部  
            cout << result_temp_din[i + j*COLS].y << endl;//虛部  
        }
    }

}

 

 

3.執行結果：