Changchun Master Li

教小学生轻松理解快速傅立叶变换 fft

2018-01-15

阅读本文您不需要掌握的知识或阅读的书籍有

  • 复变函数
  • 算法导论
  • 欧拉公式

阅读本文您需要具备

  • 全日制小学生学历及其同等学历
  • 知道复数的单位是根号负一,用 i 表示

0 概要

快速傅立叶可以在最坏时间 O(nlgn) 内计算两个多项式的乘积, 本文以

为例, 计算它们的乘积

1 最小知识

我们首先了解学习 fft 所需的最小知识 点值表示法 和 单位复根。

1.1 点值表示法

多项式 A 是由它的系数构成, 用系数来表示多项式称为系数表示法:

对于 A 来说,它是一个三次多项式,有 4 个系数,想要唯一的确定这个多项式,至少需要四个点。

所以对于有 n 个系数的多项式至少需要 n 个点才能够唯一确定。
用点 来表示多项式称为点值表示法:

1.2 单位复根

我们不加证明的给出结论, 复数域上 有 n 个解,它们被称为 单位复根。记第 i 个单位复根是

例如,当 n 等于 8 时,有 8 个单位复根:

那么问题来了,我们都知道 ,那么 是如何计算的呢?

我们可以使用小学数学的知识进行分析:



我们发现,单位复根都在复平面上以原点为圆心,半径是1的圆上。
对于单位复根 , 容易得到以下性质:

若以顺时针为正,与x轴正方向的夹角为 , 则 也是单位复根,并且它与x轴的夹角为

特别,当 n 是偶数时,

(相消引理)

2 fast fourier transform

,想要快速得到多项式 ,我们要分三步。

  1. 计算 A 和 B 在单位复根处的点值表达式;O(nlgn)
  2. 点乘得到 C 的点值表达式;O(n)
  3. 插值,将点值表达式转为系数表达式。O(nlgn)

2.1 计算点值

C 的最高次为 6 次,所以最小有6个点才能确定。

所以至少要使用 8 次单位复根,即我们刚刚已经计算好的

直接算的话,复杂度至少是O{n^2},所以我们按照每项次数的次数的奇偶性分成两组来计算

(蝴蝶操作)

单位复根代入之

由复根的性质,转化为 4 次单位复根

没错,计算规模神奇的减少了一半
同理,也在B上使用同样的方法。

开局只有一把 Python 解释器,写个程序计算一下

def fft(para=[1,2,3,4,0,0,0,0]):
    if 2 == len(para):
        [a0, a1] = para
        return [a0 + a1 * 1, a0 - a1 * 1]
        
    n = len(para)
    
    odd = para[::2]
    odd_r = compute(odd)
    
    eve = para[1::2]
    eve_r = compute(eve)
    
    res = []
    for i in range(n):
        xi = complex(math.cos(i*1.0/n * 2 * math.pi), math.sin(i*1.0/n * 2 * math.pi))
        yi = odd_r[i%(n/2)] + xi * eve_r[i%(n/2)]
        res.append(yi)
    return res

解得

A = fft([1,2,3,4,0,0,0,0])
B = fft([4,6,2,7,0,0,0,0])

In [4]: A
Out[4]: 
[(10+0j),
 (-0.41421356237309426+7.242640687119286j),
 (-2.0000000000000004-1.9999999999999996j),
 (2.4142135623730954+1.2426406871192848j),
 (-2+7.347880794884119e-16j),
 (2.4142135623730945-1.2426406871192857j),
 (-1.9999999999999991+2.0000000000000004j),
 (-0.4142135623730968-7.242640687119285j)]
In [5]: B
Out[5]: 
[(19+0j),
 (3.292893218813454+11.192388155425117j),
 (1.9999999999999991-0.9999999999999998j),
 (4.707106781186548+7.1923881554251174j),
 (-7+1.592040838891559e-15j),
 (4.707106781186546-7.1923881554251174j),
 (2.000000000000001+1.0000000000000002j),
 (3.29289321881345-11.192388155425117j)]

2.2 点乘

这个没啥好说的

C_dot = []
for a, b in zip(A, B):
    C_dot.append(a * b)
In [7]: C_dot
Out[7]: 
[(190+0j),
 (-82.42640687119285+19.213203435596448j),
 (-5.999999999999998-1.9999999999999973j),
 (2.4264068711928566+23.213203435596427j),
 (14-8.327598234202001e-15j),
 (2.4264068711928424-23.21320343559642j),
 (-6.000000000000001+2.000000000000003j),
 (-82.42640687119285-19.213203435596387j)]

可以得到

2.3 插值

将点值表达式转为系数表达式即插值。

$$



$$

即求转化矩阵 的逆矩阵

容易得证 在(j, k)处的元素为

所以插值和fft算法基本相同,直接上代码吧

def rfft():
    if 2 == len(para):
        [a0, a1] = para
        return [a0 + a1 * 1, a0 - a1 * 1]
        
    n = len(para)
    
    odd = para[::2]
    odd_r = fft(odd)
    
    eve = para[1::2]
    eve_r = fft(eve)
    
    res = []
    for i in range(n):
        xi = complex(math.cos(((-i) % n)*1.0/n * 2 * math.pi), math.sin(((-i) % n)*1.0/n * 2 * math.pi))
        yi = odd_r[(-i)%(n/2)] + xi * eve_r[(-i)%(n/2)]
        res.append(yi/n)
    return res
In [10]: rfft(C_dot)
Out[10]: 
[(4+8.118390173882291e-15j),
 (13.999999999999996+5.6847126901744435e-15j),
 (25.999999999999996+1.3762094118649165e-17j),
 (45-7.538383680317139e-15j),
 (44-1.1206293398714794e-14j),
 (29.000000000000004-3.197071506826809e-15j),
 (28.000000000000004-3.783880944511323e-15j),
 (3.552713678800501e-15+8.448827874285115e-15j)]

来手动算一下验证一下结果吧

3 reference

https://en.wikipedia.org/wiki/Discrete_Fourier_transform_%28general%29

使用支付宝打赏
使用微信打赏

若你觉得我的文章对你有帮助,欢迎点击上方按钮对我打赏

扫描二维码,分享此文章