DFT and FFT

tags: Linear Algebra

What is DFT

DFT for Discrete Fourier Transform,顾名思义，DFT是一种转换方式，那是什么的转换方式呢？从线性代数的角度，简单的来说就是利用单位根去从多项式的一种表达形式转换为另外一种表达形式。DFT可以说是FFT的一部分。

Root of Unity

1的n次单位根(root of unity)：\(\omega_n =\sqrt[n]{1}=e^{2\pi i/n}\)
这里不多做说明，相关知识在GSLA第九章。直观上可以想象复平面上的单位圆被均分为n份。

Representation of Polynomial：

（1）coefficient represent:
对于一个n阶多项式：\(A(x)=a_0+a_1x+a_2x^2+…+a_{n-1}x^{n-1}\)
我们可以用一个coefficient vector来表示:\((a_0,a_1,…,a_{n-1})\)

(2)point-value represent:
设\(y_k=A(x_k)\)，那么我们可以用n个点来确定一个n-1多项式：
\(\{(x_0,y_0),(x_1,y_1),…,(x_{n-1},y_{n-1})\}\)
这种用一组点来求对应的多项式的过程叫多项式插值（interpolation)。用这种表示方法去确定一个唯一的n-1阶多项式的条件是有n个\(x\)分量不同的点。

现证明：n个\(x\)分量不同的点可以唯一确定一个n-1阶多项式

分析：多项式的系数可以唯一确定一个多项式，所以我们把问题转换成：n个\(x\)分量不同的点可以唯一确定一个coefficient vector

首先我们有n个方程：
\(y_0=a_0+a_1x_0+a_2x_0^2+…+a_{n-1}x_0^{n-1}\)
\(y_1=a_0+a_1x_1+a_2x_1^2+…a_{n-1}x_1^{n-1}\)
.
.
\(y_{n-1}=a_0+a_1x_{n-1}+a_2x_{n-1}^2+…+a_{n-1}x_{n-1}^{n-1}\)
写成矩阵形式就是：
\(\begin{bmatrix}1&x_0\cdots&x_0^{n-1}\\1&x_1\cdots&x_1^{n-1}\\\vdots&\ddots&\vdots\\1&x_{n-1}\cdots&x_{n-1}^{n-1}\end{bmatrix}\) \(\begin{pmatrix}a_0\\a_1\\\vdots\\a_{n-1}\end{pmatrix}\)

=\(\begin{pmatrix}y_0\\y_1\\\vdots\\y_{n-1}\end{pmatrix}\)

留意到这个矩阵是之前GSLA第五章提到过的Vandermonde matrix，它的行列式为：\(\prod\limits_{0 \le j<i \le n-1}x_i-x_j\)

因为\(x_i \not =x_j\)，所以它的determinant不为0，即这是一个可逆矩阵，\(Aa=y\),有唯一解为\(a=A^{-1}y\)
所以n个不同点可以唯一确认一个n-1阶多项式
Q.E.D.

The DFT

刚才我们证明了，用n个\(x\)分量不同的点可以唯一确定一个n-1阶多项式，所以我们可以通过巧妙的选择这n个点来达到一些目的。DFT选择了用\(\omega_n^k\)来作为n个点的\(x\)分量。

将\(\omega_n^k\)代入刚才的矩阵中，我们就得到了一个n阶Fourier matrix：
\(\begin{bmatrix}1&1\cdots&1\\1&\omega_n\cdots&\omega_n^{n-1}\\\vdots&\ddots&\vdots\\1&\omega_n^{n-1}\cdots&\omega_n^{(n-1)^2}\end{bmatrix}\)

记作\(F_n\),\(F_n{ij}=\omega_n^{(i-1)(j-1)}\)
记\(y=(y_0,y_1,…,y_{n-1}）\)，其中\(y_k=A(\omega_n^k)=\sum\limits_{j=0}^{n-1} a_j\omega_n^{kj}\)
记作\(y=DFT_n(a)\),\(a=DFT_n^{-1}(y)\)

这就是DFT的过程，\(\Theta(n^2)\)的复杂度，但是通过\(\omega_n^k\)的一些性质还有分治法，我们可以将复杂度下降为\(\Theta(nlgn)\)，其实就是FFT(Fast Fourier Transform)

Cancelation Lemma

root of unity：\(\omega_n=e^{2\pi i/n}\)
The cancelation lemma：\(\omega_{dn}^{dk}=\omega_n^k\)

Proof:
\(\omega_{dn}^{dk}=(e^{2\pi i/dn})^{dk}\)=\(e^{2dk\pi i/dn}\)=\(e^{2k\pi i/n}\)=\(\omega_n^k\)

Halving Lemma

First，\(\omega_n^{n/2}=-1\)
直观上可以想象从复平面上（1，0）的位置沿单位圆转了半圈，落到（-1，0）这个点。

So，halving lemma:
\(\omega_n^{k+n/2}=\omega_n^k \omega_n^{n/2}=-\omega_n^k\)

From DFT to FFT

FFT用了分治的策略，对于多项式，分成奇偶次项是比较自然的思路,为了方便说明，我们设n是2的幂。
现在，对于原来的多项式：
\(A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}\)
分奇偶次项：
\(A^{even}(x)=a_0+a_2x^2+a_4x^4+...+a_{n-2}x^{n-2}\)
\(A^{odd}(x)=a_1x+a_3x^3+a_5x^5+...+a_{n-1}x^{n-1}\)
我们有：\(A(x)=A^{even}(x)+A^{odd}(x)\) （1）

要达到分治的目的，我们希望把刚才拆分出的两个子问题变形为与原问题有相同的形式的问题，如：
\(A^{even}(\sqrt{x})=a_0+a_2x+a_4x^2+...+a_{n-2}x^{n-2/2}\),记作\(A^{[0]}(x)\)
\(A^{even}(\sqrt{x})\)和\(A(x)\)相比，形式上除了最高幂次是\(A(x)\)的一半且系数是\(A(x)\)的偶次系数之外，没有其他的变化，我们现在想将\(A^{odd}(x)\)也转变成这样的形式。

观察发现，\(A^{odd}(x)/x\)和\(A^{even}(x)\)有类似的结构，所以我们将它变形为：
\(A^{odd}(\sqrt{x})/x=a_1+a_3x+a_5x^2+...+a_{n-1}x^{n-2/2}\),记作\(A^{[1]}(x)\)

\(\because A^{even}(\sqrt{x})=A^{[0]}(x)\)
\(\therefore A^{even}(x)=A^{[0]}(x^2)\) （2）

\(\because A^{odd}(\sqrt{x})/x=A^{[1]}(x)\)
\(\therefore A^{odd}(x)=xA^{[1]}(x^2)\) （3）

根据式子（1）：\(A(x)=A^{even}(x)+A^{odd}(x)\)
代入（2）、（3）：
\(A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)\) （4）

在DFT中我们用\(\omega_n^k\)(\(1,\omega_n,\omega_n^2,...,\omega_n^{n-1}\))这n个点来进行插值，根据式子（4），我们在\(A^{[1]}(x)\)和\(A^{[0]}(x)\)中应该用\(\omega_n^{2k}\)(\(1,\omega_n^2,\omega_n^4,..,\omega_n^{2(n-1)}\))来进行插值。

又根据cancelation lemma:\(\omega_{dn}^{dk}=\omega_n^k\),上下标同除以2:
\(1,\omega_n^2,\omega_n^4,..,\omega_n^n,...\omega_n^{2(n-1)}\)=\(1,\omega_{n/2},\omega_{n/2}^2,...,\omega_{n/2}^{n/2-1},\omega_{n/2}^{n/2},...,\omega_{n/2}^{n-1}\)

对于\(\omega_{n/2}^{n/2}\)之后的式子，都有\(\omega_{n/2}^{n/2+k}\)的形式 \((k=1,2,3,..,n/2-1)\)
\(\omega_{n/2}^{n/2+k}=\omega_{n/2}^{n/2} \omega_{n/2}^k=\omega_{n/2}^k\)

所以用\(\omega_n^{2k}\)插值实际上就是用\(\omega_{n/2}^k\)进行插值，又因为\(A^{[0]}(x)\)和\(A^{[1]}(x)\)都是\(n/2\)阶多项式，所以我们发现这其实就是两次\(n/2\)的DFT，我们成功地把原来的问题分成了两个size为原来的一半的子问题。

*Prove of Vandermonde Determinant

-TBC-