快速傅里叶变换的相关定义、原理及其递归算法

时间:2022-12-25 17:08:04

    快速傅里叶变换FFT是离散傅里叶变换DFT的一种快速算法,实际上诸如Matlab等科学计算软件都已经实现了FFT,只需调用相应的接口即可。在ACM里,FFT的典型应用就是大数的乘法或者多项式的乘法。顺便,如果题目规模不是很大,有关大数的运算推荐使用Java语言,使用java.math.BigInteger包完成;包括高精度运算,可以使用BigDecimal包完成。任何情况下,会一门外语总是很重要的!

    傅里叶变换一般指的是连续傅里叶变换,其定义如下:


快速傅里叶变换的相关定义、原理及其递归算法快速傅里叶变换的相关定义、原理及其递归算法

    傅里叶逆变换的定义如下:

快速傅里叶变换的相关定义、原理及其递归算法


    离散傅里叶变换与连续傅里叶变换“类似”,只不过把积分换成级数。DFT以及IDFT的定义式如下:


快速傅里叶变换的相关定义、原理及其递归算法


快速傅里叶变换的相关定义、原理及其递归算法


    令源序列x={0, 1, 3},则DFT的结果是X={4, exp(-i*2PI/3)+3exp(-i*4PI/4), exp(0i*4PI/3)+3exp(-i*8PI/3)}  ;

另外一个例子, 若令源序列为x={2, 4},则DFT的结果是X={6, -2}。

  

    FFT和卷积关系密切,连续卷积和离散卷积的定义分别如下:

快速傅里叶变换的相关定义、原理及其递归算法

快速傅里叶变换的相关定义、原理及其递归算法

    下面举一个简单的例子显示离散卷积的运算,令序列f={0,1,3},序列g={2,4},则其卷积的计算式如下:

快速傅里叶变换的相关定义、原理及其递归算法

所以,卷积的结果序列y={0, 2, 10, 12}


    卷积可以使用DFT完成,本质上这也是大数乘法、多项式乘法可以使用DFT完成的原理。当然在完成时,首先要将源序列的长度补齐到合适的数值。使用DFT完成卷积或者多项式乘法时,序列长度应该是2的幂,显然又不能小于最终结果的长度。上例中结果长度为4,恰好是2的幂,所以源序列补到4长度即可。于是 f = {0, 1, 3, 0},g = {2, 4, 0, 0}。NOTE:对大数乘法一定要注意结果长度的估算,本例如果表示大数乘法,则是一个3位数乘2位数,其结果最多是一个5位数,所以长度应补到8为佳。

    f和g的DFT分别是

快速傅里叶变换的相关定义、原理及其递归算法快速傅里叶变换的相关定义、原理及其递归算法

     将其对应项相乘,得到

快速傅里叶变换的相关定义、原理及其递归算法

    再对该序列做一个离散傅里叶逆变换,可以得到{0, 2, 10, 12},与原始卷积的结果相同,在本质上这也是两个多项式乘法的结果。另一方面,如果考虑三位数310和两位数42的乘法,其结果是13020,与卷积的结果同样是吻合的。注意乘法运算需要考虑进位!

    所以,给定多项式a、b,求其乘积c,可以这样实现:

  1. 将a、b高位补零到合适长度
  2. 对a、b分别作DFT得到Fa和Fb
  3. Fa、Fb对应项相乘得到Fc
  4. 对Fc做离散傅里叶逆变换IDFT即为结果c

   很显然第3步为O(N),所以关键在第2、4步。而DFT和IDFT按照定义做,同样为O(N*N),没有任何改进。这个步骤之所以可行,就是因为存在FFT算法,可以达到O(NlogN)。

    为了方便书写,令WN=exp(-i*2PI/N),其中N是下标,而WN称之为单位复根。数学习惯上应该用希腊字母Ω的小写来表示单位复根。 单位复根有以下性质:

  1. WN^k = WN^(k+N)。很显然,因为WN的N次方就是1,这个其实是周期性
  2. WN^k = - WN^(k+N/2), 当N为偶数时。这也很显然,因为WN的N/2次方就是-1,这是对称性。
  3. WN^k = W(N/2)^(k/2),当N、k均为偶数。这个根据定义很容易推算,右边实际就是exp(-i*2PI/(N/2))^(k/2),显然等于左边
    此处仅以一个8点DFT来描述FFT的原理,如果想看到更一般性、更本质的描述及证明,可以去阅读《Introduction to Algorithms》。

    考虑一个8点DFT,其原始公式是:

快速傅里叶变换的相关定义、原理及其递归算法

将其奇偶项分开书写,可写成:

快速傅里叶变换的相关定义、原理及其递归算法

将奇数项提取出一个公因子,可写成(注意下标n的变化):

快速傅里叶变换的相关定义、原理及其递归算法

注意到单位复根的脚标是8,为偶数,根据单位复根的性质3,可以做一个比例变换,得到:


快速傅里叶变换的相关定义、原理及其递归算法


于是再对求和符号的下标n做一个变换,得到:


快速傅里叶变换的相关定义、原理及其递归算法

接下来,按照k值,将上式写成2个部分:

快速傅里叶变换的相关定义、原理及其递归算法

    注意上式中的第1个式子,该式的前后两个部分都是一个4点序列的DFT,所以该式实际上就是2个4点DFT之和。至于第2个式子,首先变化k将其改写为:

快速傅里叶变换的相关定义、原理及其递归算法

然后考虑到单位复根的性质1和性质2,上式可以改写为:

快速傅里叶变换的相关定义、原理及其递归算法

    所以最终一个8点DFT变为了:

快速傅里叶变换的相关定义、原理及其递归算法

 

    该式表明一个8点DFT转为了2个4点DFT,8点DFT的结果中的前4项是这两个4点DFT之和(带一个单位复根的因子),而后4项是其差。递归下去,4点DFT可以变为2个2点DFT,而2点DFT变为2个1点DFT,最后1点DFT就是其本身。

    所以,基于此原理,FFT的递归算法非常容易实现:

快速傅里叶变换的相关定义、原理及其递归算法


    当然这个递归算法完全是纯按原理实现的,考虑到程序实现,有好几个可以优化的地方。例如红色部分,完全可以进行优化,提高代码执行效率。