本文主要是介绍快速傅里叶变换的相关定义、原理及其递归算法,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!
快速傅里叶变换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,可以这样实现:
- 将a、b高位补零到合适长度
- 对a、b分别作DFT得到Fa和Fb
- Fa、Fb对应项相乘得到Fc
- 对Fc做离散傅里叶逆变换IDFT即为结果c
很显然第3步为O(N),所以关键在第2、4步。而DFT和IDFT按照定义做,同样为O(N*N),没有任何改进。这个步骤之所以可行,就是因为存在FFT算法,可以达到O(NlogN)。
为了方便书写,令WN=exp(-i*2PI/N),其中N是下标,而WN称之为单位复根。数学习惯上应该用希腊字母Ω的小写来表示单位复根。 单位复根有以下性质:
- WN^k = WN^(k+N)。很显然,因为WN的N次方就是1,这个其实是周期性
- WN^k = - WN^(k+N/2), 当N为偶数时。这也很显然,因为WN的N/2次方就是-1,这是对称性。
- WN^k = W(N/2)^(k/2),当N、k均为偶数。这个根据定义很容易推算,右边实际就是exp(-i*2PI/(N/2))^(k/2),显然等于左边
考虑一个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的递归算法非常容易实现:
当然这个递归算法完全是纯按原理实现的,考虑到程序实现,有好几个可以优化的地方。例如红色部分,完全可以进行优化,提高代码执行效率。
这篇关于快速傅里叶变换的相关定义、原理及其递归算法的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!