多维快速傅里叶变换

本文章最后编辑于2018年5月21日。


以下内容摘自《算法导论》思考题30-3

我们可以将一维离散傅里叶变换推广到d维上,这时输入是一个d维的数组\(A=(a_{j_1,j_2,\dots,j_d})\),维数分别为\(n_1,n_2,\dots,n_d\),其中\(n_1n_2\dots n_d=n\)。定义d维离散傅里叶变换如下:

\(y_{k_1,k_2,\dots,k_d}=\sum_{j_1=0}^{n_1-1}\sum_{j_2=0}^{n_2-1}\dots\sum_{j_d=0}^{n_d-1}a_{j_1,j_2,\dots,j_d}\omega _{n_1}^{j_1k_1}\omega _{n_2}^{j_2k_2}\dots\omega _{n_d}^{j_dk_d}\)

其中

\(0\le k_1<n_1,0\le k_2<n_2,\dots,0\le k_d<n_d\)

 

a.证明:我们可以依次在每个维度上计算一维的DFT来计算一个d维的DFT。也就是说,首先沿着第1维计算\(n/n_1\)个独立的一维DFT。然后,把沿着第一维的DFT的结果作为输入,我们计算沿着第2维的\(n/n_2\)个独立的一维DFT。利用这个结果作为输入,我们计算沿着第3维的\(n/n_3\)个独立的一维DFT,如此下去,知道第d维。

b.证明:维度的次序并无影响,于是可以通过在d个维度的任意顺序中计算一维DFT来计算一个d为的DFT。

c.证明:如果采用计算快速傅里叶变换计算每个一维的DFT,那么计算一个d维的DFT的总时间是\(O(nlogn)\),与d无关。


练习题目

Five Dimensional Discrete Fourier Transform

2017年ICPC南宁赛区的G题

直接做多维DFT,复杂度\(O(TN^6)\),需要卡卡常数。

代码:

 

快速傅里叶变换学习笔记

本文章最后编辑于2018年5月10日。


一、写在前面

最近数字图像处理课正在学快速傅里叶变换,发现自己对此理解的还不是很到位。于是借此机会,对照着《算法导论》,对这部分内容啃一啃。

两个\(n\)次多项式相加的最直接方法所需的时间是\(O(n)\),但是相乘的最直接方法所需的时间为\(O(n^2)\)。用快速傅里叶变换(Fast Fourier Transform,FFT)可以使多项式相乘的时间复杂度降低为\(O(nlogn)\)。

需要的一些前置技能:复数、多项式、线性代数。


二、多项式

一个以\(x\)为变量的多项式定义在一个代数域\(F\)上,将函数\(A(x)\)表示为形式和:

\(A(x)=\sum_{j=0}^{n-1}a_jx^j\)

我们称\(a_0,a_1,\dots,a_{n-1}\)为如上多项式的系数,所有系数都属于域\(F\),典型的情形是复数集合\(C\)。

如果一个多项式\(A(x)\)的最高次的非零系数是\(a_k\),则称\(A(x)\)的次数是\(k\),记\(degree(A)=k\)。任何严格大于一个多项式次数的整数都是该多项式的次数界,因此,对于次数界为\(n\)的多项式,其次数可以是\(0\sim n-1\)之间的任何整数。

多项式加法

如果\(A(x)\)和\(B(x)\)是次数界为\(n\)的多项式,那么它们的和也是一个次数界为\(n\)的多项式\(C(x)\),对所有属于定义域的\(x\),都有\(C(x)=A(x)+B(x)\)。也就是说,

\(A(x)=\sum_{j=0}^{n-1}a_jx^j\)

\(B(x)=\sum_{j=0}^{n-1}b_jx^j\)

\(C(x)=\sum_{j=0}^{n-1}c_jx^j(c_j=a_j+b_j)\)

例如,如果有多项式\(A(x)=6x^3+7x^2-10x+9\)和\(B(x)=-2x^3+4x-5\),那么\(C(x)=4x^3+7x^2-6x+4\)。

多项式乘法

如果\(A(x)\)和\(B(x)\)是次数界为\(n\)的多项式,那么它们的乘积\(C(x)\)是一个次数界为\(2n-1\)的多项式\(C(x)\),对所有属于定义域的\(x\),都有\(C(x)=A(x)B(x)\)。方法类似还是用上一个例子,那么得到

\(C(x)=-12x^6-14x^5+44x^4-20x^3-75x^2+86x-45\)

形式化的式子有

\(C(x)=\sum_{j=0}^{2n-2}c_jx^j\)

其中

\(c_j=\sum_{k=0}^{j}a_{k}b_{j-k}\)

此时

\(degree(C)=degree(A)+degree(B)\)

多项式的表示

从某种意义上,多项式的系数表达与点值表达式等价的。

系数表达

对一个次数界为\(n\)的多项式\(A(x)=\sum_{j=0}^{n-1}a_jx^j\)而言,其系数表达是一个由系数组成的(列)向量\(a=(a_0,a_1,\dots,a_{n-1})\)。对于多项式乘法,系数向量\(c\)成为输入向量\(a\)和\(b\)的卷积,表示成\(c=a\otimes b\)。

点值表达

一个次数界为\(n\)的多项式\(A(x)\)的点值表达就是一个由\(n\)个点值对组成的集合

\(\{(x_0,y_0),(x_1,y_1),\dots,(x_{n-1},y_{n-1})\}\)

使得对\(k=0,1,\dots,n-1\),所有\(x_k\)各不相同,且\(y_k=A(x_k)\)。

一个多项式可以有很多不同的点值表达。如果采用的点都相同的话,用点值表达多项式做乘法只需\(O(n)\)的时间。

求值与插值

从一个多项式的系数表达转化为点值表达的过程是求值,其逆运算称为插值。

定理(插值多项式的唯一性):对于任意n个点值对组成的集合\(\{(x_0,y_0),(x_1,y_1),\dots,(x_{n-1},y_{n-1})\}\),其中所有的\(x_k\)都不同,那么存在唯一的次数界为n的多项式\(A(x)\),满足\(y_k=A(x_k)\)。

证明列出矩阵方程,然后结合范德蒙德矩阵的性质。

简单的求值和插值(拉格朗日插值)的时间复杂度都是\(O(n^2)\)的。

我们之后就要通过巧妙选取点来加速这两个过程,使其运行时间变为\(O(nlogn)\)。


三、单位复数根

\(n\)次单位复数根是满足\(\omega ^n=1\)的复数\(\omega\)。

\(n\)次单位复数根恰好有\(n\)个:

\(\omega _{n}^{0},\omega _{n}^{1},\dots,\omega _{n}^{n-1}\)

其中主\(n\)次单位复数根为

\(\omega _n=e^{2\pi i/n}=\cos(2\pi/n)+i\sin(2\pi/n)\)

其他\(n\)次单位复数根都是\(\omega _n\)的幂次。

消去引理:对于任何整数\(n\ge 0,k\ge 0,d>0\),有\(\omega _{dn}^{dk}=\omega _{n}^{k}\)

推论:对于任意偶数\(n>0\),有\(\omega _{n}^{n/2}=\omega _{2}=-1\)

折半引理:如果\(n>0\)为偶数,那么\(n\)个\(n\)次单位复数根的平方的集合就是\(n/2\)个\(n/2\)次单位复数根的集合。

求和引理:对任意整数\(n\geq 1\)和不能被\(n\)整除的非负整数\(k\),有\(\sum_{j=0}^{n-1}(\omega _n^k)^j=0\)。


四、快速傅里叶变换

DFT

现在我们希望计算次数界\(n\)的多项式

\(A(x)=\sum_{j=0}^{n-1}a_jx^j\)

在\(\omega _{n}^{k}\)处的值,记为\(y_k\)

\(y_k=A(\omega _{n}^{k})=\sum_{j=0}^{n-1}a_j\omega _{n}^{kj}\)

向量\(y=(y_0,y_1,\dots,y_{n-1})\)就是系数向量\(a=(a_0,a_1,\dots,a_{n-1})\)的离散傅里叶变换(DFT),记为\(y=DFT_n(a)\)。

FFT

快速傅里叶变换(FFT)利用复数单位根的特殊性质,可以在\(O(nlogn)\)时间内计算出\(DFT_n(a)\)。首先通篇假设\(n\)恰好是2的整数幂。

FFT利用了分治策略,采用\(A(x)\)中偶数下标的系数与奇数下标的系数,分别定义两个新的次数界为\(n/2\)的多项式\(A^{[0]}(x)\)和\(A^{[1]}(x)\):

\(A^{[0]}(x)=a_{0}+a_{2}x+a_{4}x^2+\dots+a_{n-2}x^{n/2-1}\)

\(A^{[1]}(x)=a_{1}+a_{3}x+a_{5}x^2+\dots+a_{n-1}x^{n/2-1}\)

于是有

\(A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)\)

所以,求\(A(x)\)在\(\omega _{n}^{0},\omega _{n}^{1},\dots,\omega _{n}^{n-1}\)处的值转换为求次数界为\(n/2\)的多项式\(A^{[0]}(x)\)和\(A^{[1]}(x)\)在点\((\omega _{n}^{0})^2,(\omega _{n}^{1})^2,\dots,(\omega _{n}^{n-1})^2\)的值。可以发现其实是\(n/2\)个\(n/2\)次单位复数根,且每个根恰好出现两次。

IDFT

将点值表达的多项式转换回系数表达,是相似的过程。

我们把DFT写成矩阵乘积\(y=V_{n}a\)。

其中\(V_{n}\)是一个范德蒙德矩阵,在\((k,j)\)处的元素为\(\omega _{n}^{kj}\)。

对于逆运算\(a=DFT_{n}^{-1}(y)\),我们把\(y\)乘以\(V_{n}\)的逆矩阵来处理。

定理:对\(j,k=0,1,\dots,n-1\),\(V_{n}^{-1}\)在\((j,k)\)元素为\(\omega _{n}^{-kj}/n\)。

证明\(V_{n}^{-1}V_{n}=I_n\)时用求和引理即可,注意使用条件。

所以可以推导出\(DFT_{n}^{-1}(y)\):

\(a_j=\frac{1}{n}\sum_{k=0}^{n-1}y_{k}\omega _n^{-kj}\)

可以看出只需将单位根取倒数,做一次FFT,最后将结果都除以n,就做完逆变换了。


五、代码实现

首先是手写复数类,也可以用std::complex<T>,没有对比过优劣。

递归实现

需要保证n是2的幂:

合并过程的推导:

对于\(0 \le k < n/2\)

\(y_k=A(\omega _{n}^{k})\).

\(=A^{[0]}(\omega _{n}^{2k})+\omega _{n}^{k}A^{[1]}(\omega _{n}^{2k})\).

\(=A^{[0]}(\omega _{n/2}^{k})+\omega _{n}^{k}A^{[1]}(\omega _{n/2}^{k})\).

\(=y_k^{[0]}+\omega _{n}^{k}y_k^{[1]}\).

前半段没什么问题,再来看后半段

\(y_{k+(n/2)}=A(\omega _{n}^{k+(n/2)})\).

\(=A^{[0]}(\omega _{n}^{2k+n})+\omega _{n}^{k+(n/2)}A^{[1]}(\omega _{n}^{2k+n})\).

\(=A^{[0]}(\omega _{n}^{2k})-\omega _{n}^{k}A^{[1]}(\omega _{n}^{2k})\).

\(=A^{[0]}(\omega _{n/2}^{k})-\omega _{n}^{k}A^{[1]}(\omega _{n/2}^{k})\).

\(=y_k^{[0]}-\omega _{n}^{k}y_k^{[1]}\).

迭代实现

递归实际运行起来常数很大,我们需要更高效的实现方法。

先来观察一下递归过程中输入向量的下标变化,以\(n=8\)举例,可以将这个过程自行脑补成一个完全二叉树的样子:

0 1 2 3 4 5 6 7 

0 2 4 6-1 3 5 7

0 4-2 6-1 5-3 7

0-4-2-6-1-5-3-7

如果观察二进制的话会发现对应的下标是反转二进制位得到的,比如“011”变成“110”,即下标3变成了6。

代码实现举例两种

一种是直接求出对应位置反转二进制位后的数,然后交换:

第二种是从高位模拟二进制加一,用经典的摊还分析保证复杂度还是\(O(nlogn)\):

第二种的常数应该比第一种小一些(测试+分析),不过并不是很重要。

之后我们再考虑自底向上的合并,在之前的递归版本中,有一个公用子表达式\(\omega _{n}^{k}y_k^{[1]}\)计算了两次,我们可以只计算一次乘积,存放在临时变量\(t\)里,然后从\(y_k^{[0]}\)中增加或者减去\(t\),这一系列操作称为一个蝴蝶操作

代码:

on取值1或-1,on为-1代表逆变换。实际上,预处理单位根可能会更好一点。


六、模板

多项式乘法


七、结语

啊,终于弄完快速傅里叶变换了,撒花!

实际上,关于FFT还有很多东西没有讨论到,先到这里吧。

不定期更新~