离散傅里叶变换

引言

设函数$f(x)$在区间$[0,2\pi)$上$n$个等分点$x_k=2\pi k/n$（$k=0,1,\dots,n-1$）上的值$f_k=f(x_k)=f(2\pi k/n)$（一般为复数）为已知。在本文中，我们令$i^2=-1$。

插值问题：用周期函数

$$\{e^{ijx}=\cos jx+i\sin jx\}_{j=0,1,\dots,n-1}$$
的线性组合
$$P(x)=\sum^{n-1}_{j=0}c_je^{ijx}$$
作为 $f(x)$在这区间上的三角函数插值，即求系数 $c_j$，使得
$$f_k=\sum^{n-1}_{j=0}c_je^{i\frac{2jk\pi}{n}},\quad k=0,1,\dots,n-1$$

---

离散傅里叶变换

以上插值问题在理论上容易求解，因为它等价于求解线性方程组

$$\left(\begin{array}{cccccc} 1 & 1 & \cdots & 1 & \cdots & 1 \\ 1 & e^{i\frac{2\pi}{n}} & \cdots & e^{i\frac{2j\pi}{n}} & \cdots & e^{i\frac{2(n-1)\pi}{n}} \\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots\\ 1 & e^{i\frac{2k\pi}{n}} & \cdots & e^{i\frac{2jk\pi}{n}} & \cdots & e^{i\frac{2(n-1)k\pi}{n}} \\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots\\ 1 & e^{i\frac{2(n-1)\pi}{n}} & \cdots & e^{i\frac{2j(n-1)\pi}{n}} & \cdots & e^{i\frac{2(n-1)^2\pi}{n}} \\ \end{array}\right)\left(\begin{array}{c} c_0 \\ c_1 \\ \vdots \\ c_j \\ \vdots \\ c_{n-1} \end{array}\right)=\left(\begin{array}{c} f_0 \\ f_1 \\ \vdots \\ f_k \\ \vdots \\ f_{n-1} \end{array}\right)$$
记系数矩阵为
$$F^{-1}=\left(\begin{array}{cccccc} 1 & 1 & \cdots & 1 & \cdots & 1 \\ 1 & e^{i\frac{2\pi}{n}} & \cdots & e^{i\frac{2j\pi}{n}} & \cdots & e^{i\frac{2(n-1)\pi}{n}} \\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots\\ 1 & e^{i\frac{2k\pi}{n}} & \cdots & e^{i\frac{2jk\pi}{n}} & \cdots & e^{i\frac{2(n-1)k\pi}{n}} \\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots\\ 1 & e^{i\frac{2(n-1)\pi}{n}} & \cdots & e^{i\frac{2j(n-1)\pi}{n}} & \cdots & e^{i\frac{2(n-1)^2\pi}{n}} \\ \end{array}\right)$$
（注：由于这个矩阵对应于傅里叶逆变换，因此记为逆的形式。）显然，系数矩阵 $F^{-1}$是一个对称的范德蒙德(Vandermonde)矩阵。

性质1 以上系数矩阵的行（或列）向量具有如下的正交性：

$$\sum^{n-1}_{j=0}e^{i\frac{2jk\pi}{n}}\cdot e^{-i\frac{2jl\pi}{n}}=\begin{cases}n, & \mbox{当}k=l\\0, & \mbox{当}k\neq l\\\end{cases}$$

由性质1可知，

$$F=\frac{1}{n}\left(\begin{array}{cccccc} 1 & 1 & \cdots & 1 & \cdots & 1 \\ 1 & e^{-i\frac{2\pi}{n}} & \cdots & e^{-i\frac{2k\pi}{n}} & \cdots & e^{-i\frac{2(n-1)\pi}{n}} \\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots\\ 1 & e^{-i\frac{2j\pi}{n}} & \cdots & e^{-i\frac{2kj\pi}{n}} & \cdots & e^{-i\frac{2(n-1)j\pi}{n}} \\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots\\ 1 & e^{-i\frac{2(n-1)\pi}{n}} & \cdots & e^{-i\frac{2k(n-1)\pi}{n}} & \cdots & e^{-i\frac{2(n-1)^2\pi}{n}} \\ \end{array}\right)$$
从而得到
$$c_j=\frac{1}{n}\sum^{n-1}_{k=0}f_ke^{-i\frac{2jk\pi}{n}},\quad j=0,1,\dots,n-1$$

习惯上，由$\{f_k\}$求$\{c_j\}$的过程称为$f(x)$的离散傅里叶变换（也称$\{c_j\}$是$\{f_k\}$的离散傅里叶变换），而由$\{c_j\}$求$\{f_k\}$的过程称为离散傅里叶逆变换。用上面公式计算傅里叶变换的计算复杂度为$O(n^2)$。另外，应该注意到，$F$（或$F^{-1}$）的$n^2$个元素$e^{-i\frac{2kj\pi}{n}}$（或$e^{i\frac{2jk\pi}{n}}$）中，实际上只有$n$个不同的值。这是下面进行快速傅里叶变换的基础。

性质2（Parseval等式）

$$\sum^{n-1}_{j=0}|c_j|^2=\frac{1}{n}\sum^{n-1}_{k=0}|f_k|^2$$

---

快速傅里叶变换

在这一节，我们假设$n=2^t$，$t$为正整数。令$\omega=e^{-i\frac{2\pi}{n}}$，则$\omega^n=1$。离散傅里叶变换需要计算

$$\left(\begin{array}{c} c_0 \\ c_1 \\ \vdots \\ c_{n-1} \end{array}\right)=F\left(\begin{array}{c} f_0 \\ f_1 \\ \vdots \\ f_{n-1} \end{array}\right)$$
即
$$c_j=\frac{1}{n}\sum^{n-1}_{k=0}f_k\omega^{jk},\quad j=0,1,\dots,n-1$$

基本原理：根据$\omega$的周期性（即$\omega^n=1$），合并上式中的同类项，将$n$个乘积的和转化为$n/2$个乘积的和，从而提高计算的速度。快速傅里叶变换的计算复杂度为$O(n\log_2 n)$。

为了简化后面的讨论，我们记$a_k=\frac{f_k}{n}$（$k=0,1,\dots,n-1$），以及

$$F_{k,\omega}=\left(\begin{array}{cccccc} 1 & 1 & \cdots & 1 & \cdots & 1 \\ 1 & \omega^{1} & \cdots & \omega^{k} & \cdots & \omega^{n-1} \\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots\\ 1 & \omega^{j} & \cdots & \omega^{jk} & \cdots & \omega^{j(n-1)}\\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots\\ 1 & \omega^{(n-1)} & \cdots & \omega^{(n-1)k} & \cdots & \omega^{(n-1)^2} \\ \end{array}\right)$$
则离散傅里叶变换转化为计算
$$\left(\begin{array}{c} c_0 \\ c_1 \\ \vdots \\ c_{n-1} \end{array}\right)=F_{k,\omega}\left(\begin{array}{c} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{array}\right)$$

计算过程

分别考虑$c_j$的下标$j$是偶数和奇数的情况：

• 当$j=2l$时
  

  $$\begin{eqnarray} c_{2l}&=&\sum^{n-1}_{k=0}a_k\omega^{2lk}=\sum^{\frac{n}{2}-1}_{k=0}a_k\omega^{2lk}+\sum^{n-1}_{k=\frac{n}{2}}a_k\omega^{2lk}\\ &=&\sum^{\frac{n}{2}-1}_{k=0}a_k\omega^{2lk}+\sum^{\frac{n}{2}-1}_{k=0}a_{\frac{n}{2}+k}\omega^{2l(\frac{n}{2}+k)} \end{eqnarray}$$
  由 $\omega^{2l(\frac{n}{2}+k)}=\omega^{nl}\cdot\omega^{2lk}=1\cdot \omega^{2lk}=\omega^{2lk}$，上式转化为
  $$c_{2l}=\sum^{\frac{n}{2}-1}_{k=0}a_k\omega^{2lk}+\sum^{\frac{n}{2}-1}_{k=0}a_{\frac{n}{2}+k}\omega^{2lk}=\sum^{\frac{n}{2}-1}_{k=0}(a_k+a_{\frac{n}{2}+k})\omega^{2lk}$$

  • 当$j=2l+1$时
    

    $$c_{2l+1}=\sum^{n-1}_{k=0}a_k\omega^{(2l+1)k}=\sum^{n-1}_{k=0}(a_k\omega^{k})\omega^{2lk}$$
    我们将上式中的$a_k\omega^{k}$看作一个整体，则转化为偶数情况。从而得到
    $$c_{2l+1}=\sum^{\frac{n}{2}-1}_{k=0}(a_k\omega^{k}+a_{\frac{n}{2}+k}\omega^{\frac{n}{2}+k})\omega^{2lk}$$
    由$\omega^{\frac{n}{2}+k}=\omega^{\frac{n}{2}}\cdot\omega^{k}=-1\cdot\omega^{k}=-\omega^{k}$，上式转化为
    $$c_{2l+1}=\sum^{\frac{n}{2}-1}_{k=0}(a_k\omega^{k}-a_{k}\omega^{k})\omega^{2lk}=\sum^{\frac{n}{2}-1}_{k=0}[(a_k-a_{\frac{n}{2}+k})\omega^{k}]\cdot\omega^{2lk}$$

  • 综上，我们得到
    

    $$\begin{eqnarray} c_{2l}&=&\sum^{\frac{n}{2}-1}_{k=0}(a_k+a_{\frac{n}{2}+k})\omega^{2lk}\\ c_{2l+1}&=&\sum^{\frac{n}{2}-1}_{k=0}[(a_k-a_{\frac{n}{2}+k})\omega^{k}]\cdot\omega^{2lk} \end{eqnarray}$$

    令$\omega_1=\omega^2$，则$\omega^{\frac{n}{2}}_1=1$。记$b^{(0)}_k=a_k+a_{\frac{n}{2}+k}$，$b^{(1)}_k=(a_k-a_{\frac{n}{2}+k})\omega^{k}$。则上面的两个式子化为
    

    $$c_{2l}=\sum^{\frac{n}{2}-1}_{k=0}b^{(0)}_k\omega^{lk}_1,\quad c_{2l+1}=\sum^{\frac{n}{2}-1}_{k=0}b^{(1)}_k\cdot\omega^{lk}_1,\quad l=0,1,\dots,\frac{n}{2}-1$$
    从而得到两个离散傅里叶变换。记
    $$F_{k-1,\omega_1}=\left(\begin{array}{cccc} 1 & 1 & \cdots & 1 \\ 1 & \omega^{1}_1 & \cdots & \omega^{\frac{n}{2}-1}_1\\ \vdots & \vdots & \ddots & \vdots \\ 1 & \omega^{(\frac{n}{2}-1)}_1 & \cdots & \omega^{(\frac{n}{2}-1)^2} \\ \end{array}\right)$$
    则两个两个离散傅里叶变换分别为
    $$\begin{eqnarray} \left(\begin{array}{l} c_0 \\ c_2 \\ \vdots \\ c_{2(\frac{n}{2}-1)} \end{array}\right)&=&F_{k-1,\omega_1}\left(\begin{array}{l} b^{(0)}_0 \\ b^{(0)}_1 \\ \vdots \\ b^{(0)}_{\frac{n}{2}-1} \end{array}\right)\\ \left(\begin{array}{l} c_1 \\ c_3 \\ \vdots \\ c_{2(\frac{n}{2}-1)+1} \end{array}\right)&=&F_{k-1,\omega_1}\left(\begin{array}{l} b^{(1)}_0 \\ b^{(1)}_1 \\ \vdots \\ b^{(1)}_{\frac{n}{2}-1} \end{array}\right) \end{eqnarray}$$
    由于 $\frac{n}{2}=2^k/2=2^{k-1}$，我们可以递归地使用前面的过程完成计算。这种计算 $c_j$的过程称为 快速傅里叶变换。

    计算复杂度

    假设离散傅里叶变换需要$M_k$次乘法和$A_k$次加法。从$2^k$个数据的傅里叶变换转换为两个$2^{k-1}$个数据的傅里叶变换，需要额外的$2^{k-1}$个乘法和$2^k$个加法，因此，我们得到递推关系
    

    $$\begin{eqnarray} M_k&=&2M_{k-1}+2^{k-1}\\ A_k&=&2A_{k-1}+2^k \end{eqnarray}$$
    由 $M_0=0$和 $A_0=0$，得到
    $$M_k=k\cdot 2^{k-1}=\frac{1}{2}n\log_2 n,\quad A_k=k\cdot 2^{k}=n\log_2 n$$

    关于$\omega^k$的计算
    在作快速傅里叶变换时，（在额外的乘法运算中）要用到
    

    $$\omega^k=e^{-i\frac{2k\pi}{n}}=\cos\frac{2k\pi}{n}-i\sin\frac{2k\pi}{n},\quad l=0,1,\dots,\frac{n}{2}-1$$
    我们可以用如下步骤计算 $c_k=\cos\frac{2k\pi}{n}$和 $s_k=\sin\frac{2k\pi}{n}$。
    步骤1：令 $c_0=1$和 $s_0=0$；
    步骤2：用标准库函数计算出 $c_1=\cos\frac{2\pi}{n}$和 $s_1=\sin\frac{2\pi}{n}$；令 $k=1$；
    步骤3：计算 $c_{k+1}=c_{k}c_1-s_ks_1$和 $s_{k+1}=s_kc_1-c_ks_1$；令k=k+1，重复步骤3。

    实验代码（待补充）