CUDA并行算法系列之FFT快速卷積

卷積定義

在維基百科上，卷積定義為：

離散卷積定義為：

[ 0, 1, 2, 3]和[0, 1, 2]的卷積例子如下圖所示：

Python實現（直接卷積）

根據離散卷積的定義，用Python實現：

def conv(a, b):N = len(a)M = len(b)YN = N + M - 1y = [0.0 for i in range(YN)]for n in range(YN):for m in range(M):if 0 <= n - m and n - m < N:y[n] += a[n - m] * b[m]return y

把數組b逆序，則可以不交叉計算卷積（使用numpy的array[::-1]即可實現逆序）：

import numpy as np def conv2(a, b):N = len(a)M = len(b)YN = N + M - 1y = [0.0 for i in range(YN)]b = np.array(b)[::-1] # 逆序for n in range(YN):for m in range(M):k = n - M + m + 1;if 0 <= k and k < N:y[n] += a[k] * b[m]return y

測試

可以利用numpy.convolve來檢驗計算結果的正確性：

if __name__ == '__main__':a = [ 0, 1, 2, 3 ]b = [ 0, 1, 2 ]print(conv2(a, b))print(np.convolve(a, b))

完整代碼可以在Github上找到。

利用FFT快速卷積

時域的卷積和頻域的乘法是等價的，同時時域的乘法和頻域的卷積也是等價的。基于這個這個前提，可以把待卷積的數組進行FFT變換，在頻域做乘法，然后再進行IFFT變換即可得到卷積結果。算法流程描述如下：

設N=len(a), M = len(b), 其中a, b為待卷積的數組，將長度增加到L>= N+M-1, L=2^n, n 屬于 Z，即 L=2^（log_2^(N + M - 1) +1）。

增加a, b的長度到L，后面補零。

分別計算afft = fft(a)，bfft=fft(b)。

abfft = afft × bfft。

用IFFT計算abaft的FFT逆變換，取前（N + M - 1）個值即為卷積結果。

FFT快速卷積Python代碼如下：

def convfft(a, b):N = len(a)M = len(b)YN = N + M - 1FFT_N = 2 ** (int(np.log2(YN)) + 1)afft = np.fft.fft(a, FFT_N)bfft = np.fft.fft(b, FFT_N)abfft = afft * bffty = np.fft.ifft(abfft).real[:YN]return y

測試

對比直接卷積、FFT卷積、numpy的卷積結果：

if __name__ == '__main__':a = [ 0, 1, 2, 3 ]b = [ 0, 1, 2 ]print(conv2(a, b))print(convfft(a, b))print(np.convolve(a, b))

可以看到，3個版本的計算結果是一致的。完整代碼可以在Github上找到。

性能分析

復雜度分析

直接卷積的時間復雜度為o(MN)，即o(n^2)。
FFT的時間復雜度為o(nlogn)，FFT卷積復雜度為3次FFT+L次乘法，3o(nlogn)+o(n)=o(nlogn)，及o(nlogn)。
在實際應用中，卷積核（b）被提前計算，則只需2次FFT變換。

運行測試

分別測試3個版本在數組長度為n * 1000 + 10, n=0,1,…,9的運行時間，并繪制運行時間曲線，編寫如下測試代碼：

def time_test():import timeimport matplotlib.pyplot as pltdef run(func, a, b):n = 1start = time.clock()for j in range(n):func(a, b)end = time.clock()run_time = end - startreturn run_time / nn_list = []t1_list = []t2_list = []t3_list = []for i in range(10):count = i * 1000 + 10print(count)a = np.ones(count)b = np.ones(count)t1 = run(conv, a, b) # 直接卷積t2 = run(conv2, a, b)t3 = run(convfft, a, b) # FFT卷積n_list.append(count)t1_list.append(t1)t2_list.append(t2)t3_list.append(t3)# plotplt.plot(n_list, t1_list, label='conv')plt.plot(n_list, t2_list, label='conv2')plt.plot(n_list, t3_list, label='convfft')plt.legend()plt.title(u"convolve times")plt.ylabel(u"run times(ms/point)")plt.xlabel(u"length")plt.show()

運行得到的曲線圖如下：

從圖中可知，FFT卷積比直接卷積速度要快很多。完整代碼可以在Github上找到。

CUDA實現

直接卷積

只需要把外層循環并行化就可以在CUDA上實現卷積，代碼如下：

// 直接計算卷積 __global__ void conv_kernel(const float *ina, const float *inb, float *out, size_t len_a, size_t len_b, size_t len_out) {const int tid = blockIdx.x * blockDim.x + threadIdx.x;if (tid >= len_out){return;}float sum = 0.0f;for (int m = 0; m < len_b; ++m){int k = tid - m;if (0 <= k && k < len_a){sum += ina[k] * inb[m];}}out[tid] = sum; }

當然，可以使用共享內存和常量內存（卷積核存入常量內存）進行優化，優化的代碼請查看Github。

cuFFT卷積

使用CUDA的cuFFT可以方便的進行快速傅里葉變換，cuFFT的詳細說明可以查看NVIDIA的官方文檔。本文主要使用到一下兩個函數：

• cufftExecR2C：實數到復數的快速傅里葉變換（FFT）
• cufftExecC2R：復數到實數的快速傅里葉逆變換（IFFT）

基于cuFFT的實數到復數的快速傅里葉變換代碼如下：

void fft(float *in, Complex *out, size_t size) {cufftHandle plan;cufftPlan1d(&plan, size, CUFFT_R2C, 1);cufftExecR2C(plan, in, out);cufftDestroy(plan); }

基于cuFFT的復數到實數的快速傅里葉逆變換代碼如下：

void ifft(Complex *in, float *out, size_t size) {cufftHandle plan;cufftPlan1d(&plan, size, CUFFT_C2R, 1);cufftExecC2R(plan, in, out);cufftDestroy(plan); }

其中Complex被定義為float2，typedef float2 Complex;

有了FFT，那么基于CUDA的卷積代碼可如下編寫：

void convfft( float *ina, float *inb, float *out, size_t len_out, size_t L, size_t numThreads) {thrust::device_vector<Complex> d_a_fft(L);thrust::device_vector<Complex> d_b_fft(L);thrust::device_vector<Complex> d_c_fft(L);Complex *raw_point_a_fft = thrust::raw_pointer_cast(&d_a_fft[0]);Complex *raw_point_b_fft = thrust::raw_pointer_cast(&d_b_fft[0]);Complex *raw_point_c_fft = thrust::raw_pointer_cast(&d_c_fft[0]);fft(ina, raw_point_a_fft, L);fft(inb, raw_point_b_fft, L);// 計算 d_c_fft = d_a_fft * d_b_fft;ifft(raw_point_c_fft, out, L); }

最后只剩下乘法運算了，可以自己編寫一個復數乘法的內核，也可以使用thrush的transform。使用thrush實現復數乘法，首先定義一個復數乘法操作函數（可以參考Transformations）：

struct complex_multiplies_functor {const int N;complex_multiplies_functor(int _n) : N(_n) {}__host__ __device__ Complex operator()(const Complex &a, const Complex &b) const{Complex c;c.x = (a.x * b.x - a.y * b.y) / N;c.y = (a.x * b.y + a.y * b.x) / N;return c;} };

然后使用thrush::transform即可完成計算：

// 計算 d_c_fft = d_a_fft * d_b_fft; thrust::transform(d_a_fft.begin(), d_a_fft.end(), d_b_fft.begin(), d_c_fft.begin(), complex_multiplies_functor(L));

結語

本文首先簡要介紹了卷積運算，然后使用Python實現了卷積運行的代碼，接著討論了基于FFT的快速卷積算法，并使用Python實現了FFT卷積，接著對直接卷積和基于FFT的快速卷積算法的性能進行了分析，從實驗結果可以看出，FFT卷積相比直接卷積具有更快的運行速度。最后，基于CUDA實現了直接卷積算法，并且使用cuFFT和thrush在CUDA平臺實現了基于FFT的快速卷積算法。

本文完整代碼可在Github上下載。

參考文獻

維基百科.卷積.https://zh.wikipedia.org/zh/%E5%8D%B7%E7%A7%AF

百度文庫.利用FFT計算卷積.http://wenku.baidu.com/view/5606967101f69e3143329407.html

用Python做科學計算.FFT卷積的速度比較.http://old.sebug.net/paper/books/scipydoc/example_spectrum_fft_convolve_timeit.html

NVIDIA.cuFFT.https://developer.nvidia.com/cufft

thrust. https://github.com/thrust/thrust/tree/master/thrust

分類: CUDA&GPGPU, 算法, CUDA并行算法

總結

以上是生活随笔為你收集整理的CUDA并行算法系列之FFT快速卷积的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。