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

卷積定義

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

離散卷積定義為：

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

Python實(shí)現(xiàn)（直接卷積）

根據(jù)離散卷積的定義，用Python實(shí)現(xiàn)：

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

把數(shù)組b逆序，則可以不交叉計(jì)算卷積（使用numpy的array[::-1]即可實(shí)現(xiàn)逆序）：

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來檢驗(yàn)計(jì)算結(jié)果的正確性：

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

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

利用FFT快速卷積

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

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

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

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

abfft = afft × bfft。

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

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的卷積結(jié)果：

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個版本的計(jì)算結(jié)果是一致的。完整代碼可以在Github上找到。

性能分析

復(fù)雜度分析

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

運(yùn)行測試

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

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

運(yùn)行得到的曲線圖如下：

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

CUDA實(shí)現(xiàn)

直接卷積

只需要把外層循環(huán)并行化就可以在CUDA上實(shí)現(xiàn)卷積，代碼如下：

// 直接計(jì)算卷積 __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; }

當(dāng)然，可以使用共享內(nèi)存和常量內(nèi)存（卷積核存入常量內(nèi)存）進(jìn)行優(yōu)化，優(yōu)化的代碼請查看Github。

cuFFT卷積

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

• cufftExecR2C：實(shí)數(shù)到復(fù)數(shù)的快速傅里葉變換（FFT）
• cufftExecC2R：復(fù)數(shù)到實(shí)數(shù)的快速傅里葉逆變換（IFFT）

基于cuFFT的實(shí)數(shù)到復(fù)數(shù)的快速傅里葉變換代碼如下：

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

基于cuFFT的復(fù)數(shù)到實(shí)數(shù)的快速傅里葉逆變換代碼如下：

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);// 計(jì)算 d_c_fft = d_a_fft * d_b_fft;ifft(raw_point_c_fft, out, L); }

最后只剩下乘法運(yùn)算了，可以自己編寫一個復(fù)數(shù)乘法的內(nèi)核，也可以使用thrush的transform。使用thrush實(shí)現(xiàn)復(fù)數(shù)乘法，首先定義一個復(fù)數(shù)乘法操作函數(shù)（可以參考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即可完成計(jì)算：

// 計(jì)算 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));

結(jié)語

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

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

參考文獻(xiàn)

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

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

用Python做科學(xué)計(jì)算.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并行算法

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。