CUDA并行算法系列之FFT快速卷积
卷积定义
在维基百科上,卷积定义为:
离散卷积定义为:
[ 0, 1, 2, 3]
和[0, 1, 2]
的卷积例子如下图所示:
Python实现(直接卷积)
根据离散卷积的定义,用Python实现:
def conv(a, b):
N = len(a)
M = len(b)
YN = N + M - 1
y = [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 - 1
y = [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=2n,n∈Z ,即