0%

傅立叶级数,变换以及快速算法

按照惯例,开坑先写序。

写此文的主要原因还是因为以前花了钞票买的SDR吃了灰,(一直想着把FFT弄明白之后再好好的玩一下SDR)。
果然花钱的好处(尤其是花很多钱),总是在某个深夜冒出的愧疚之心能够鞭策你稍微写点什么东西。
前段时间,没事的时候又看了一下FFT代码实现。之前也顶多就是顺手拿来用,没有很好的去理解其中的原理。而后就又去看了很多关于的FFT的推导文,以及3Blue1Brown的视频,(傅立叶变换,初等群论简介)。
数学最美妙的地方莫过于推导和应用。前者让你仿佛感觉自己稍稍的触及了宇宙冰山的一角,而后者更让你感受到了大学的学费没有白花。
所以写此文章做一个推导流程的复盘和记录,顺带理一下其中一些混淆不清的东西。
顺便增加一下Latex的熟练度。

FFT-源起多项式


FFT(FAST FOURIER TRANSFORM),全称是快速傅立叶变换,其包括傅立叶变换对的两个部分:傅立叶变换和傅立叶逆变换。
看大多数的文章在讲FFT的时候,通篇下来大多数都和傅立叶变换或者傅立叶级数八杆子都打不着一块。因为FFT这个东西,如果不看名字的话,这货其实是用来优化多项式乘法问题的。
而实际上,傅立叶级数其实又可以看作是多项式的一种,DTFT(离散傅立叶变换)也可以使用FFT进行加速。
既然万物起源多项式,那就先从多项式开始说起。

多项式表达方式


一个n次的n+1项多项式的两种表达方式:一种是系数表达式,一种是点值表达式

  • 系数表达式
    有如下多项式:
    f(x)=a0+a1x+a2x2+a3x3+...+an+1xn;f(x) = a_0 + a_1 \cdot x + a_2 \cdot x^2 + a_3 \cdot x^3 + ... + a_{n+1} \cdot x^n;
    通常情况下,从自打学习到函数开始,我们看到的多项式基本都是这样。
    在系数表达式中,最重要的就是多项式的系数。
    在抽象向量空间,可以用一个向量来描述这个多项式:
    f(x)=[a0a1a2a3an+1]f(x) = \begin{bmatrix} a_0 & a_1 & a_2 & a_3 & \cdots & a_{n+1} \end{bmatrix}
  • 点值表示法
    假设我们在选取n+1个{x,y}点,那么,可以通过这n+1个点,来唯一表达出一个最高次数为n次项的多项式。
    为森么呢?
    原因:因为对于多项式来说,确定的系数,决定了唯一的多项式。
    我们把点值带入多项式中,来求出这个多项式系数矩阵,相当于是求一个n+1元1次方程。
    当带入点数为n+1次的时候,这个方程有唯一解。
    用矩阵来表示一下这个求解方程
    [1x0x02x03x0n1x1x12x13x1n1x2x22x23x2n1x3x32x33x3n1xnxn2xn3xnn][a0a1a2a3an]=[f(x0)f(x1)f(x2)f(x3)f(xn)]\begin{bmatrix}1 & x_0 & x_0^2 & x_0^3 & \cdots & x_0^n\\1 & x_1 & x_1^2 & x_1^3 & \cdots & x_1^n\\1 & x_2 & x_2^2 & x_2^3 & \cdots & x_2^n\\1 & x_3 & x_3^2 & x_3^3 & \cdots & x_3^n\\\vdots & \cdots & \cdots& \cdots & \vdots\\1 & x_n & x_n^2 & x_n^3 & \cdots & x_n^n\\\end{bmatrix}\begin{bmatrix} a_0\\a_1\\a_2\\a_3\\\vdots\\a_n\end{bmatrix}=\begin{bmatrix} f(x_0)\\f(x_1)\\f(x_2)\\f(x_3)\\\vdots\\f(x_n)\end{bmatrix}
    因此,由于n+1个固定点能表达一个最高次数为n次项的多项式。
    我们可以记做:
    f(x)={{x0,f(x0)},{x1,f(x1)},{x2,f(x2)},,{xn,f(xn)}}f(x) = \{\{x_0,f(x_0)\}, \{x_1,f(x_1)\}, \{x_2,f(x_2)\}, \cdots, \{x_n,f(x_n)\}\}

以上两种表达方式,是可以互相转换的。
如果要把系数表达式转换成点值表达式,直接带入n+1个点就行了。
如果要把点值表达式转换成系数表达式,直接根据上面那个矩阵,求解n+1元一次方程即可。

粗粗一看,这个转换方式和傅立叶变换对冥冥之中就是有所联系。

不同表达式下多项式乘法的区别


上一小节在简单说了多项式的两种表达方式,之前有提到FFT这个东西的根本目的是用来优化多项式乘法的运算速度(引申,也就是说只要是能抽象成多项式乘法的问题都能用FFT来加速)。
因此以下要说的就是两种不同表达方式下多项式乘法的区别。

系数表示法下的多项式乘法


想象两个n次的多项式相乘,生成一个新的多项式。
在系数表达式的情况下,其复杂度为O(n2)O(n^2),此处的复杂度很好理解,因为一个多项式的每一项要和另外一个多项式的每一项进行相乘。

ps:吐槽一下,之前一直没搞清楚这个复杂度为啥是这个,其实是搞混了多项式和多项式结果的区别,因此,下面直接用矩阵把这一运算过程表示出来。

假定一个多项式的系数矩阵是 A=[a0an]A=\begin{bmatrix}a_0 & \cdots & a_n \end{bmatrix}B=[b0bn]B=\begin{bmatrix}b_0 & \cdots & b_n \end{bmatrix} 生成的新的多项式 C=[c0c2n]C=\begin{bmatrix}c_0 & \cdots & c_{2n} \end{bmatrix}

  • 两个n次多项式相乘生成最高2n次多项式这一过程可以表达成:
    [a00a1a00a2a1a00anan1an2a00anan1a100ana20an][b0b1b2b3b4b5b6bn]=[c0c1c2c3c4c5c6c2n]\begin{bmatrix}a_0 & \cdots & \cdots & \cdots & 0\\a_1 & a_0 & \cdots & \cdots & 0\\a_2 & a_1 & a_0 & \cdots & 0\\\vdots & \ddots & \ddots & \ddots & \vdots\\a_n & a_{n-1} & a_{n-2} & \ddots & a_0\\0 & a_{n} & a_{n-1} & \ddots & a_1\\0 & 0 & a_{n} & \ddots & a_2\\\vdots & \ddots & \ddots & \ddots & \vdots\\0 & \cdots & \cdots & \cdots & a_n\end{bmatrix}\begin{bmatrix} b_0\\ b_1\\b_2\\b_3\\b_4\\b_5\\b_6\\\vdots\\b_n\end{bmatrix}=\begin{bmatrix} c_0\\c_1\\c_2\\c_3\\c_4\\c_5\\c_6\\\vdots\\c_{2n}\end{bmatrix}
仔细看左边这个矩阵,是不是很像卷积操作,数学的神奇和迷人之处就在于万物之间仿佛都存在那摸不着看不见但却又切实存在的千丝万缕的联系

点值表示法下的多项式乘法


再来看点值表示法下的多项式乘法,简单到怀疑人生。

  • 假设两个最高n次的多项式,点值表示法情况分别如下,带入到n+1个x均相同:
    f(x)={{x0,f(x0)},{x1,f(x1)},{x2,f(x2)},,{xn,f(xn)}}f(x) = \{\{x_0,f(x_0)\}, \{x_1,f(x_1)\}, \{x_2,f(x_2)\}, \cdots, \{x_n,f(x_n)\}\}
    g(x)={{x0,g(x0)},{x1,g(x1)},{x2,g(x2)},,{xn,g(xn)}}g(x) = \{\{x_0,g(x_0)\}, \{x_1,g(x_1)\}, \{x_2,g(x_2)\}, \cdots, \{x_n,g(x_n)\}\}
  • 那么这个多项式相乘的算计结果(记做Q(x)Q(x),就是直接把对应的点值相乘即可。
    Q(x)={{x0,f(x0)g(x0)},{x1,f(x1)g(x1)},,{xn,f(xn)g(xn)}}Q(x) = \{\{x_0,f(x_0)*g(x_0)\}, \{x_1,f(x_1)*g(x_1)\}, \cdots, \{x_n,f(x_n)*g(x_n)\}\}

理解起来很简单,如果两个函数f(x)g(x)f(x)*g(x)相乘,那么其再x0x_0处的值,肯定是等于f(x0)g(x0)f(x_0)*g(x_0)
那么点值表达式下的多项式乘法的复杂度,就是O(n)O(n)

想一下,信号系统里,时域下的卷积等于频域下的乘积,是否意味着同一信号的时域和频率,其实就是这个信号的系数表示法和点值表示法?点又代表里什么?采样?

结论


既然复杂度差距那么大,这里面是否可以做文章呢?
初步一想,直接那点值表示法的多项式做乘法,得到的新的点值表达式,其值表示的也就是离散的点。
如果说想要任意x下的多项式值,还是需要使用系数表达式。
所以说我们可以把多项式转换成点值表达式,然后完成多项式相乘运算之后,再转换回系数表示法。

可是在这种情况下,复杂度会在转换步骤,重新升回O(n2)O(n^2)
仿佛就像是物理学中的能量守恒,令人沮丧。

然而FFT的出现,降低里变换对的复杂度,从而使整个过程的复杂度,降低到O(nlog2n)O(nlog_2 n)

嗯,之前的熵还不够大

复数,复数域


在做FFT的原理说明之前,需要对复数和复数域做一个简单的说明。
如果说,多项式的系数和点值表达式之间的变换动作,是取值代入计算和解方程组。
那么如果说要优化这一步,那么实数内的点已经没有办法做到里,这里要找的点,存在于复数域

欧拉公式和初等群论


复数以前就是个让人很头痛的东西。
以前读书的时候,复数这个概念大概最早也是从求解x2=1x^2 = -1的虚数根引入的。
但是实际上,一开始接受虚数概念的时候,完全没有办法理解x2=1x=ix^2 = -1 x=i的原因(或者说不能理解为什么i2=1i^2=-1)
按课本上入门方式,实在是太陡峭,根本没办法很好的去了解这东西到底是个啥,有什么用。

此处通过简单引入初等群论中的部分概念,来解释为什么x2=1x^2 = -1以及进一步的了解复数的应用。

对,借助复数乘法群,你可以很自然的理解x2=1x^2 = -1,就像理解1+1=21+1=2一样自然。

初等群论


学识有限,以下为简单不严谨的概括一些群论是干啥的:

  • 就是某个数学对象对象上对称作用的集合
    • 数学对象,可以理解成为任意的一个集合。
    • 对称作用相对难理解,是指对此集合内对象操作之后,依然得到是这个集合内对象。对于集合图形相对好理解,比如正方形绕中轴旋转180度依然是正方形。
  • 群论主要是研究群的代数结构(简单理解就是某种作用等效于那些作用的相继进行)
    群论的好处就是,在研究对称作用的一些性质后,可以进行更广义抽象意义上的推导。

比如说:

  • 复数加法群:其对象就是整个复域,其对称作用就是复数加法
  • 复数乘法群:其对象就是整个复域,其对称作用就是复数乘法

对于复数域中的:1+i1+i11这两个数字

  • 如果用复数加法群中的对称作用,可以了解到,1+i1+i可以通过在11上作用一个+i+i的移动来得到。
  • 如果用复数乘法群中的对称作用,可以了解到,1+i1+i可以通过在11上先拉伸2\sqrt{2}然后旋转45度角来实现(也可以从极坐标角度理解)

因此同一种结果的映射,在不同的群中可以由不同的变换来构成。

  • 不同群中相同效果的作用叫做同态
举个不恰当的例子:对于早期不支持fpu浮点运算的处理器,其处理实数的乘法操作也可以通过重复累加操作实现,其减法操作,也是通过和补码的加法操作实现,也就是说这些对称操作实际是通过加法操作的某些相继作用来等效实现的。

复数域的乘法计算特点和欧拉公式


按照群的映射,对于一个复数加法群对称操作集合,可以通过f(x)=exf(x)=e^x映射到复数乘法群的操作集合:

  • 对于在复数轴上的加法操作,可以映射为旋转变换x每变化i,旋转角转换1弧度(因为幂的底数为e)
  • 对于在实数轴上的加法操作,可以映射为拉伸变换x每变换1,模长拉伸e倍

后面解释的时候要用到这个操作的映射关系。

从另外一个角度理解:

  1. 根据三角函数泰勒展开证明,可以得到cosθ+isinθ{cos\theta\\+ isin\theta\\}的另一等价表示方法eiθ{e^{i\theta\\}}
  2. 对于cosθ+isinθ{cos\theta\\+isin\theta\\},按照复数乘法群中的对称作用,可以知道其可以从11旋转θ\theta\\度而来。

使用eiθe^{i\theta\\}来进行复数域的数值表示而非a+ia+i的好处就是能简化复数域乘法。
因为对于指数表达式来说复数乘法就是指数相加,模值相乘,也就是复数乘法群来说分别就是旋转拉伸变换的数值表达。
根据这一性质,通过复数乘法可以很好的表达旋转和拉伸动作。

回到i2=1i^2=-1上,通过旋转就很好理解了。
对于复数i,其本身乘以i,可以理解成e12πie12πi=e12πi+12πie^{\frac{1}{2}\pi i}*e^{\frac{1}{2}\pi i}=e^{\frac{1}{2}\pi i+\frac{1}{2}\pi i}
而之前说过,复数加法群和复数乘法群通过指数函数进行对称操作的映射。
增加了12πi\frac{1}{2}\pi i可以等效为也就是旋转90度(变化12π\frac{1}{2}\pi个弧度)
结合做坐标系,不难理解为什么i2=1i^2=-1就像理解1+11+1一样自然

就像证明了一遍为什么1+1等于2 扯的有点远了,似乎已经逐渐忘记标题。

复数的单位根


时隔一周,在逐渐忘记标题之前有必要梳理一下我们现在的递归堆栈:

  1. 为了了解FFT是如何优化多项式乘法的,我们了解了多项式的表达方式
  2. 对比了多项式表达方式后,需要了解如何去优化多项式两种表达方式的互相转换
  3. FFT优化多项式,选用了复数域的特定点来构成多项式的点值表达式
  4. 为了搞清楚为什么使用复数域的特定点,我们了解复数的一些计算特点。

好了,继续单位根。
对于f(x)=zn=1f(x)=z^n=1的所有解,称做单位根。
之所以叫单位根,因为这些根都分布在模值为1都单位元上。

对于n次幂,可以理解为复数乘法群到单位元(1)在经过n次z到变换后,依然是单位元。
因此,变换但模值必然是1了,旋转弧度可以是2πkn,k{0,1,,n1}\frac{2\pi k}{n},k\in\{0,1,\dots,n-1\}
(证明很简单)

单位根有2个好特性:

  1. 能简化多项式计算: 其kn,k{N}kn,k\in \{N\}次幂结果为1
    • 是不是可以简化多项式中但一部分指数运算。
  2. 解要多少有多少:解个数和n次一致,n次就有n个解
    • 想一想系数表达式转换成点值表达式,对于最高n-1次但多项式来说,需要n个点,是不是可以直接用n次的单位根

单位根的特殊性质


我想说最重要但就是单位根的特殊性质了,FFT的优化关键所在。
此处只需要记住即可,证明也很简单
令主n次单位根(即n次下的1次单位根)ωn=e2πin\omega\\_n=e^{\frac{2\pi i}{n}}

  1. ωdndk=ωnk\omega\\_{dn}^{dk}=\omega\\_n^k(直接展开即可证明)
  2. ω2nk+n=ω2nkω2nn=ω2nk1=ω2nk\omega\\_{2n}^{k+n}=\omega\\_{2n}^k*\omega\\_{2n}^n=\omega\\_{2n}^k*{-1}=-\omega\\_{2n}^k
  3. ωnk+n2=ωnk\omega\\_{n}^{k+\frac{n}{2}} = -\omega\\_n^k
  4. (ωnk+n2)2=ωn2k+n=ωn2k=(ωnk)2(\omega\\_n^{k+\frac{n}{2}})^2=\omega\\_n^{2k+n}=\omega\\_n^{2k}=(\omega\\_n^k)^2(可从第三条得证)

FFT快速傅立叶变换


多项式能进行FFT优化的前提


一个能进行FFT的多项式,必须满足其项数必须是2的幂(或者其最高次+1的值必须得是2的幂),切勿将项数和多项式的最高次数搞混,以下说明中n指代多项式一共有n项,那么其最高次为n-1。

那么是否代表此优化方式具有局限性?当然不是,FFT被泛用于各个领域,对于项数不满足要求的多项式来说,我们只需要补充上系数为0的项,均可以满足条件)

FFT如何将问题分治和递归的


FFT真正的操作从这里开始。

ok,已经有了神奇的点了,那么可以直接把点丢到多项式里看一下计算点值表达式的时候,有没有奇迹发生了。
先来一个项数为n的多项式:简单记作G(x)=AG(x)=A,AA为系数向量。
现在计算矩阵可以表达成

[1ω n0(ω n0)2(ω n0)3(ω n0)n11ω n1(ω n1)2(ω n1)3(ω n1)n11ω n2(ω n2)2(ω n2)3(ω n2)n11ω n3(ω n3)2(ω n3)3(ω n3)n11ω nn1(ω nn1)2(ω nn1)3(ω nn1)n1][a0a1a2a3an1]=[f(ω n0)f(ω n1)f(ω n2)f(ω n3)f(ω nn1)]\begin{bmatrix}1 & \omega\ _n^0 & (\omega\ _n^0)^2 & (\omega\ _n^0)^3 & \cdots & (\omega\ _n^0)^{n-1}\\1 & \omega\ _n^1 & (\omega\ _n^1)^2 & (\omega\ _n^1)^3 & \cdots & (\omega\ _n^1)^{n-1}\\1 & \omega\ _n^2 & (\omega\ _n^2)^2 & (\omega\ _n^2)^3 & \cdots & (\omega\ _n^2)^{n-1}\\1 & \omega\ _n^3 & (\omega\ _n^3)^2 & (\omega\ _n^3)^3 & \cdots & (\omega\ _n^3)^{n-1}\\\vdots & \cdots & \cdots & \cdots & \ddots & \vdots\\1 & \omega\ _n^{n-1} & (\omega\ _n^{n-1})^2 & (\omega\ _n^{n-1})^3 & \cdots & (\omega\ _n^{n-1})^{n-1}\\\end{bmatrix}\begin{bmatrix} a_0\\a_1\\a_2\\a_3\\\vdots\\a_{n-1}\end{bmatrix}=\begin{bmatrix} f(\omega\ _n^0)\\f(\omega\ _n^1)\\f(\omega\ _n^2)\\f(\omega\ _n^3)\\\vdots\\f(\omega\ _n^{n-1})\end{bmatrix}

当然直接朴素计算时间复杂度还是O(n2)O(n^2)的。

奇偶项拆分


首先挑出矩阵的一行来看:
A(x)=a0+a1x1++an1xn1,k[0,n)A(x)=a_0+a_1x^1+\dots +a_{n-1}x^{n-1},k\in[0,n)

单独挑出奇数项和偶数项,重新组成2个项数一致的多项式(别忘记大前提,多项式项数必须是2的幂所以一定能二等分):
A(x)=A1(x)+xA2(x)A(x)=A_1(x)+x*A_2(x)

  • 奇数项多项式:
    A1(x)=a0+a2x2+a4x4++an2xn2A_1(x)=a_0+a_2x^2+a_4x^4+\dots+a_{n-2}x^{n-2}
  • 偶数项多项式:
    xA2(x)=x(a1+a3x2+a5x4++an1xn2)x*A_2(x)=x(a_1+a_3x^2+a_5x^4+\dots+a_{n-1}x^{n-2})
    A2(x)=a1x1+a3x2+a5x4++an1xn2A_2(x)=a_1x^1+a_3x^2+a_5x^4+\dots+a_{n-1}x^{n-2}

子串A1A_1以及A2A_2,变成了x的指数为等差偶数的多项式。
此处可以看到,由于分成奇偶项之后,导致其指数值不连续。如果此处能将指数值调整为连续,那么就可以用分治法继续用同样的方法去处理这俩个奇偶多项式子串。
因此想要调整为和最初多项式类似的连续指数形式:

注意以下式子中的k,n为常数(且n为奇数,k和矩阵的行数有关,n和单位根个数有关),令:

  • x=t12x=t^{\frac{1}{2}}带入A(x)A(x)

得到:

  • A1(t)=a0+a2t+a4t2++an2tn22A_1(t)=a_0+a_2t+a_4t^2+\dots+a_{n-2}t^{\frac{n-2}{2}}
  • A2(t)=a1+a3t+a5t2++an1tn22A_2(t)=a_1+a_3t+a_5t^2+\dots+a_{n-1}t^{\frac{n-2}{2}}

我们取x为单位根,那么将按照t=x2t=x^2,带单位根到多项式A(x)A(x)中:

  • A(ω nk)=A1(ω n2k)+ω nkA2(ω n2k)A(\omega\ _n^k)=A_1(\omega\ _n^{2k})+\omega\ _n^k*A_2(\omega\ _n^{2k})

从最后一个式子就可以很明显的看出来如何处理了。
利用单位根性质:

  • ωdndk=ωnk\omega\\_{dn}^{dk}=\omega\\_n^k

消除指数上的2,处理最后一个式子可以得到:

  • A(ω nk)=a0+a2ω n2k++an1(ω n2k)n22+ω nk(a1+a3ω n2k++an(ω n2k)n22),0<=k<n2A(\omega\ _n^k)=a_0+a_2\omega\ _{\frac{n}{2}}^k+\dots+a_{n-1}(\omega\ _{\frac{n}{2}}^k)^{\frac{n-2}{2}}+\omega\ _n^k(a_1+a_3\omega\ _{\frac{n}{2}}^k+\dots+a_n(\omega\ _{\frac{n}{2}}^k)^{\frac{n-2}{2}}),0<=k<\frac{n}{2}
  • A(ω nk)=A1(ω n2k)+ω nkA2(ω n2k),0<=k<n2A(\omega\ _n^k)=A_1(\omega\ _{\frac{n}{2}}^k)+\omega\ _n^k*A_2(\omega\ _{\frac{n}{2}}^k),0<=k<\frac{n}{2}

A1A_1A2A_2又变成了两个指数连续的指数级数,同理,其可以被继续二分。

通过以上的步骤,我们把同一行的点值计算,按照奇偶项拆分成2个多项式,长度减少了一半,由于n为2的幂,所以能够一直递归拆分到A1A2A_1 A_2只剩1项为止。

矩阵折半


到上节所述为止,好像没什么问题,但是在分治迭代的时候,多了一个条件k<n2k<\frac{n}{2},因为作为单位根来说,我们限制了k的范围,不会超过其单位根的总数。
显然,只能解决一半的单位根点值求解…
那么,现在我们需要看一下n>k>=n2n>k>=\frac{n}{2}的时候,需要如何解决。

很显然,我们希望k能回到[0,n2)[0,\frac{n}{2})的范围,这样我们还是可以用上面的方法进行处理。
令:
n>k>=n2n>k>=\frac{n}{2},令k=t+n2k=t+\frac{n}{2},那么n2>t>=0\frac{n}{2}>t>=0
我们直接把k进行代换,直接快进到刚刚做指数连续化处理前这一步,可得:

  • A(ω nt+n2)=A1((ω nt+n2)2)+ω nt+n2A2((ω nt+n2)2),0<=t<n2A(\omega\ _n^{t+\frac{n}{2}})=A_1((\omega\ _{n}^{t+\frac{n}{2}})^2)+\omega\ _n^{t+\frac{n}{2}}*A_2((\omega\ _n^{t+\frac{n}{2}})^2),0<=t<\frac{n}{2}

我们期望使用t而不是k来进行处理,这样,在进行下一步
利用单位根性质(ωdndk=ωnk\omega\\_{dn}^{dk}=\omega\\_n^k)进行指数连续化处理的时候
我们才能保证变量范围是合理的。
下面利用单位根的第二条性质:

  • ωnk+n2=ωnk\omega\\_{n}^{k+\frac{n}{2}} = -\omega\\_n^k
    于是,原式子可以变化为:

  • A(ω nt+n2)=A1((ω nt)2)ω ntA2((ω nt)2),0<=t<n2A(\omega\ _n^{t+\frac{n}{2}})=A_1((-\omega\ _n^t)^2)-\omega\ _n^t*A_2((-\omega\ _n^t)^2),0<=t<\frac{n}{2}

  • A(ω nk)=A(ω nt)=A1(ω n2t)ω ntA2(ω n2t),0<=t<n2,k=t+n2,n>k>=n2A(\omega\ _n^k)=A(-\omega\ _n^t)=A_1(\omega\ _{\frac{n}{2}}^t)-\omega\ _n^t*A_2(\omega\ _{\frac{n}{2}}^t),0<=t<\frac{n}{2},k=t+\frac{n}{2},n>k>=\frac{n}{2}
    去掉中间步骤,可得以下结论:

  • A(ω nt)=A1(ω n2t)ω ntA2(ω n2t),0<=t<n2,t=kn2,n>k>=n2A(\omega\ _n^t)=A_1(\omega\ _{\frac{n}{2}}^t)-\omega\ _n^t*A_2(\omega\ _{\frac{n}{2}}^t),0<=t<\frac{n}{2},t=k-\frac{n}{2},n>k>=\frac{n}{2}

  • A(ω nk)=A1(ω n2k)+ω nkA2(ω n2k),0<=k<=n2A(\omega\ _n^k)=A_1(\omega\ _{\frac{n}{2}}^k)+\omega\ _n^k*A_2(\omega\ _{\frac{n}{2}}^k),0<=k<=\frac{n}{2}
    以上结论,那一个具体的例子说明:

  • 假设有一个4项的多项式(n=4),其最高次为3次,k的范围为[0,3]:
    那么当带入的点为单位1次根(k=1)和单位3根(k=3)的时候,这两个多项式,均能拆分成2个相同的A1A_1A2A_2,很明显,我们只需要使用递归计算一半点值即可推导出剩下的那一半。

蝶形算法


那么为森么要叫蝶式算法,讲道理这名字还是挺好听的。
现在有前2节的结论后,可以看一下,多项式计算是如何被FFT优化的。

首先,假设计算的多项式A(k)A(k)是一个项数为8,带入第k次单位根的多项式结果(即n=8,k[0,7]k\in[0,7])。
我们以xn(k)x_n(k)代指这个多项式中当每一项计算结果(即已经带入单位根做指数运算并乘前面当系数)
ω (k)\omega\ (k)表示奇偶分治当时候当偶数队列提取出当公共系数变量。

  1. 理论上,我们需要计算出A(0),A(1),A(2),A(3),A(4),A(5),A(6),A(7)。
    但是按照之前的推导,我们只需要将A(0),A(1),A(2),A(3),分治奇数项和偶数项,即可得到A(4),A(5),A(6)和A(7):
    A(0)=B1(0)+ω (0)B2(0)A(0)=B_1(0)+\omega\ (0)*B_2(0)
    A(1)=B1(1)+ω (1)B2(1)A(1)=B_1(1)+\omega\ (1)*B_2(1)
    A(2)=B1(2)+ω (2)B2(2)A(2)=B_1(2)+\omega\ (2)*B_2(2)
    A(3)=B1(3)+ω (3)B2(3)A(3)=B_1(3)+\omega\ (3)*B_2(3)
    A(4)=B1(0)ω (0)B2(0)A(4)=B_1(0)-\omega\ (0)*B_2(0)
    A(5)=B1(1)ω (1)B2(1)A(5)=B_1(1)-\omega\ (1)*B_2(1)
    A(6)=B1(2)ω (2)B2(2)A(6)=B_1(2)-\omega\ (2)*B_2(2)
    A(7)=B1(3)ω (3)B2(3)A(7)=B_1(3)-\omega\ (3)*B_2(3)

  2. 观察上述子序列,我们需要求解B1(k)B_1(k)B2(k),k[0,3]B_2(k),k\in[0,3]:
    +B1(0)=C1(0)+ω (0)C2(0)B_1(0)=C_1^{'}(0)+\omega\ (0)*C_2^{'}(0)
    B1(1)=C1(1)+ω (1)C2(1)B_1(1)=C_1^{'}(1)+\omega\ (1)*C_2^{'}(1)
    B1(2)=C1(0)ω (0)C2(0)B_1(2)=C_1^{'}(0)-\omega\ (0)*C_2^{'}(0)
    B1(3)=C1(1)ω (1)C2(1)B_1(3)=C_1^{'}(1)-\omega\ (1)*C_2^{'}(1)\\
    B2(0)=D1(0)+ω (0)D2(0)B_2(0)=D_1^{'}(0)+\omega\ (0)*D_2^{'}(0)
    B2(1)=D1(1)+ω (1)D2(1)B_2(1)=D_1^{'}(1)+\omega\ (1)*D_2^{'}(1)
    B2(2)=D1(0)ω (0)D2(0)B_2(2)=D_1^{'}(0)-\omega\ (0)*D_2^{'}(0)
    B2(3)=D1(1)ω (1)D2(1)B_2(3)=D_1^{'}(1)-\omega\ (1)*D_2^{'}(1)

  3. 重复以上步骤,继续分治子序列:
    C1(0)=E1(0)+ω (0)E2(0)C_1^{'}(0)=E_1^{''}(0)+\omega\ (0)*E_2^{''}(0)
    C2(1)=F1(0)ω (0)F2(0)C_2^{'}(1)=F_1^{''}(0)-\omega\ (0)*F_2^{''}(0)\\
    D1(0)=G2(0)+ω (0)G2(0)D_1^{'}(0)=G_2^{''}(0)+\omega\ (0)*G_2^{''}(0)
    D2(1)=H2(0)ω (0)H2(0)D_2^{'}(1)=H_2^{''}(0)-\omega\ (0)*H_2^{''}(0)

  4. 到此已经分治结束,E1,E2,F1,F2,G1,G2,H1,H2E_1,E_2,F_1,F_2,G_1,G_2,H_1,H_2都已经只剩下一项,即:
    E1=x0(0)E2=x4(0)F1=x2(0)F2=x6(0)G1=x1(0)G2=x5(0)H1=x3(0)H2=x7(0)E_1=x_0(0)\\E_2=x_4(0)\\F_1=x_2(0)\\F_2=x_6(0)\\G_1=x_1(0)\\G_2=x_5(0)\\H_1=x_3(0)\\H_2=x_7(0)

总结:

  • 每一次分治,都会减少一半的计算项长度,但会分治出两个子问题。由n项计算n次变成2个n2\frac{n}{2}项计算n2\frac{n}{2}次,算上合成时候的时间,因此时间复杂度变成了O(nlogn)O(n\log n)
再仔细想一下时间复杂度的问题,从上面例子看,8项分成2个4项,2个4项又分成2个2项,2个两项又分成2个1项,总的计算量由 8*8=>2*(4*4)=>2*(2*(1*1))*(2*(1*1))似乎并没有减少。 实际上原朴素计算方法的计算量是计算n项的n个点位置的值。(复杂度n*n) 而蝶形算法中,总的来说代入单位根的A(0)的n项是必然都要算的(复杂度n) 无非通过代入一次单位根后的多项式结果能推导出代入其他单位根后的多项式结果(即A(1)到A(n-1))(这一过程的复杂度为logn)

傅立叶逆变换


有了蝶形算法,可以很快的计算单位根对应的出点值向量
[f(ω n0)f(ω n1)f(ω n2)f(ω n3)f(ω nn1)] \begin{bmatrix} f(\omega\ _n^0)\\ f(\omega\ _n^1)\\ f(\omega\ _n^2)\\ f(\omega\ _n^3)\\ \vdots\\ f(\omega\ _n^{n-1}) \end{bmatrix}
O(nlogn)O(n\log n)么做逆变换的时候,根据之前正变换的矩阵,我们求出矩阵的逆,再和f(X)f(X)做矩阵乘法,就可以重新转换出参数矩阵。
(之前的参数矩阵为范德蒙方阵,一个范德蒙方阵必然是一个满秩的矩阵,因此其拥有逆矩阵。)
以上为朴素求解法,其时间复杂度都会远远超出O(nlogn)O(n\log n)

快速傅立叶逆变换


傅立叶变换的本质是拿多项式的系数向量,左乘单位根矩阵(单位根矩阵实际上是一个范德蒙矩阵),也就是多项式求值。
而傅立叶逆变换的本质是拿多项式的点值向量,左乘单位根逆矩阵,也就是多项式插值。

由上可以看出正变换和逆变换从数学角度上看还是非常相似的,所以也希望也可以使用正变换的FFT优化方法来优化逆变换。

那么首要问题就是如何推导出逆变换矩阵。

单位根矩阵属于范德蒙矩阵,范德蒙矩阵的逆矩阵,可以根据其特定性质归纳的公式直接求解。
以下阐述其他的几种求解方式。

由此傅立叶逆变换,又转变成了一个FFT的问题,不禁让人想到计算机中的减法实现其实也是由被减数加上减数的补码实现的。其实很多问题并不是对称,而是相同的变换,处于不同的呈现。

逆变换矩阵推导-拉格朗日插值法


之前说了,傅立叶逆变换,也是属于多项式插值。那么是否可以使用拉格朗日插值法来导出这个逆矩阵?
拉格朗日插值公式如下:
f(x)=iyijixxjxixjf(x)=\sum_{i}y_i\prod_{j\ne i} {x-x_j \over x_i-x_j}

其基本思想是由开关方法(在xix_i处时,只有一项其值为yiy_i,其他项都为0,故求和后自然能得到yiy_i),来得到一个最小次且满足所有点插值的多项式(应该和傅立叶逆变换解算出的多项式是同一个)

VA=YV*A = Y(傅立叶变换)
V1VA=V1YV^{-1}*V*A = V^{-1}*Y
A=V1YA = V^{-1}*Y(傅立叶逆变换)

根据拉格朗日插值公式,其可以类比为以下矩阵乘法:
[0xxjx0xj1x1x12x13x1n1x2x22x23x2n1x3x32x33x3n1xnxn2xn3xnn][y0y1y2y3yi]\begin{bmatrix} 0 & {x-x_j \over x_0-x_j} \\ 1 & x_1 & x_1^2 & x_1^3 & \cdots & x_1^n\\ 1 & x_2 & x_2^2 & x_2^3 & \cdots & x_2^n\\ 1 & x_3 & x_3^2 & x_3^3 & \cdots & x_3^n\\ \vdots & \cdots & \cdots& \cdots & \vdots\\ 1 & x_n & x_n^2 & x_n^3 & \cdots & x_n^n\\ \end{bmatrix} \begin{bmatrix} y_0\\ y_1\\ y_2\\ y_3\\ \vdots\\ y_i \end{bmatrix}

对应到傅立叶逆变换需要求解的逆矩阵中的每一项中,显然可以得到:
(V1)ij(V^{-1})_{ij}处,对应的值为kixxkxixk\prod\limits_{k\ne i} {x-x_k \over x_i-x_k}

傅立叶级数,变换和FFT的其他应用场景


傅立叶级数


傅立叶变换


快速数论变换