4,265 热度

快速傅里叶变换FFT

转载请注明原创

似乎,我们集训队对FFT的研究者并没有太多深入,那么,我就来做一做这个东西吧。

涉及名词:

FFT(Fast Fourier Transform): 快速傅里叶变换

DFT(Discrete Fourier Transform): 离散傅里叶变换

FT(Fourier Transform): 连续傅里叶变换

卷积公式:通过两个函数f 和g 生成第三个函数

单位复根:满足{\varpi}^n=1的复数\varpi

欧拉公式:e^{ix} = \cos{x} + i\sin{x}

欧拉公式推导如下:

e^z = 1 + z + \frac{z^2}{2!} + ...+ \frac{z^n}{n!} + ...
z = iy
e^{iy} = 1 + {iy} + \frac{{(iy)}}{2!} + \frac{(iy)^3}{3!} + \frac{(iy)^4}{4!} + ... + \frac{(iy)^{2n}}{(2n)!} + \frac{(iy)^{2n+1}}{(2n+1)!} + ...
e^{iy} = [1 - {\frac{y^2}{2!}} + {\frac{y^4}{4!}} - ... + (-1)^n \cdot {\frac{y^{2n}}{(2n)!}}]\\{\hspace{5mm}}+i[y - \frac{y^3}{3!}+\frac{y^5}{5!} - ... + (-1)^n \cdot {\frac{y^{2n+1}}{(2n+1)!}}]
e^{ix} = \cos{x} + i\sin{x}

 

相消引理:对于任意偶数n > 0,有\varpi_{n}^{n/2} = \varpi_2 = -1

\varpi_{n}^{n} = \varpi_{n}^{0} = 1

快速傅里叶解决的问题类型:1、多项式相乘(大整数相乘) 2、求解偏微分方程  3、数字信号处理(这玩意跟我扯不上关系。。)等等

对于多项式乘法来说,朴素算法是O(n^2)。其表现形式为两个多项式相加。即A(x) = \sum_{i=0}^{i=n-1}a_i x^iB(x) = \sum_{i=0}^{i=n-1}b_i x^iC(x) = A(x) + B(x)

多项式的表现方式一般为两种:

  1. 系数表示法
  2. 点值表示法

系数表示法是我们正常使用的方法,所以他的计算复杂度为O(n^2),而点值表示法为O(n \log(n))

点值表示法表示对于 A(x) 我们采用了n个不同的x值以及他对应的值形成的点对组成。即\{(x_0, y_0), (x_1, y_1), ..., (x_{n-1}, y_{n-1})\},则该多项式通过下面的形式就可以唯一的表示出来:

\left[\begin{array}{ccc} 1 & x_0 & x_0^2 & ... & x_0^{n-1}\\ 1 & x_1 & x_1^2 & ... & x_1^{n-1}\\\vdots & \vdots& \vdots& \ddots& \vdots\\ 1 & x_{n-1} & x_{n-1}^2 & ... & x_{n-1}^{n-1}\end{array}\right]\left[\begin{array}{ccc} a_0 \\ a_1\\\vdots\\ a_{n-1} \end{array}\right]=\left[\begin{array}{ccc}y_0 \\ y_1\\\vdots\\ y_{n-1}\end{array}\right]

通过点值表达式,我们可以很快地计算出C(x_k) = A(x_k) + B(x_k)C(x_k) = A(x_k) * B(x_k)。而多项式C的次数界是2n的。所以,需要把A和B都扩充到2n的次数级。

在计算的过程中,我们需要先明白两个概念:

  1. 求值 :系数=>点值,即将多项式形式转化为点值表达式。
  2. 插值 :点值=>系数,即将点值转化为多项式形式。

而通过O(n\log(n))的时间复杂度实现多项式乘法的步骤由以下几步完成。

  1. 使次数界增加一倍,使多项式A和B的次数级增加到2^n。(复杂度为O(n))
  2. 求值:通过对A和B分别应用2^n阶的FFT计算出A(x)B(x)的长度为2n的点值表示。这两个点值表示中包含了两个多项式在2^n次单位根处的值。(复杂度为O(n\log(n))
  3. 点乘:把A(x)的值和B(x)的值逐点相乘,就可以计算出多项式C(x) = A(x)B(x)的点值表示,这个表示中包含了C(x)在每个2^n次单位根处的值。(复杂度为O(n))
  4. 插值:只要对2n个点值对应用一次FFT以计算出其逆DFT,就可以构造出多项式C(x)的系数表示。(复杂度为O(n\log(n)))

DFT

现在对于多项式n

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

我们取n个值{\varpi}_n^0,{\varpi}_n^1,{\varpi}_n^2, ..., {\varpi}_n^{n-1}来作为n的点值。

那么,我们定义的y如下

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

则成向量y = (y_0, y_1, y_2, ..., y_{n-1})是系数向量a = (a_0, a_1, a_2, ..., a_{n-1})离散傅立叶变换,也写作y = {DFT}_n(a)

接下来,我们就要利用FFT的方式来进行计算了。我们假设n是2的幂。也就是求值的过程。

A(x) = a_0 + a_{1}x + a_{2}x^2 + ... + a_{n-1}x^{n-1}

A(x) = a_{0}x^0 + a_{2}x^2 + a_{4}x^4 + ... + a_{n-2}x^{n-2} \\ {\hspace{11mm}} + a_{1}x^1 + a_{3}x^3 + a_{5}x^5 + ... + a_{n-1}x^{n-1}

A(x) = a_{0}x^0 + a_{2}x^2 + a_{4}x^4 + ... + a_{n-2}x^{n-2} \\ {\hspace{9mm}} + x(a_{1}x^0 + a_{3}x^2 + a_{5}x^4 + ... + a_{n-1}x^{n-2})

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

所以我们可以求出A^{[0]}A^{[1]}在点(\varpi_n^0)^2, (\varpi_n^1)^2, ... , (\varpi_n^{n-1})^2

其中,

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

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

为了更好的理解,我们看成下面

y_0 = A^{[0]}((\varpi_n^0)^2) + (\varpi_n^0)A^{[1]}((\varpi_n^0)^2)

y_1 = A^{[0]}((\varpi_n^1)^2) + (\varpi_n^1)A^{[1]}((\varpi_n^1)^2)

y_2 = A^{[0]}((\varpi_n^2)^2) + (\varpi_n^2)A^{[1]}((\varpi_n^2)^2)

\vdots

y_{n/2 - 1} = A^{[0]}((\varpi_n^{n/2 - 1})^2) + (\varpi_n^{n/2 - 1})A^{[1]}((\varpi_n^{n/2 - 1})^2)

y_{n/2 + 0} = A^{[0]}((\varpi_n^{n/2 + 0})^2) + (\varpi_n^{n/2 + 0})A^{[1]}((\varpi_n^{n/2 + 0})^2)

y_{n/2 + 1} = A^{[0]}((\varpi_n^{n/2 + 1})^2) + (\varpi_n^{n/2 + 1})A^{[1]}((\varpi_n^{n/2 + 1})^2)

y_{n/2 + 2} = A^{[0]}((\varpi_n^{n/2 + 2})^2) + (\varpi_n^{n/2 + 2})A^{[1]}((\varpi_n^{n/2 + 2})^2)

\vdots

y_{n/2 + n/2 - 1} = A^{[0]}((\varpi_n^{n/2 + n/2 - 1})^2) + (\varpi_n^{n/2 + n/2 - 1})A^{[1]}((\varpi_n^{n/2 + n/2 - 1})^2)

我们把式子在转化一下,因为\varpi_{n}^{n/2} = \varpi_2 = -1

y_0 = A^{[0]}((\varpi_n^0)^2) + (\varpi_n^0)A^{[1]}((\varpi_n^0)^2)

y_1 = A^{[0]}((\varpi_n^1)^2) + (\varpi_n^1)A^{[1]}((\varpi_n^1)^2)

y_2 = A^{[0]}((\varpi_n^2)^2) + (\varpi_n^2)A^{[1]}((\varpi_n^2)^2)

\vdots

y_{n/2 - 1} = A^{[0]}((\varpi_n^{n/2 - 1})^2) + (\varpi_n^{n/2 - 1})A^{[1]}((\varpi_n^{n/2 - 1})^2)

y_{n/2 + 0} = A^{[0]}((\varpi_n^{0})^2) - (\varpi_n^{0})A^{[1]}((\varpi_n^{0})^2)

y_{n/2 + 1} = A^{[0]}((\varpi_n^{1})^2) - (\varpi_n^{1})A^{[1]}((\varpi_n^{1})^2)

y_{n/2 + 2} = A^{[0]}((\varpi_n^{2})^2) - (\varpi_n^{2})A^{[1]}((\varpi_n^{2})^2)

\vdots

y_{n/2+n/2 - 1} = A^{[0]}((\varpi_n^{n/2 - 1})^2) - (\varpi_n^{n/2 - 1})A^{[1]}((\varpi_n^{n/2 - 1})^2)

在转化下

y_0 = A^{[0]}(\varpi_n^0) + (\varpi_n^0)A^{[1]}(\varpi_n^0)

y_1 = A^{[0]}(\varpi_n^2) + (\varpi_n^1)A^{[1]}(\varpi_n^2)

y_2 = A^{[0]}(\varpi_n^4) + (\varpi_n^2)A^{[1]}(\varpi_n^4)

\vdots

y_{n/2 - 1} = A^{[0]}(\varpi_n^{n - 2}) + (\varpi_n^{n/2 - 1})A^{[1]}(\varpi_n^{n - 2})

y_{n/2 + 0} = A^{[0]}(\varpi_n^{0}) - (\varpi_n^{0})A^{[1]}(\varpi_n^{0})

y_{n/2 + 1} = A^{[0]}(\varpi_n^{2}) - (\varpi_n^{1})A^{[1]}(\varpi_n^{2})

y_{n/2 + 2} = A^{[0]}(\varpi_n^{4}) - (\varpi_n^{2})A^{[1]}(\varpi_n^{4})

\vdots

y_{n/2+n/2 - 1} = A^{[0]}(\varpi_n^{n - 2}) - (\varpi_n^{n/2 - 1})A^{[1]}(\varpi_n^{n - 2})

所以,就可以利用二分折半的方法,来缩小规模。这样的话,计算n个点值,就可以通过O(n\log(n))的时间算出来。

在计算的时候,我们要引入复数型。

在完成了求值了之后,就是点乘。这个也很好求的。直接将两个相乘就可以了。

然后是最后一步,插值。我们把DFT写成y = V_{n}a。矩阵如下:

\left[\begin{array}{ccc}y_0 \\ y_1\\ y_2 \\ y_3 \\ \vdots\\ y_{n-1}\end{array}\right]=\left[\begin{array}{ccc} 1 & 1 & 1 & 1 & ... & 1\\ 1 & {\varpi_n} & {\varpi_n^2}& {\varpi_n^3} & ... & {\varpi_n^{n-1}}\\ 1 &{\varpi_n^2}&{\varpi_n^4}&{\varpi_n^6}& ... &{\varpi_n^{2(n-1)}}\\1 &{\varpi_n^3}&{\varpi_n^6}&{\varpi_n^9}& ... &{\varpi_n^{3(n-1)}}\\ \vdots & \vdots & \vdots& \vdots& \ddots& \vdots\\ 1 & \varpi_n^{n-1} & \varpi_n^{2(n-1)}& \varpi_n^{3(n-1)} & ... & \varpi_n^{(n-1)(n-1)}\end{array}\right]\left[\begin{array}{ccc} a_0 \\ a_1\\ a_2 \\ a_3 \\ \vdots\\ a_{n-1} \end{array}\right]

所以,接下来,这要乘以一个逆矩阵,就可以变成V_n^{-1} y_n = a_n,公式如下(这个求的过程,我就不列出来,要写好多。。。):

\left[\begin{array}{ccc} 1 / n & 1 / n & 1 / n & ... & 1 / n \\ 1 / n & \varpi_n^{-1}/n & \varpi_n^{-2}/n & ... & \varpi_n^{-(n-1)}/n\\ 1 / n & \varpi_n^{-2}/n & \varpi_n^{-4} / n & ... & \varpi_n^{-2(n-1)}/n \\ \vdots & \vdots & \vdots & \ddots\ & \vdots \\ 1 / n & \varpi_n^{-(n-1)}/n & \varpi_n^{-2(n-1)}/n & ... & \varpi_n^{-(n-1)(n-1)}/n\end{array}\right]\left[\begin{array}{ccc} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{n-1} \end{array}\right]

所以,求得a数组等于下式:

a_j = \frac{1}{n}\sum^{n-1}_{k=0}y_{k} \varpi_n^{-kj}

我们发现与之前求y数组一样。只是多了一些常量系数。这些可以等算出后面的,在出去这些常量系数。

在编写模式,会发现,FFT的整个过程类似于归并排序,其空间复杂度很大。所以,其常量因子还是蛮大的。

在100个数据,n <= 1000时的时间对比:

1. N^2的复杂度

n^2

 

2. NlogN的复杂度

NlogN

在10个数据,n <= 10000时的时间对比:

1. N^2的复杂度

N^2

2. NlogN的复杂度

NlogN

在5个数据, n <= 100000时的时间对比:

1. N^2的复杂度

N^2

2. NlogN的复杂度

nlog

可以发现,在系数比较小的情况,FFT的求解方法还是比较慢的,在n越来越大的时候,就会体现出其优越性。当然由于需要像归并排序,开空间,结果,导致,时间没有体现出那么明显。虽然,我估计我的写法还是有待优化,才开始,我使用了Vector来作传递,发现太慢了,然后,看上交的模板,又换成了指针作为传递的东西,才使效率提升了很多。

想必应该还是有很多优化的地方。

以下为我的FFT代码,参考了其他方式,增加了非递归的方式。

由于本篇文章,写了太多latex指令,所以,为了保证文章的访问速度。做题的代码,我打算放在另一篇文章里。希望大家见谅。我会在后面补上。这篇文章写得一般,希望大神们指教不足。

 

6 thoughts on “快速傅里叶变换FFT

  1. Pingback: 【转载】快速傅里叶变换FFT | Shangke7788

Leave a Reply

Your email address will not be published. Required fields are marked *