(快速)傅里叶变换
Published:
笔者在算法与数据结构、图像处理等领域都多次接触到傅里叶变换和快速傅里叶变换，但却一直未能窥其全貌，决定趁此机会好好整理。
由于Fourier变换常常有不同常数的差异，因此可能在某些方面稍有差异，但只要忽略其常数，仅关注其本质，足以窥得其要点。
Descrete Fourier Transform(DFT)
我们考虑从周期函数的Fourier级数展开出发，推导离散Fourier变换（DFT）
对于周期函数，可以用基$1,\cos x,\sin x, \cos 2x, \sin 2x,…$ 将其展开为Fourier级数
使用Euler公式，我们知道$\text{span} (\sin nx, \cos nx) = \text{span}(e^{inx}, e^{-inx})$
因此，可以选择$e^{inx}$作为基展开，得到$f(x) = \sum_{n = -\infty}^{\infty} c_ne^{inx}$
由于这组及$e^{inx}$满足正交性，也即$\langle e^{inx} , e^{imx} \rangle = \int_0^{2 \pi } e^{i(n-m)x} = 2 \pi I(n =m)$ ,其中$I$表示示性函数
根据基的正交性，$c_n = \frac{1}{2 \pi} \int_0^{2 \pi} f(x) e^{-inx}$
如果使用数值积分近似，则可以得到离散Fourier变换的公式，
\[F_n = \frac{1}{N} \sum_k f_k e^{-\frac{i 2 \pi kn}{N}}\]使用$N$次单位根$w_N = e^{-\frac{i2 \pi }{N}}$构成的矩阵$F$，可以将离散Fourier变换表达为矩阵形式，
\[F_n = DFT(x) = F x\]容易根据单位根的性质，$w_n^n = (1-w_n) (1+w_n +w_n^2 +…+ w_n^{n-1}) = 1$验证，$F^{\star} F = N I_N$。
因此可以定义离散逆Fourier变换（IDFT）为：
\[\begin{align} f_n = \sum_k F_k e^{\frac{i 2 \pi}{N} kn} \\ x = \frac{1}{N} F^{\star} \hat x \end{align}\]对于非周期函数，令周期$T \rightarrow \infty$, 则可以得到连续形式的Fourier变换(FT),其可以适用于任意函数，
\[\begin{align*} F(u) = \int_{-\infty}^{\infty} x(t) e^{-itu} dt \end{align*}\]详细的推导需要用到复变函数的知识，此处可以先作为一个定义理解。
关于FT和DFT的关系，如何从DFT推导出FT的公式，后文将展开分析。
Fitting and Approximation
利用Fourier级数或者Fourier变换，也可以进行函数拟合或者插值。
在插值 和 拟合 中，笔者详细讨论了使用多项式进行函数插值和拟合的手段，但并没有讨论基于Fourier变换的三角拟合或三角插值。
Approximation
由于Fourier级数中选取了无限个正交基，考虑在有限个正交基下$\text{span}(e^{-i Nx} , …, e^{iNx})$的最小均方逼近，得到带限的Fourier级数(Truncated Fourier Series):
\[p(x) = \frac{1}{2 \pi } \sum_{k = -N}^{N} e^{ikx} \int_0^{2 \pi} e^{ -ik t} f(t) dt\]求得的多项式$p(x)$也即要求的拟合函数。
Fitting
如果在上式中考虑用离散的数值积分近似连续的积分，得到，
\[\begin{align} p(x) &= \frac{1}{2 N +1 } \sum_{k = -N}^{N} \sum_{j=0}^{2N} e^{ ik x - \frac{ik 2j \pi }{2N+1}} f(\frac{k2j \pi }{2N+1}) dt \\ &= \frac{1}{2 N +1 } \sum_{k = -N}^{N} \sum_{j=0}^{2N} e^{ ik x - ik x_j} f(x_j) dt \text{ (Let } x_j = \frac{k2j \pi }{2N+1})) \\ \end{align}\]上式其实相当于同时做了一个DFT和IDFT，因此如果将插值结点$x_k$代入，可以得到$p(x_k) = f(x_k)$，此时则得到了相应的插值多项式，该插值方法也称为三角插值。
Fast Fourier Transform(FFT)
快速Fourier变换(FFT)实现在$O(n \log n)$的时间内实现DFT。
观察DFT的矩阵表达式，
\[DFT(x) = F x\]本质上相当于根据多项式$f(x) = \sum_i a_i x^i$ 在单位根$w_n$这些点上的取值$f(w_n)$求解系数$a_i$即为DFT的结果。
而单位根具有折半定理，也即，$w_n^2 = w_{n / 2}$，将多项式分为奇偶项展开，
\[\begin{align} f(w_n) &= \sum_i^n a_i x^i \\ &= \sum_i^{n /2} a_{2i} w_n^{2i} +\sum_i^{n /2} a_{2i+1} w_n^{2i+1} \\ &= \sum_i^{n /2} a_{2i} w_{n/2}^{i} + w_n \sum_i^{n /2} a_{2i+1} w_{n /2}^{i} \end{align}\]也即相当于求解在$w_{n /2}$这$n/2$个结点上的系数，然后在$O(n)$时间内对奇偶项的系数进行排列得到最终的系数$a_n$.
时间复杂度可以用递归式表示为，
\[T(n) = T(n/2) + \Theta(n)\]根据算法中著名的主定理，可以知道$T(n) = \Theta(n \log n)$.
Convolution
卷积定理刻画了Fourier变换核心性质，可以说是Fourier变换中最重要的定理。
Continuous Convolution
卷积在信号处理、概率论等中各个领域都是极为重要的基础操作，其定义为 \(\begin{align*} (u \ast v)(s) = \int_{-\infty}^{\infty} u(t) v(s- t) dt \end{align*}\)
对于FT，可以证明其满足如下的卷积定理：
\[FT(u \ast v) = FT(u) \cdot FT(v)\]证明利用积分的交换性即可：
\[\begin{align} FT(u \ast v) &= \int e^{-i sx} dx\int u(y) v(x-y) dy \\ &= \int u(y) dy\int e^{-i sx} v(x-y) dx \\ &= \int u(y)e^{-i sy} dy\int e^{-i s(x-y)} v(x-y) dx \\ &= \int u(y)e^{-i sy} dy\int e^{-i sz} v(z) dz \text{ (Let } z = x-y )\\ &= FT(u) \cdot FT(v) \end{align}\]对于IDFT，也可以有类似的性质，但根据不同FT的定义，其可能用乘以一个常数。
根据卷积定理，可以得到Fourier变换几乎不改变能量，首先由 \(\begin{align*} FT(u \cdot \bar u) \Big{\vert}_{s=0} = \int \vert u(t) \vert^2 \text{d} t \\ \end{align*}\)
再根据卷积定理， \(\begin{align*} FT(u \cdot \bar u) \Big{\vert}_{s=0} &= \frac{1}{2 \pi}\text{FT}(u) \ast \text{FT}(\bar u) \Big{\vert}_{s=0}\\ &= \frac{1}{2 \pi} F(s) \ast \bar F(-s) \Big{\vert}_{s=0} \\ &= \frac{1}{2 \pi} \int \vert F(u) \vert^2 \text{d}u \end{align*}\)
因此我们得到了，
\[\begin{align*} \int \vert u(t) \vert^2 \text{d} t = \frac{1}{2 \pi} \int \vert F(u) \vert^2 \text{d}u . \end{align*}\]在离散情况下，也即对于DFT，上述结论是显然的，因为$F F^{\ast} = N I_N$.
因此，$\Vert DFT(u) \Vert_2 = \Vert Fu \Vert_2 = \sqrt{N} \Vert u \Vert_2$
Circular Descrete Convolution
对于DFT，如果定义循环卷积，$u \ast v = \sum_j u_j v_{j-k}$，相当于假设$u,v$均为周期函数。
容易验证，在上述定义下，DFT也满足卷积定理，证明可以类似FT卷积定理的方法。
\[DFT(u \ast v) = DFT(u) \cdot DFT(v)\]利用循环矩阵，也可以根据矩阵方法证明循环卷积下DFT的卷积定理，可以参见 知乎‘离散傅里叶变换、循环卷积与卷积定理
General Descrete Convolution
更常用的是离散卷积，其定义为，$u \ast v = \sum_{k \le j} u_j v_{j-k}$
此时可以处理$u,v$长度不相等的情况，该方法定义的卷积在图像处理、卷积神经网络、计算数学等都有重要的应用。
该方法定义的卷积，只要对$u,v$填充足够多的零，则可以转化为循环卷积的形式，因此仍然可以利用循环卷积的卷积定理。
考虑两个多项式$p(x),q(x)$的乘法，设其系数分别为$a_n,b_n$，而结果$p(x) + q(x)$的系数为$c_n$, 此时$c_n = \sum_{k \le j} a_j v_{j-k}$
将系数$c_n$表示为卷积的形式后，根据卷积定理可以先用FFT转化为两个函数的逐点相乘，再使用逆FFT得到结果。
因此，使用FFT也可以将卷积操作的复杂度降低为$\Theta(n \log n)$,在实际中应用广泛。
Derivative
由于差分算子可以看作一个卷积核，因此Fourier变换于导数密切相关，本节探究其关系，参考了 知乎’函数导数的傅立叶变换
根据IDFT的公式，
\(\begin{align} FT(x) &= X \\ x(t) &= \frac{1}{2 \pi} \int X(u) e^{i ts} dt \\ \end{align}\) 求两边同时求导可以得到，
\[\begin{align} \frac{dx(t)}{dt} &= iu \frac{1}{2 \pi} \int X(u) e^{i tu} dt \\ FT(x') &= iu FT(x) \end{align}\]另一种方法是根据分部积分公式，
\[\begin{align} FT(x') &= \int x'(t) e^{-iut} dt\\ &=\int e^{-iut} dx(t) \\ &= e^{-iut} dx(t) \Big{\vert}_{-\infty}^{\infty} - \int x(t) d e^{-iut}\\ &= iu \int x(t) e^{-iut} d t\\ &= iu FT(x) \end{align}\]该公式对于图像处理中Laplace锐化等操作都至关重要，也可以应用在Posisson方程的求解等。
From FT to DFT
从连续Fourier变换推到到离散Fourier变换的关键是取样，取样通常基于冲激函数$\delta(t)$，其满足
\[\begin{align} \delta(t)&=0,t \ne 0 \\ \int_{-\infty}^{\infty} \delta(t) dt &= 1 \end{align}\]本质上，冲激函数是广义函数的一种，根据定义，冲激函数满足取样性质，也即使用$\delta(t)$可以对函数进行采样，
\[\int_{-\infty}^{\infty} f(t)\delta(t-t_0) dt = f(t_0)\]考虑常用的等距采样，可以用冲激串$s_{\Delta t}(t)$进行, 表示以$\Delta t$为间隔等距采样，
\[\begin{align} s_{\Delta t}(t) &= \sum_{n = -{\infty}}^{\infty} \delta(t - n \Delta t) \\ \int_{\infty}^{\infty} f(t) s_{\Delta t}(t) &= \sum_{n = -{\infty}}^{\infty} f_n \end{align}\]由于冲激串很重要，我们首先推导其Fourier变换，
由于冲激串为周期函数，将其用Fourier级数展开，
\[s_{\Delta t}(t) = \sum_{n=-\infty}^{\infty} c_n e^{\frac{i 2 \pi n t}{\Delta t}}\]上式中的展开系数由下面公式得到，
\[c_n = \frac{1}{\Delta t} \int s_{\Delta t}(t) e^{\frac{-i 2 \pi n t}{\Delta t}} =\frac{1}{\Delta t} \int \delta(t) e^{\frac{-i 2 \pi n t}{\Delta t}} =\frac{1}{\Delta t}\]下面对于Fourier级数展开后的$s_{\Delta t}(t)$,
\[\begin{align} FT(s_{\Delta t}(t)) &= FT(\frac{1}{\Delta t}\sum_{n= -\infty}^{\infty }e^{\frac{i 2 \pi n t}{\Delta t}})\\ &= \frac{1}{\Delta t}\sum_{n= -\infty}^{\infty } FT(e^{\frac{i 2 \pi n t}{\Delta t}})\\ \end{align}\]由于我们知道，
\[\begin{align} FT^{-1}(\delta(u - u_0)) &= \int \delta(u - u_0)e^{itu} du = e^{itu_0} \end{align}\]令$u_0 = {\frac{2 \pi n }{\Delta t}}$ ,可以得到
\[\begin{align} FT(s_{\Delta t}(t)) &= \frac{1}{\Delta t}\sum_{n = \infty}^{\infty} \delta(t- \frac{2 \pi n }{\Delta t})\\ &= \frac{1}{\Delta t} s_{\Delta t'}(t) \text{ (With} \Delta t' = \frac{2 \pi n }{\Delta t}) \end{align}\]因此，冲激串的Fourier变换仍为冲激串。
使用冲激串$s_{\Delta t}(t)$对函数$f(t)$进行采样得到，
\[\tilde f(t) = f(t) s_{\Delta t} (t)\]根据卷积定理，令$FT(f(t)) = F(u)$,
\[\begin{align} FT(\tilde f(t)) &= FT(f(t) s_{\Delta t} (t)) \\ &= \frac{1}{2 \pi} FT(f(t)) \ast FT( s_{\Delta t} (t)) \\ &= \frac{1}{ 2 \pi \Delta t} F(u) \ast \frac{1}{\Delta t}\sum_{n = \infty}^{\infty} \delta(u- \frac{2 \pi n }{\Delta t}) \\ &= \frac{1}{ 2 \pi \Delta t} \sum_{n= -\infty}^{\infty} \int F(t) \delta(u-t - \frac{2 \pi n }{\Delta t}) dt \\ &= \frac{1}{2 \pi \Delta t} \sum_{n= -\infty}^{\infty}F(u - \frac{2 \pi n }{\Delta t}) \end{align}\]由此，采样后的函数$\tilde f(t)$的Fourier变换是$F(u)$的一个拷贝的无限周期序列。
定义$\tilde F(u)=FT(\tilde f(t))$,
\[\begin{align} \tilde F(u) &= \int_{-\infty}^{\infty} \tilde f(t) e^{-i ut} dt \\ &= \int_{-\infty}^{\infty} \sum_{n = -\infty}^{\infty } f(t) \delta(t - n \Delta t) e^{-i ut} dt \\ &= \sum_{n = -\infty}^{\infty } \int_{-\infty}^{\infty} f(t) \delta(t - n \Delta t) e^{-i ut} dt \\ &= \sum_{n = -\infty}^{\infty} f_n e^{-i un\Delta t} \end{align}\]此时$\tilde F(u)$为连续函数，由于在DFT中我们需要离散的数据点，考虑对$\tilde F(u)$ 在 $1/\Delta T$ 内进行等距采样得到$F_m$, 可以得到DFT的公式，
\[F_m = \sum_{n=0}^{N-1} \exp \left({-\frac{i 2 \pi mn}{N} }\right)\]