Page 211 - 《软件学报》2020年第10期

赵玉文 等:申威 26010 众核处理器上一维 FFT 实现与优化 3187

法 [38] ;文献[39]主要是对 CAM(community atmosphere model)进行了重构和优化.

3 算法背景
离散傅里叶变换(discrete Fourier transform,简称 DFT)是傅里叶变换在时域和频域上都呈离散的形式,将信
号的时域采样变换为其离散时间傅里叶变换的频域采样.使用离散傅里叶变换可以将自然科学和工程技术中
连续而复杂的问题转化为离散且简单的运算.对于一维长度为 N 的输入序列 x = [, ,...,x x 1 x N − 1 ], 其 DFT 的计算公
0
式为
N − 1
() = ∑
Xk x ( )n ω nk ,k = 0,1,..., N − 1 (1)
N
n= 0
2π
i − nk
其中, ω 为旋转因子, ω N nk = e N ,e = ix cos xi+ sin ,x i = − 1 是虚数单位,具有对称性和周期性.离散傅里叶变换
nk
N
和其逆变换(IDFT)在一段离散信号序列时域表达 x(n)和频域表达 X(k)间建立了线性变换关系.IDFT 的计算公
式为
() = ∑
xn 1 N − 1 X ( ) k ω − nk ,n = 0,1,..., N − 1 (2)
N k = 0 N
DFT 可以表示成矩阵向量乘的形式 X = F x T , 其中,
T
N
⎡ ω 0 ω 0 ω 0 " ω 0 ⎤
N
⎢ N 0 N 1 N 2 N − 1 ⎥
⎢ ω N ω N ω N " ω N ⎥
F = ⎢ ω 0 ω 2 ω 4 " ω 2(N − 1) ⎥ (3)
N
⎢ N N N N ⎥
⎢ # # # # ⎥ #
⎣ ω N 0 ω N N − ⎢ 1 ω N 2(N − 1) ω " (N − N 1)(N − 1)⎥ ⎦
2
2
( )(k =
nk
DFT 的算法复杂度为 O(N ),通常,x(n)和 X(k)是复数,当 ω 已知时,计算 Xk 0,1,..., N − 1) 需要 N 次复
N
2
数乘法和 N(N–1)次复数加法,共需要 8N 次浮点运算.
快速傅里叶变换 FFT 是离散傅里叶变换的快速算法,其中,最经典的是 1965 年由 Cooley 和 Tukey 发表的
2
Cooley-Tukey FFT 算法 [40] ,首次将离散傅里叶变换的计算复杂度从 O(N )减少到 O(N log N),该算法主要依据旋
转因子的对称性和周期性,采用分而治之的策略递归地进行计算.
Cooley-Tukey FFT 算法.对于输入长度为 N 的序列,若 N 可以分解为 N 1 和 N 2 的乘积,即 N=N 1 ×N 2 ,则可以
把输入数据按行优先自然地映射成 N 1 ×N 2 二维矩阵的形式.根据索引映射,令 n = (, )n n = 1 2 n N⋅ 1 2 + n 2 , k =
(, )kk = k + k ⋅ N , 则式(1)可表示为
2 1 1 2 1
N 2 1 N− 1 1−
() =
Xk Xk + ( k ⋅ N ) = ∑∑ ( x n N⋅ + n )ω (nN⋅ 1 2 n+ 2 )(k 1 k N+ 2 ⋅ 1 ) (4)
1 2 1 1 2 2 N
n 2 0 n= 1 0=
nk
根据 ω 的周期性和可约性,式(4)可等价变换为
N
⎡ ⎞ ⎤⎛
N 2 1 N− 1 1− N 2 1− N 1 1− ⎢ ⎟ ⎥⎜
21 ⎥⎜
⎢
⋅
nk
⋅
X (k + 1 k ⋅ 2 N 1 ) = ∑∑ ( xn N + 2 n 2 )ω n N 11 k ω nk ω nk = ∑ ∑ ( xn N + 2 n 2 )ω n N 11 k ⋅ ω⎟ N ω ⋅ nk (5)
2 2
21
2 2
N
1
N
N
1
N
n = 2 0 n = 1 0 1 2 n = 2 0 ⎢ n = 1 0 1 ⎟ ⎥⎜ 2
⎜ ⎟ ⎢ 旋转因子 ⎥
⎝
⎦
⎠ ⎣
1 FFT
N 点

N 2 点 FFT
即输出数据按行优先映射到一个 N 2 ×N 1 的二维矩阵形式.因此,N 点一维 FFT 问题的计算可以分解成 N 1 点
和 N 2 点的小规模一维 FFT 问题来完成,具体计算步骤如图 2 所示.(1) N 2 个 N 1 点一维 FFT 计算;(2) 每个元素乘
旋转因子;(3) N 1 个 N 2 点一维 FFT 计算;(4) 二维数组转置,得到正确且有序的计算结果.
如果严格按照上述 4 个步骤来执行,总共需要对主存中数组读写 4 次,即总访存量为 8N(不统计访问旋转因
子所增加的访存量).因而在实现中为了减少访存量,通常把乘旋转因子和全局数组转置分别合并到 N 1 点 FFT