FFT

note : 本文仅简要地介绍FFT以及它的简单应用, 并不会过多的进行数学推导.

Background

首先我们回忆一下傅里叶级数.

实数上的情形

对于一个函数\(f(x)\), 我们现在关注它在\([-\pi,\pi]\) 这个区间上的表现.

它能够写成(在一定条件下) \[ f(x) = \frac{a_0}{2} + \sum_{k=1}^{\infty} a_kcos(kx)+b_ksin(kx) \] 这里 \[ \begin{align*} a_k = \frac{1}{\pi}\int_{-\pi}^{\pi} f(x)cos(kx)dx \\\\ b_k = \frac{1}{\pi}\int_{-\pi}^{\pi} f(x)sin(kx)dx \end{align*} \] 虽然这里我们考虑的是在区间\([-\pi,\pi]\) 上, 但现在 \(f\) 现在在整个\(\textbf{R}\) 具有了周期性, 周期是\(2\pi\).

另一方面, \(f(x)\) 要满足"平方可积", 详细来说, 所有平方可积的函数构成一个线性空间. 更进一步的有关希尔伯特空间等泛函分析的内容.

那这个线性空间自然具有一个无限维的基. 在傅里叶级数中, 选择的基是 \[ \\{ 1, cos(x), sin(x), cos(2x), sin(2x), \dots \\} \] 在确定内积 \[ \left \langle f,g\right \rangle = \int_{-\pi}^{\pi}f(x)\overline{g(x)}dx \] 后,这组由 \(cos(kx),sin(kx)\dots\) 构成的基成为了一组正交基. 当然你可以选择规范化, 去除以每个的范数.

复数上的情形

把实数上的情形推广到复数上, 这里是说\(f(x)\) 现在是复数函数, (x依然是实数).

我们有了更一般一点的傅里叶级数 \[ f(x) = \sum_{k=-\infty}^{+\infty} c_k e^{ikx} \] 这里 \[ c_k = \frac{\left \langle f(x), e^{ikx} \right\rangle}{\left \|e^{ikx} \right \|^2} = \frac{1}{2\pi}\int_{-\pi}^{\pi} f(x)e^{-ikx}dx \] (注: Euler公式 \(e^{ix}=cos(x) + i sin(x) ,\quad i=\sqrt{-1}\))

仍然是用许多的\(cos(kx), sin(kx)\)去表达\(f(x)\), 只不过变成了复数.

前面一直讲的是\(f(x)\)在\([-\pi,\pi]\)上的一些事情, 实际上这并不是唯一的.

现在换成区间\([-L,L]\) 做个映射, \(x \to \pi x/L\)

记 \[ \psi_k(x)=e^{i\pi kx/L} \]

有 \[ c_k = \frac{\left \langle f(x), \psi_k(x) \right\rangle}{\left \|\psi_k(x) \right \|^2} = \frac{1}{2L}\int_{-L}^{L} f(x)\psi_k(x)dx \]

所以 \[ f(x) = \sum_{k=-\infty}^{\infty}c_k\psi_k(x)=\frac{1}{2L}\sum_{k=-\infty}^{\infty}\left \langle f(x), \psi_k(x) \right\rangle \psi_k(x) \]

这里观察这个级数, 它实际是在将\(f(x)\) 分解成一系列\(\psi_k(x)\) 的线性组合. 或者说一系列\(cos(2\pi kx/L), sin(2\pi kx/L)\) 的结合. 这里面\(k\) 蕴含着不同的正弦函数的频率.

记 \[ \omega_k = \pi k/L \\\\ \Delta \omega = \pi/L \] 那么 \[ \begin{align*} c_k = \frac{1}{2L} \int_{-L}^{L} f(x)& e^{-ik{\Delta\omega}x}dx = \frac{1}{2L}\int_{-L}^{L} f(x)e^{-i\omega_kx}dx \\\\ f(x)=& \sum_{k=-\infty}^{+\infty}c_k e^{i\omega_k x}=\sum_{k=-\infty}^{+\infty}c_k e^{ik\Delta \omega x} \end{align*} \]

每个\(c_k\)都对应着一个\(\omega_k\), 表示这个频率的函数所占的比重.

当\(L \to +\infty\), 也就是说\(f(x)\) 在\(\textbf{R}\) 上非周期. 我们就得到了傅里叶变换. \[ \begin{align*} \hat f(\omega) & = \mathcal{F(f(x))}=\int_{-\infty}^{\infty} f(x)e^{-i\omega x}dx\\\\ f(x) =& \mathcal{F^{-1}(\hat f(\omega))} = \frac{1}{2\pi}\int_{-\infty}^{\infty} \hat f(\omega)e^{i\omega x}d\omega \end{align*} \]

有趣的性质:

\[ \begin{align*} & \mathcal{F(\frac{d}{dx}f(x)) = i\omega F(f(x))}\\\\\\\\ & \mathcal{F(g\times f) = F(g)F(f)} \end{align*} \]

DFT

对于计算机, 我们既不能方便的处理无穷级数, 也不好处理积分. 连续的情况对计算机来说比较困难. 所以出现了离散的傅里叶变换. 他的大致含义还是一样的.

对一个\(f(x)\)进行采样, 在0, 1, 2, ... , n - 1处的值, 分别为\([f_0, f_1,f_2,...,f_{n-1}]\) , 这个向量代替了原来的函数

经过傅里叶变换后的函数变成了 \([\hat f_0, \hat f_1, \hat f_2...\hat f_{n-1}]\)

对于 \(\omega\) 的选取也有原先连续的情况变成离散. 将 \(2\pi\) 平均分成了 n 份. 每一份长度是 \(\frac{2\pi}{n}\) 作为一个基本的单元.

之前的函数 \(\psi_k(x)\) 成为了 \([(e^{2\pi ik/n})^j], j=0,1,2,...,n-1\)

简单来说, DFT的形式就是 \[ \begin{align*} \hat f_k=\sum_{j=0}^{n-1}f_j e^{-2\pi ijk/n}\\\\ f_k = \frac{1}{n}\sum_{j=0}^{n-1}\hat f_k e^{2\pi ijk/n} \end{align*} \]

注意 \(e^{2\pi ik/n},\quad k=0,1,2...,n-1\) 就是 \(x^{n}=1\) 的解

\(\psi_k = \left[ \begin{matrix} 1 \\\\ e^{2\pi ik/n} \\\\ (e^{2\pi ik/n})^2 \\\\ \vdots \\\\ (e^{2\pi ik/n})^{n-1} \end{matrix} \right]\) 可以看作是不同频率的函数, 那么\(\hat f_k\) 便可以看作是对应的\(\psi_k\) 所占的大小.

首先举个简单的例子去直观的理解这个变换.

我们在函数\(f(x)=cos(2\pi \times 50 x)+cos(2\pi \times 75x)\) 上进行采样, 同时我们加上一个随机的噪声,范围是\([-4,4]\)

黑色的是原图像, 蓝色的部分是加入噪声后的图象. 当然这里的数据都是离散的, 只不过画图的时候被连在了一起.看上去像连续的.

现在我们得到了采样后的一个\(f(x)\)的向量, \([f_0,f_1,f_2,...,f_{n-1}]\)

经过DFT后, 我们能得到一个\(\hat f\) 向量, 由于是复数, 这里只展示向量中每个复数的模长的平方. 可以看作是一种能量.

下面的图是经过DFT后的样子, 明显可以看到, 里面有两个非常高的点. 对应着原始的没有加入噪声的\(cos(2\pi \times 50 x)\) 和 \(cos(2\pi \times 75x)\) , 如果我们设置一个阈值, 小于100的置位0, 这样就能去除噪声.

中间的图就是对\(\hat f\) 进行过滤之后再经过DFT的逆操作的得到的图象, 会发现它完全去除了噪声, 还原了最初的\(f(x)\) .

计算DFT

下面的问题是如何去计算DFT呢, \[ \begin{align*} \hat f_k=\sum_{j=0}^{n-1}f_j e^{-2\pi ijk/n}\\\\ f_k = \frac{1}{n}\sum_{j=0}^{n-1}\hat f_k e^{2\pi ijk/n} \end{align*} \] 容易看到这实质上是一个矩阵乘法.

记 \(\omega_n = e^{-2\pi i/n}\) \[ \hat f=\left [ \begin{matrix} \hat f_0 \\\\ \hat f_1 \\\\ \hat f_2 \\\\ \vdots \\\\ \hat f_{n-1} \end{matrix} \right ] = \left [ \begin{matrix} 1 & 1 & 1 & 1 & \dots & 1 \\\\ 1 & \omega_n & \omega_n^{2} & \omega_n^{3} & \dots & \omega_n^{n-1}\\\\ 1 & \omega_n^{2} & \omega_n^{4} & \omega_n^{6} & \dots & \omega_n^{2(n-1)} \\\\ & & \vdots & & \dots &\vdots \\\\ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \omega_n^{3(n-1)} & \dots & \omega_n^{(n-1)(n-1)} \end{matrix}\right] \left[ \begin{matrix} f_0 \\\\ f_1 \\\\ f_2 \\\\ \vdots \\\\ f_{n-1} \end{matrix}\right] \] 记这个矩阵为\(\mathcal{F_n}\) , 或者称为DFT矩阵.

实质上\(\mathcal{F}\) 乘上它的共轭转置 \(\mathcal{F^*}\) 等于 \(nE_n\)， 所以它的逆矩阵是很好求的，DFT的逆变换也很容易得出.

如果直接去求这个矩阵，复杂度是比较高的\(O(n^2)\)

但是如果考虑到 \(\omega_n\) 的特殊的性质，便能得到非常高效的\(O(nlogn)\) 的算法，也就是FFT

FFT

FFT就是一种计算DFT的高效的算法，它最初是由高斯提出.

简单来说，我们想要快速的计算矩阵\(\mathcal{F_n}\) ，通过将他分解为更小的子问题，转而去求解\(\mathcal{F_{\frac{n}{2}}}\). 这里我们现在只考虑 n 是2的幂的情况. 如果n不是2的幂的话， 我们可以简单的扩充成2的幂， 并不会影响算法的效率.

不同的矩阵\(\mathcal{F_n}\) 内所使用的\(\omega_n\) 是不同的， 但对于n 和 n/2 我们有 \(\omega_n^2=\omega_{n/2}\) 可以轻松的进行转化。

核心的想法就是对\([f_0, f_1,...,f_{n-1}]\) 进行重新排序， 按照下标的奇偶进行分组.

对于多项式 \[ \begin{align*} p(x) =& a_0+a_1x+a_2x^2+\dots+a_{n-1}x^{n-1} \\\\ =& (a_0 + a_2x^2 + a_4x^4 + \dots + a_{n-2}x^{n-2})\\\\ &+ x(a_1 + a_3x^2+a_5x^4+\dots a_{n-1}x^{n-2})\\\\ =& E(x) +xO(x) \end{align*} \]

同理，对于这里的DFT， 我们就会的到 \[ \begin{align*} \hat f_k =& \sum_{j=0}^{n-1}f_j(\omega_n^k)^j\\\\ =& \sum_{j=0}^{n/2-1}f_{2j}(\omega_n^k)^{2j} +\omega_n^k\sum_{j=0}^{n/2-1}f_{2j+1}(\omega_n^k)^{2j}\\\\ =& \sum_{j=0}^{n/2-1}f_{2j}(\omega_{n/2}^k)^{j} +\omega_n^k\sum_{j=0}^{n/2-1}f_{2j+1}(\omega_{n/2}^k)^{j}\\\\ =& E + \omega_n^k O \end{align*} \]

\(E\) 和 \(O\) 便是我们的两个相同的子问题.

如果我们分别利用 \(f_{even}\) 和 \(f_{odd}\) 计算出了\(E\) 和 \(O\) 我们便能很快的根据上面的式子计算出所有的\(\hat f_k\)

不仅如此，我们可以利用 \[ \omega_n^{k+n/2} = -\omega_n^{k} \] 进一步的简化. \[ \begin{align*} \hat f_k=E+\omega_n^kO\\\\ \hat f_{k+n/2}=E-\omega_n^kO \end{align*} \] 这样只需要一半的循环，就能计算出全部的\(\hat f_k\) \[ \begin{align*} \hat f=\left [ \begin{matrix} \hat f_0 \\\\ \hat f_1 \\\\ \hat f_2 \\\\ \vdots \\\\ \hat f_{n-1} \end{matrix} \right ] =& \mathcal{F_n} \left[ \begin{matrix} f_0 \\\\ f_1 \\\\ f_2 \\\\ \vdots \\\\ f_{n-1} \end{matrix}\right]\\\\ =&\left[ \begin{matrix} E_{n/2} & D_{n/2} \\\\\\\\ E_{n/2} & -D_{n/2} \end{matrix}\right] \left[ \begin{matrix} \mathcal{F_{n/2}} & O \\\\\\\\ O & \mathcal{F_{n/2}} \end{matrix}\right] \left[ \begin{matrix} f_{even}\\\\\\\\f_{odd} \end{matrix}\right] \end{align*} \]

如此算法的复杂度达到了\(O(nlogn)\) ，如果利用并行计算，可以达到更快.

FFT的应用

• 像前面做过的去除噪声便是一个很常见的应用
• 图象压缩
• 多项式乘法或者大整数乘法