FFT简陋入门

FFT入门

FFT的用途

在\(\,\Theta(n\log{n})\,\)的时间内计算离散傅里叶变化(DFT),通常用来计算多项式乘法

点值表达式

引理1:任何\(\,n-1\,\)次多项式可以由其在\(\,n\,\)个点的取值唯一确定

考虑反证,设\(\,n\,\)个点\(\,a_1,a_2\cdots a_n\)同时被两个\(\,n-1\,\)次多项式函数\(\,A,B\,\)经过,令

\[C(x)=A(x)-B(x)\rightarrow \forall x_0\in a\,,\,C(x_0)=A(x_0)-B(x_0)=0 \]

与代数基本定理矛盾,假设不成立,得证

那么,可以用点值表达式唯一确定一个多项式,而点值表达式的乘积可以在\(\,\Theta(n)\,\)的时间内算出,就是对应点的点值相乘

关于复数

如果您学过复数,请略过

我们定义复数为形如\(\,a+bi\,\)的数,其中\(\,a\,\)叫实部,\(\,b\,\)叫虚部

定义复数的运算:

\[(a+bi)+(c+di)=(a+b)+(c+d)i \]

\[(a+bi)-(c+di)=(a-c)+(b-d)i \]

\[(a+bi)\times(c+di)=(ac-bd)+(ad+bc)i \]

\[\frac{a+bi}{c+di}=\frac{ac+bd}{c^2+d^2}+\frac{bc-ad}{c^2+d^2}i \]

不难发现这些运算符合交换律结合律分配律

定义复数的模长\(\,|z|\,\)和辐角\(\,\theta\,\):

\[z=a+bi\rightarrow |z|=\sqrt{a^2+b^2} \]

\[sin\theta=\frac{b}{\sqrt{a^2+b^2}}\,\,,\,\,cos\theta=\frac{a}{\sqrt{a^2+b^2}} \]

定义\(\,z\,\)的共轭复数\(\,z^\prime\,\):

\[z=a+bi\rightarrow z^\prime=a-bi \]

结合定义,有

\[|z|=|z^\prime|,z·z^\prime=a^2+b^2 \]

把复数\(\,z=a+bi\,\)映射为平面上的点\(\,(a,b)\,\),可以发现复数\(\,a+bi\,\)和向量\(\,(a,b)\,\)一一对应

考虑两者之间的关联,发现加减规则相同,向量内积的几何意义是平行四边形的对角线长度,而复数乘意为模长相乘,辐角相加,这一点结合复数的运算不难推出

单位根

定义\(\,n\,\)次单位根\(\,\omega\,\)为满足\(\,\omega^{n}=1\,\)的复数\(\,\omega\),发现\(\,n\,\)次单位根对应平面上单位圆的\(\,n\,\)等分点

设\(\,\omega_{n}^{k}\,\)表示辐角第\(\,k\,\)大的\(\,n\,\)次单位根,联想复数乘法的几何意义,不难发现:

\[w_{n}^k=w_{n}^{k\,\,mod\,\,n}\,\,,\,\,w_n^{a+b}=w_n^a+w_n^b\,\,,\,\,w_n^n=w_n^0=1 \]

运用欧拉公式\(\,e^{i\theta}=isin\theta+cos\theta\,\),可得\(w_n^k=e^{2\pi\frac{k}{n}i}\),由此可以得到下方的两条引理:

引理2(消去引理):\(w_{mn}^{mk}=w_{n}^k\).

引理3(折半引理):\(w_n^{k+\frac{n}{2}}=-w_n^k\).

引理2证明如下:

\[w_{mn}^{mk}=e^{2\pi\frac{mk}{mn}i}=e^{2\pi\frac{k}{n}i}=w_n^k \]

引理3证明如下:

\[w_{n}^{k+\frac{n}{2}}=e^{2\pi\frac{k+\frac{n}{2}}{n}i}=e^{2\pi\frac{k}{n}i}e^{2\pi\frac{1}{2}i}=-e^{2\pi\frac{k}{n}i}=-w_n^k \]

单位根与FFT

考虑到FFT有两个难点,一是系数表达式转点值表达式,二是反向转化

正向变换

先考虑用单位根优化第一步,令偶次多项式\(\,A,B\,\)为

\[A(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}=\sum_{i=0}^{n-1}a_ix^i \]

\[B(x)=b_0+b_1x+b_2x^2+\cdots+b_{n-1}x^{n-1}=\sum_{i=0}^{n-1}b_ix^i \]

这里如果次数不足,系数用\(\,0\,\)补足

不妨讨论\(\,A\,\),令

\[A1(x)=a_0+a_2x+a_4x^2+\cdots+a_{n-2}^{\frac{n-2}{2}}=\sum_{i=0}^{\frac{n-2}{2}}a_{2i}x^i \]

\[A2(x)=a_1+a_3x+a_5x^2+\cdots+a_{n-1}x^{\frac{n-2}{2}}=\sum_{i=0}^{\frac{n-2}{2}}a_{2i+1}x^i \]

于是有:

\[A(x)=A1(x^2)+xA2(x^2) \]

将\(\,n\,\)次单位根一一带入\(\,A\,\)有

若\(\,0\leq k \leq \frac{n}{2}-1\),则

\[A(w_n^k)=A1(w_n^{2k})+w_n^kA2(w_n^{2k})=A1(w_{\frac{n}{2}}^k)+w_n^kA2(w_{\frac{n}{2}}^k) \]

若\(\,\frac{n}{2}\leq k^\prime \leq n-1\),记\(\,k^\prime=k+\frac{n}{2}\,\)则

\[A(w_n^{k+\frac {n}{2}})=A1(w_n^{2k+n})+w_n^{k+\frac {n}{2}}A2(w_n^{2k+n})=A1(w_{\frac{n}{2}}^k)-w_n^kA2(w_{\frac{n}{2}}^k) \]

显然,知道\(\,A1,A2\,\)的值之后,我们可以\(\,\Theta(1)\,\)求出\(\,A\,\)的取值,这样操作\(\,n\,\)次,我们就得到了\(\,A\,\)的点值表达式

\(A1,A2\)是两个\(\frac{n}{2}\)次的多项式,可以递归求解,于是复杂度为:

\[T(n)=2T(\frac{n}{2})+\Theta(n)=\Theta(n\log n) \]

反向变换

再考虑将点值表达式转化为系数表达式,这个过程叫做离散傅里叶反变换

不妨记原多项式函数为\(\,F(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}=\sum_{i=0}^{n-1}a_ix^i\)

令\(\,A(x)=\sum_{i=0}^{n-1}d_ix^i\,\)其中\(\,d\,\)是\(\,a\,\)离散傅里叶变换的结果,即我们对\(\,a\,\)进行正向变换从而得到\(\,d\,\)

\[c_k=A(w_n^{-k})=\sum_{i=0}^{n-1}d_i(w_n^{-k})^i \]

考虑将\(\,d\,\)展开,得

\[c_k=\sum_{i=1}^{n-1}(\sum_{j=0}^{n-1}a_j(w_n^i)^j)(w_n^{-k})^i \]

为了利用单位根的性质,考虑交换求和号

\[c_k=\sum_{i=1}^{n-1}(\sum_{j=0}^{n-1}a_j(w_n^i)^j)(w_n^{-k})^i=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(w_n^i)^j(w_n^{-k})^i=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(w_n^i)^{j-k} \]

\[S(j,k)=\sum_{i=0}^{n-1}(w_n^i)^{j-k} \]

不难看出这是一个等比数列,而且有\(S(i,i)=n\),接着考虑\(j\neq k\)的情况

直接运用等比数列求和公式:

\[S(j,k)=\sum_{i=0}^{n-1}(w_n^i)^{j-k}=\frac{w_n^0(w_n^{n(j-k)}-1)}{w_n^{j-k}-1}=0 \]

可以注意到上式中\(\,w_n^n\equiv1\,\),那么便不难计算

整理得:

\[S(j,k)=n\times[j\neq k] \]

将\(\,S(j,k)\,\)带回\(\,c_k\,\)的表达式,得

\[c_k=\sum_{j=0}^{n-1}a_jS(j,k)=na_k \]

考虑用\(\,c\,\)来表达\(F(x)\),于是有

\[F(x)=\frac{1}{n}\sum_{i=0}^{n-1}c_ix^i \]

此时就得到了大致的思路:

正向变换时已经得到了\(\,d\),就是点值两两相乘的结果

由\(c_k=A(w_n^{-k})=\sum_{i=0}^{n-1}d_i(w_n^{-k})^i\)不难看出\(\,c\,\)是\(\,d\,\)进行一次离散傅里叶变换的结果

最后有

\[a_k=\frac{c_k}{n} \]

综上,我们实现了在\(\,\Theta(n\log n)\,\)的时间内计算多项式乘法

代码实现

//等我会写了再上传

上一篇:Js - 垃圾回收


下一篇:MATLAB矩阵选取某行(列)数据