MPI程序实例：FFT算法及应用

news2024/10/11 1:24:59

一、一维串行FFT算法

二、二维串行FFT算法

三、并行FFT算法

四、应用示例

4.1、多项式相乘

4.2 循环矩阵方程组的求解

1965年，两位美国科学家J.W.Cooley和J.W.Tukey发明了一种有效计算傅氏变换的方法，被称为FFT（Fast Fourier Transform，快速傅里叶变换），该算法在众多科学与工程计算中起着至关重要的作用，是20世纪计算科学的重要贡献之一。这里我们将对FFT串行算法、并行实现及FFT的应用进行介绍。

FFT快速计算：

$y_{k}=\sum_{j=0}^{n-1}x_{j}e^{-\frac{2\pi ijk}{n}},k=0,1,\cdots,n-1.\;\;\;\;\;\;\;\;\;(1)$

其中 $i^{2}=-1$ 。记 $w(n)=e^{-\frac{2\pi i}{n}}$ ，则 $w(n)^{k}$ 是方程 $x^{n}=1$ 的根。下面的性质成立：

$(1)\;\;(w(n)^{k})^{n}=1\\ (2)\;\;w(n)^{2}=w(n/2)\\ (3)\;\;w(n)^{n/2}=-1$

记 $Y=(y_{0},y_{1},\cdots,y_{n-1})^{T},X=(x_{0},x_{1},\cdots,x_{n-1})^{T},\Omega_{kj}=w(n)^{kj}$ ，公式（1）的计算过程可以写成矩阵乘向量的形式 $Y=\Omega X$ 。因此，直接计算公式（1）需要 $O(n^{2})$ 个浮点运算。由于 $w(n)$ 具有特殊性，假设 $n=2m$ ，则公式（1）可以分为以下的计算过程

$\left\{\begin{matrix} y_{k}=\sum_{j=0}^{m-1}x_{2j}w(m)^{kj}+w(n)^{k}\sum_{j=0}^{m-1}x_{2j+1}w(m)^{kj}\\ y_{k+m}=\sum_{j=0}^{m-1}x_{2j}w(m)^{kj}-w(n)^{k}\sum_{j=0}^{m-1}x_{2j+1}w(m)^{kj}\;\;\;\;\;\;\;\;\;\;\;\;(2)\\ k=0,1,\cdots,m-1 \end{matrix}\right.$

假设 $T_{n}$ 是计算所有 $y_{k}$ 的计算量，由公式（2）有 $T_{n}=2T_{n/2}+3/2n$ ，即 $T_{n}=3/2nlog_{2}n+n$ 。在此讨论的算法中，假定数据的长度是 $n=2^{m}$ 。

一、一维串行FFT算法

公式（2）是FFT的一种计算方法基础，可以通过递推的方式来完成其计算任务。在计算过程中，需要大量的数据移动，从而使得计算效率有所降低。为有效实现FFT方法，需要对数据进行重排序，使得新的数据顺序适合于计算过程。以 $n=8$ 为例，计算过程的数据依赖关系如下图。图中可知，在进行FFT算法的执行过程中，需要对原始数据进行重排序，以利于计算。数据重新排列的规则是按位倒置（bit reverse）方式进行的，以 $n=16$ 为例，按位倒置变换如下表所示。

按位倒置变换
原始顺序	第一次	第二次	第三次	十进制
0000	0000	0000	0000	0
0001	1000	1000	1000	8
0010	0001	0100	0100	4
0011	1001	1100	1100	12
0100	0010	0001	0010	2
0101	1010	1001	1010	10
0110	0011	0101	0110	6
0111	1011	1101	1110	14
1000	0100	0010	0001	1
1001	1100	1010	1001	9
1010	0101	0110	0101	5
1011	1101	1110	1101	13
1100	0110	0011	0011	3
1101	1110	1011	1011	11
1110	0111	0111	0111	7
1111	1111	1111	1111	15

假设 $n=2^{m}$ ，所有 $x_{j}$ 按照按位倒置规则进行重新排序，记为 $y_{j}$ 。令 $D(2k)=diag(w(2k)^{j}),j=0,\cdots,k-1$ ，它是k阶对角矩阵。按照公式（2），在时间上进行大幅度减少（decimation in time，DIT）的FFT算法如下：