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并行算法