游戏开发中的数学：傅里叶变换(FT)篇

洗尽铅华

零、摘要与目录
之前在大学时学的是电子信息科学与技术专业，同时也参加过 ACM 竞赛，因此对傅里叶变换（后都用 FT 替代）就已有一部门了解，而如今转眼间已经在游戏行业打工三年，没想到 FT 在游戏开发中依旧有应用，可谓是贯穿始终了，但是无论是电科专业中的 FT、还是 ACM 竞赛算法中的 FT、又或者是游戏开发中的 FT，三者学习与了解时的侧重点、思路、应用目的都有所分歧（当然从需要掌握的深度与难度来讲，要求最高的还是电子信息技术）
因此，时隔4年半筹备专门再写一篇文章：此次会站在一个游戏开发者，准确来讲是图形衬着入门开发者角度，再次对 FT 做一个简要的介绍与解析，和 ACM 的目的是 AC 标题问题一样，最终也会展示一个简单的应用的例子，毕竟抛开需求谈技术就是耍地痞，一切目的都是为了实现
文章难度：科普及入门，和类于 Heinrich：傅里叶分析之掐死教程（完整版）更新于2014.06.06 分歧，这篇文章还是会涉及一些公式推导的，当然都很简单，只需要对积分和三角函数有最简单的了解即可
如有错误，欢迎指出，万分感激
0.1 总目录

由于篇幅较长，总计会拆分成四篇文章

• 一、Games101 中呈现的傅里叶变换(FT)的简单推导
  
  
  
  • 1.1 前置1：三角函数系
    
  • 1.2 前置2：三角函数的指数暗示
    
  • 1.3 傅里叶级数的推导
    
    
    
    • 1.3.1 任意周期函数的傅里叶级数
      
    • 1.3.2 傅里叶级数的复数形式
      
    
    
  • 1.4 非周期函数傅里叶变换推导
    
    
    
    • 1.4.1 时域与频域
      
    • 1.4.2 傅里叶变换的最终形式
      
    • 1.4.3 一个形象的解释，为什么梳状函数的频域图像仍然是个梳状函数
      
    
    
  
  
• 二、“虚拟世界”中的傅里叶变换：DFT
  
  
  
  • 2.1 采样与冲激函数(Dirac δ 函数)
    
  • 2.2 时域与频域的离散化计算
    
  
  
• 三、快速傅里叶变换(FFT)
  
  
  
  • 3.1 预备性质与数学道理
    
    
    
    • 3.1.1 单元复 N 次方根
      
    • 3.1.2  FFT 蝶形算法
      
    
    
  • 3.2 速解多项式乘法：从另一个角度看 FFT
    
    
    
    • 3.2.1 何为卷积
      
    • 3.2.2 nlogn 的卷积计算
      
    
    
  • 3.3 FFT 的简单法式实现
    
  
  
• 四、FT 在图像措置中的应用：二维傅里叶变换
  
  
  
  • 4.1 平面正弦波与二维频率域(K-Space)
    
  • 4.2 低通滤波的含义
    
  • 4.3 卷积核(kernel)与滤波算子
    
  • 4.4 纹理采样与抗锯齿
    
    
    
    • 4.4.1 频域混叠现象
      
    • 4.4.2 奈奎斯特频率(Nyquist frequency)
      
    
    
  • 4.5 附录1：滤波技术在 ShadowMapping 软暗影中的应用
    
  • 4.6 附录2：高斯模糊的本质
    
  • 4.7 附录3：球谐光照
    
  
  
• 五、FT 在图形衬着中的应用：基于 FFT 的波浪模拟
  
  
  
  • 5.1 FFT 海洋公式：二维 IDFT
    
    
    
    • 5.1.1 二维 IFFT
      
    
    
  • 5.2 波浪统计模型
    
    
    
    • 5.2.1 [2004] 菲利浦频谱(Tenssendorf Phillips spectrum)
      
    • 5.2.2 [2015] TMA JONSWAP 频谱与 Hasselman 标的目的扩展
      
    • 5.2.3 海洋初始状态
      
    • 5.2.4 海洋公式最终形式
      
    
    
  • 5.3 基于 GPGPU 的 Ocean IFFT 计算
    
    
    
    • 5.3.1 初始频谱计算
      
    • 5.3.2 任意时刻频谱计算
      
    • 5.3.3 海洋法线
      
    • 5.3.4 IFFT 蝶形图生成
      
    • 5.3.5 海洋高度场、梯度与程度偏移量
      
    
    
  • 5.4 Extra：基于 LOD 的无限大海洋的频谱填充
    
  • 5.5 Extra：一种移动平台的 FFT 可行方案：离线 FFT
    
  
  

• 第一篇：Games101 中呈现的傅里叶变换(FT)的简单推导
  
• 第二篇：离散傅里叶变换(DFT)
  
• 第三篇：FT 在图像措置中的应用
  
• 第四篇：FT 在图形衬着中的应用：基于 FFT 的波浪模拟
  

<hr/>一、Games101 中呈现的傅里叶变换(FT)的简单推导

到底什么是傅里叶变换：它的物理意义是什么，公式又从何而来？

以下的内容呈此刻 Games101 中的第六章：光栅化（深度测试与抗锯齿）


此中，课程中这一部门大致介绍的内容是
任何周期函数 f(x) 都可以分化为无数个分歧频率、分歧幅值的正、余弦信号。和泰勒展开一样，这也是一种函数展开/逼近方案：即傅里叶级数展开，此中左边那一张 PPT 就是一个简单的例子
右侧则为傅里叶变换及其逆变换：它和傅里叶级数展开的概念长短常相似的
如果不了解傅里叶变换(FT)，可能会对视频中右边突如其来的内容非常陌生，只知道它用于将一个函数变成另一个函数，但是这个式子是怎么推导的，具体是什么变成什么可能就没有太多头绪
而这一章的内容，主要就是对上面内容的答疑解惑

1.1 前置1：三角函数系

三角函数系即包含所有分歧周期的三角函数的调集
\{0,1, \sin x, \cos x, \sin 2 x, \cos 2 x, \cdots, \sin n x, \cos n x\}
此中三角函数拥有一个非常重要的性质：正交性

• \int_{-\pi}^\pi \cos n x \cos m x d x=0 \quad n \neq m
  
• \int_{-\pi}^\pi \sin n x \sin m x d x=0 \quad n \neq m
  
• \int_{-\pi}^\pi \sin n x \cos m x d x=0
  

和向量垂直内积为0一样，两个函数内积 \int_{x_0}^{x_1} f(x) g(x) d x=0 即两个函数正交，而对于上面的两个特例：即 m = n 的情况，其积分成果为 \pi
正是依赖这个特性，后续才可以任意函数中分手出指定频率的波形

1.2 前置2：三角函数的指数暗示

3Blue1Brown 有一个视频很好的介绍了 e^i  到底表达了什么，如果有兴趣的话可以直接去看：【官方双语】微分方程概论-第五章：在3.14分钟内理解e^iπ，搞懂这个之后，这一部门就可以直接跳过了

虚数 i 这个概念基本上只要上过高中的都知道，只不外大大都人对它的印象只剩下了一个公式 i^2 ：其余的就都不记得了，那么乘以虚数的意义又是什么呢？


如上图有三个线段线段，线段 AB 的长度是1：当它乘以2的时候，它的长度发生了变化，拉伸成了线段 EF，而当它乘以 -1 的时候，就变成了线段 CD，或者说代表着原线段在数轴上围绕原点旋转了 180°
既然乘以 -1 相当于旋转 180°，考虑到 i^2，那么是不是就可以理解为：乘以 i 就是旋转了 90°


是的，不外到此如果只考虑一维的实数轴就不够用了，因此我们可以加一条垂直的虚数轴已构成一个复平面，这样我们就可以把虚数乘法也暗示出来：即乘以实数表拉伸，乘以虚数表旋转
但是到这里还没有结束，要知道本章之所以要提到虚数 i，主要是为了解释后面呈现的 e^i ，以及经典的欧拉公式 e^{i x}=\cos x+i \sin x ，可能就算你理解了虚数乘法，也未必能理解虚数作为指数意味着什么
一个事实是：其实很多时候某些东西单独看是没有意义的，只不外为了实现一些目的，我们后来的人们手动赋予了它意义，尽管你完全理解不了：是的它就是一个定义，数学家就是这么理解它的（当然其必然是有依据的，而不是乱来）
而虚数作为指数就是这么被“定义”出来的，其实它已经和正常我们理解的指数完全不是一个概念了，也就是说我们不能说 e^i 就是持续 i 个 e 相乘，它是错误且没有意义的，事实上，它表达的是下面这个东西：


即 e^i 暗示为单元线段旋转必然角度，这个角度对应的弧长为1，这个定义是按照复导数 \frac{d}{d t} e^{i t}=i \cdot e^{i t} 反推出来的，当然具体细节就不再展开，有兴趣的话可以本节看上面引用的视频
知道这个后，两个著名的公式其实就不难理解了，你也可以瞬间化为数学家去把它们推出来：

• e^{i \pi}+1=0
  
• e^{i x}=\cos x+i \sin x
  

后面主要会用到的，其实就是的下面的公式

1.3 傅里叶级数的推导

好了，可以回来了前面的内容：任何周期函数 f(x) 都可以分化为无数个分歧频率、分歧幅值的正、余弦信号，不妨可以先考虑最简单的情况，也就是先保证 f(x) 周期 T = 2\pi ，即 f(x)=f(x+2 \pi)
那么对于上面部门的数学暗示，就有：f(x)=\sum_{n=0}^{\infty} a_n \cos n x+\sum_{n=0}^{\infty} b_n \sin n x，此中 a_n 和 b_n 都为常数项，然后有意思的就来了：有了这个级数暗示，我们能不能得到一个左侧为 a_n 的等式呢？
当然是可以的，就操作上面的三角函数的正交性求积分即可
先计算最简单的 a_0 ：
f(x)=\sum_{n=0}^{\infty} a_n \cos n x+\sum_{n=0}^{\infty} b_n \sin n x = \\ a_0+\sum_{n=1}^{\infty} a_n \cos n x+\sum_{n=1}^{\infty} b_n \sin n x
两边积分有：
\int_{-\pi}^\pi f(x) d x=\int_{-\pi}^\pi a_0 d x+\int_{-\pi}^\pi \sum_{n=1}^{\infty} a_n \cos n x d x+\int_{-\pi}^\pi \sum_{n=1}^{\infty} b_n \sin n x d x
按照三角函数的正交性，可以知道后面两项看上去复杂，其实都为0，也就有：
\int_{-\pi}^\pi f(x) d x=a_0 \int_{-\pi}^\pi d x=\left.a_0 x\right|_{-\pi} ^\pi=2 \pi a_0
好了，到此成果显而易见了
接下来用同样的方式计算 a_n：不外这一次我们在两边同时积分之前，先全部乘上 cosmx
\int_{-\pi}^\pi f(x) \cos m x d x= \\ \int_{-\pi}^\pi \frac{a_0}{2} \cos m x d x+\int_{-\pi}^\pi \sum_{n=1}^{\infty} a_n \cos n x \cos m x d x+\int_{-\pi }^{\pi} \sum_{n=1}^{\infty} b_n \sin n x \cos mx d x \\一样按照三角函数的正交性，第一项和第三项为0：就有
\int_{-\pi}^\pi f(x) \cos m x d x=\int_{-\pi}^\pi \sum_{n=1}^{\infty} a_n \cos n x \cos m x d x
而且当前仅当 n = m 时，右侧的积分成果不为0，也就是说
\int_{-\pi}^\pi f(n) \cos n x d x=\int_{-\pi}^\pi \operatorname{a_ncos} n x \cos n x d x=\operatorname{a_n} \int_{-\pi}^\pi \cos ^2 nx d x
a_n=\frac{1}{\pi} \int_{-\pi}^\pi f(x) \cos n x d x
游戏结束，关于 b_n 的计算同理
好了，到此就可以得到一个最基础的傅里叶级数形式：
对于满足周期 T = 2\pi 的 f(x) ， f(x)=f(x+2 \pi) ，有：
f(x)=\frac{a_0}{2}+\sum_{n=1}^{\infty}\left(a_n \cos n x+b_n \sin n x\right)
此中：
\begin{aligned} &a_n=\frac{1}{\pi} \int_{-\pi}^\pi f(x) \cos n x d x \\ &b_n=\frac{1}{\pi} \int_{-\pi}^\pi f(x) \sin n x d x \end{aligned}
1.3.1 任意周期函数的傅里叶级数

前面考虑的f(x) 周期满足最简单的 T = 2\pi ，此刻在理解了上面内容的基础值之上，可以考虑去推导任意周期的通式
考虑周期 T = 2L ，L 为任意值，也就是 f(x)=f(x+2L) ，那么就可以用换元法：
f(t)=f\left(\frac{L}{\pi} x\right) \triangleq g(x) ，此中 g(x)=g(x+2 \pi) ，满足前面的傅里叶变换的级数形式
那么此时就只需要把 x=\frac{\pi}{L} t = wt  带入到前面的级数推导中，就可以得出任意周期函数的傅里叶级数形式：
f(t)=\frac{a_0}{2}+\sum_{n=1}^{\infty}\left(a_n \cos n \omega t+b_n \sin n \omega t\right)
此中：
\begin{aligned} &a_n=\frac{2}{T} \int_0^T f(t) \cos n \omega dt \\ &b_n=\frac{2}{T} \int_0^T f(t) \sin n \omega dt \end{aligned}
1.3.2 傅里叶级数的复数形式

好了转折点来了，此刻我们把前面第二节的欧拉公式拿过来： e^{i x}=\cos x+i \sin x
这个依旧是可以逆推得到（因为离本章的主题斗劲远，这一部门的详细推导就略了）
\cos \theta=\frac{1}{2}\left(e^{i \theta}+e^{-i \theta}\right) \\ \sin \theta=-\frac{1}{2} i\left(e^{i \theta}-e^{-i \theta}\right)
将其带入前面的傅里叶级数 f(t)=\frac{a_0}{2}+\sum_{n=1}^{\infty}\left(a_n \cos n \omega t+b_n \sin n \omega t\right) 有
f(t)=\frac{a_0}{2}+\sum_{n=1}^{\infty}\left[\frac{1}{2}a_n\left(e^{i n \omega t}+e^{-i nw t}\right)-\frac{1}{2} ib_n\left(e^{i n \omega t}-e^{-i n \omega t}\right)\right] \\  =\frac{a_0}{2}+\sum_{n=1}^{\infty}\left[\frac{a_n-i b_n}{2} e^{i n \omega t}+\frac{a_n+i b_n}{2} e^{-i n \omega t}\right] \\ =\sum_{n=0}^0 \frac{a_0}{2} e^{i n \omega t}+\sum_{n=1}^{\infty} \frac{a_{n} -ib_n}{2} e^{i n \omega t} + \sum_{n=-\infty}^{-1} \frac{a_{-n} + ib_{-n}}{2} e^{i n \omega t} \\
到这你会发现，其实这三个求和公式对应的区间加在一起正是负无穷到正无穷，而且它们的项都是 e^{i n \omega t} ，只是系数分歧而已，因此可以直接将这三项整合起来得到：
$$f(t)=\sum_{-\infty}^{\infty} c_n e^{i n \omega t}$$
此中：
c_n= \begin{cases}\frac{a_0}{2}, & n=0 \\ \frac{a_n-i b_n}{2}, & n=1,2,3,4 \cdots \\ \frac{a_{-n}+ib_{-n}}{2} & n=-1,-2,-3,-4 \cdots\end{cases}  
好了，到此就只剩下最后一步了：那就是把前面 a_n 和 b_n 的成果代进去，就有
c_n=\frac{a_{n}-ib_{n}}{2}  =\frac{1}{2}\left(\frac{2}{T} \int_0^T f(t) \cos n \omega t-i\frac{2}{T} \int_0^T f(t) \sin n \omega t\right) \\ =\frac{1}{T} \int_0^T f(t)(\cos n \omega t-i \sin n \omega t) d t
仔细看，右边像不像欧拉公式 e^{i x}=\cos x+i \sin x ？此时配合奇函数偶函数的积分性质，就可以得出 c_n =\frac{1}{T} \int_0^T f(t) e^{-i n \omega t} d t
而对于其它情况例如 c_n =\frac{a_{-n}+ib_{-n}}{2} 的推导，你也会得到一个完全不异的成果，可以说长短常的妙，到此为止，我们就得到了一个完美的傅里叶级数复数形式：
对于周期为 T 的函数 f(t)  ，有：

• c_n =\frac{1}{T} \int_0^T f(t) e^{-i n \omega t} d t , w_0 = \frac{2\pi}{T}
  
• f(t)=\sum_{-\infty}^{\infty} c_n e^{i n \omega t}
  

通过等式①，我们可以对一段函数进行滤波，把指定频率的波段/信号给筛出来

1.4 非周期函数傅里叶变换推导

休息一会再来，前面对于 f(t) ，我们求得了他为周期函数情况下的傅里叶级数展开，此刻测验考试求出真正的“通解”：即任意 f(t) 的傅里叶变换
照应开头，主角要登场了！
1.4.1 时域与频域

这里强烈保举一篇文章：傅里叶分析之掐死教程（完整版），里面关于频域与时域的介绍非常形象，而且和本文长篇的公式推导分歧，里面基本不含任何复杂的数学公式，如果看多了 LaTex 换个表情也是可以的

还是阿谁周期为 T 的函数 f(t) ，如果我们画出对应的函数图像，横轴为 t，纵轴为 f(t) ，那么这个图像就是尺度的时域图像，此中横轴变量 t 就代表着时间(Time)


然后就是傅里叶级数：

• c_n =\frac{1}{T} \int_0^T f(t) e^{-i n \omega t} d t , w_0 = \frac{2\pi}{T}
  
• f(t)=\sum_{-\infty}^{\infty} c_n e^{i n \omega t}
  

对于第②个式子，真正能决定最终函数成果的其实是此中的系数 c_n ，尔后面的 e^{i n \omega t}  其实是一种法则，代表着分歧频率的波
因此，有另一个函数图像也可以用于暗示 f(t) ，其横轴为 n \omega ，纵轴为对应的常量 c_n ，这便是 f(t) 的尺度频域图像，该函数也被称为频谱函数，正好是上面的第①个式子，此中横轴 n \omega 对应着频率，而纵轴 c_n 就对应着该频率的波的幅值

如果要更准确的说法，常量系数 c_n 应该是一个复数，因为不止考虑幅值，还应考虑相位，不外我感觉这篇文章作为科普性质的话，讲这个稍微有点超纲了，所以只需要测验考试去理解下就好，后面就默认只考量幅值

对于周期 T 有限的情况下，n \omega 的成果是离散的，因此对于周期函数的频域图像，其更像是一个冲激函数图像，即只可能会在横轴坐标为 \{0,w, 2w, 3w, \cdots, nw\}  的位置下呈现非0的 c_n 成果
将频域图像和时域图像联系起来，就是下面这个样子（该 GIF 来源于维基百科）：


从频域图中也可以看出，周期函数 T 越大，对应的 w 越小，冲激函数的图像就越加密集，而对于非周期函数 T，其实也可以理解为 T 无穷大，此时对应的频率 w 无穷小，频域图像就看上去是持续的了：


同理，如果拉伸 f(t) 的时域图像，那么其频域图像则会被压扁
1.4.2 傅里叶变换的最终形式

基频 w_0 = \frac{2\pi}{T} 相当于周期函数的傅里叶级数中两个相邻频率的差值 (n+1) \omega_0-n \omega_0 ，当周期 T 无穷大时，我们可以把它记作 d\omega 或 \Delta \omega ，这样就得到了针对非周期函数的频谱函数：
c_n=\frac{\Delta \omega}{2 \pi} \int_{-\infty}^{+\infty} f(t) e^{-i w t} d t
此时 c_n  趋近于0，意义不大，但是我们可以将其带入到上面②式 f(t)=\sum_{-\infty}^{\infty} c_n e^{i n \omega t} 就有
f(t)=\sum_{n=-\infty}^{+\infty}\left(\frac{\Delta \omega}{2 \pi} \int_{-\infty}^{+\infty} f(t) e^{-i \omega t} d t\right) e^{in \omega t}= \frac{1}{2 \pi}\int_{-\infty}^{+\infty}\left(\int_{-\infty}^{+\infty} f(t) e^{-i \omega t} d t\right) e^{i \omega t} d \omega
此中，傅里叶变换(FT) 即
F(\omega)=\int_{-\infty}^{+\infty} f(t) e^{-i \omega t} d t ，用频率 f 来描述就是 F(f)=\int_{-\infty}^{+\infty} f(t) e^{-i 2{\pi}f t} d t  
对于 f(t)=\frac{1}{2 \pi} \int_{-\infty}^{+\infty} F(\omega) e^{i \omega t} d \omega ，则为傅里叶变换的逆变换
是的，其实本质上 c_n=\frac{\Delta \omega}{2 \pi} F(\omega) ， F(\omega) 即信号幅值 c_n 在频域中的分布密度
如果文字和公式不好理解，那就上图：


对于 f(t) 的部门图像如上，那么该部门转成复数形式的 f(t) e^{-i \omega t} 图像就为下（ w=\frac{2\pi}{5} ，其取值越小，整体复平面下的图像就越“拥挤”，即这里的 w 对应的是下图中的“转速”（0.2圈/秒），而非函数 f(t) 的频率）


而很明显，我们最后对上面成果求得的积分，就是对这张图像无限延伸的版本中的无数个点对应的向量做加和，想象一下最终加和成果是什么样子的：
由于 w=\frac{2\pi}{5} ，对应转速只有时域周期满足 T=5 的波形，最终得出的傅里叶变换成果不为0（左下图），而且会随着积分上下限的提高，得出的成果也越来越大，其余波形成的缠绕图像城市在复平面上均匀展开，从而得出0或者接近0的成果（右下图），从这也可以看出傅里叶变换滤波的特性


好了，到此就完成了最前面右边那张 PPT 中的内容的推导，是不是还是挺简单的
当然也有可能看懂了推导过程，最后还是有疑问：折腾了这么久，傅立叶变换到底有什么意义和实际应用呢？不着急，到这里还只是一个开始，路还很长很长
1.4.3 一个形象的解释，为什么梳状函数的频域图像仍然是个梳状函数

以下图片同样节选自 Games101 中的第八章：光栅化（深度测试与抗锯齿）


左侧为梳状函数的时域图像
考虑到对于如上的冲激函数序列（即梳状函数）画出其指数空间下的图像


这个图像的性质就是：所有图像中的点城市遍布在单元圆上，而且脉冲信号越密或者采样时间越长，点在圆上的摆列也越密，很都雅出来：颠末无穷时间后，所有的点就会汇集成一个尺度的圆或者在圆上等距摆列，此时做采样积分，最后得出的成果必然是无限接近于0的，但是只有一种情况特殊，那就是转速(s/circle)刚好和脉冲间隔 f 完全一致或者说后者是前者的整数倍：此时所有的点城市落到图中的 (1, 0) 坐标上
因此，梳状函数的频域图像就也会是一段段脉冲信号的反复