两年做过fft海面，当时精力放在代码实现上了，并没有细想背后原理。这些天回顾了一下，把之前得过且过的地方想明白了，当然，仍有盲区，但整个逻辑链条是清晰的了。本系列将完整梳理一下。

当初调试时的视频：

00:15

00:13

本篇先讲基础理论和模型。

一，频谱，傅里叶变换与逆变换

自变数为x的函数f(x)可看作一个随空间变化的信号（空域信号），只要满足一定条件，它就可以表示为一堆不同频率的随空间变化的正弦信号（空域正弦信号）的线性组合（或积分）。

注：是求和还是积分，取决于信号是否具有周期性。

这些正弦信号（称为基底）的频率构成频域。如果将频域作为定义域，相应频率基底的振幅和相位作为函数值，则得到一个新函数F(ω)。称为信号f(x)的频谱。

注：频域（Frequency Domain），空域（Space Domain），时域（Time Domain）。其中空域和时域在数学上是一个意思，不同仅在于自变数用x还是t。对于本文应用场景而言，用空域说法更恰当。

知道了频谱F(ω)也就知道了f(x)，反之亦然，二者是等价的，是同一个信号的两种不同表示方法。

由f(x)求频谱F(ω)，称为傅里叶变换。

由频谱F(ω)求f(x)，称为傅里叶逆变换。

实际傅里叶变换/逆变换不是以正弦信号为基底，而是以复指数信号 为基底，频谱也相应地变成了复数。

二，从三角形式到复数形式

本节对文章目标意义不大，可略过，我只是因为看到傅里叶变换带复数比较懵逼，所以尝试捋了一下。

信号分解为基信号的线性组合（即傅里叶逆变换）有无数种方法，取决于基底形式的选择。不同基底形式导致不同的频谱形式。

（一）以带相位正弦信号作为基底：

此时频谱需要给出振幅和相位，形式为 。

（二）以不带相位余弦和正弦信号联合作为基底：

使用三角恒等式 对（一）作变形，得：

注：可以看出， 取任意有限值都不影响结果。

此时频谱需要给出两个振幅A1和A2，形式为 。

（三）以复指数信号作为基底：

使用欧拉恒等式 对（二）作变形，得：

上式还可以化简：

考虑到 ，也就是说频率只使用了正半轴，函数 和 也仅在正半轴上有定义，负半轴成了三不管地带我们可以为所欲为随便定义，所谓延拓。

我们将 延拓成偶函数，将 延拓成奇函数。

注意：奇函数必过原点，所以 需满足 ，因为前面已经说了「 取任意有限值都不影响结果」，所以这是能做到的。

延拓之后，有：

（由于 )

（定义新函数 ）

（令 ）

化简完成。

需要注意，与前面两种形式不同，此时 的取值范围（即频域）已经拓展到负半轴。

此时的频谱只需给出B(ω)，故频谱形式为F(ω)=B(ω)。

也就是说此时B(ω)恰好就是频谱。那不如干脆写成：

此即复数形式的傅里叶逆变换，是不是很简洁。

由于 是复数，所以频谱F(ω)也是复数。

这可以理解，因为频谱需要能解码出振幅和相位，所以必定需要两个分量。

三，傅里叶逆变换与海面模拟是如何扯上关系的

如果将开放海域的波浪高度看作定义在XZ平面上的空域信号 ，根据经验，这个信号天然地就很接近大量正弦信号的叠加。如果我们知道了其频谱 ，则使用傅里叶逆变换，就可以求出 ，即得到海面的高度场。

由于计算机不能处理连续事物，所以海面模拟里用的傅里叶逆变换是离散傅里叶逆变换（Inverse Discrete Fourier Transform（IDFT)）。

三，海面的IDFT模型

海面高度由下面IDFT给出：

其中：

为空域坐标。

为频域坐标。kx,kz均为频率，相当于前文中的 和 。

为频谱。

这里频谱 较前文多了个参数t，表示此频谱会随时间变化，相应地高度函数就也变成随时间变化的了，所以也加参数t。

另外注意 中与是点乘，即 ，表示：固定z，只让x变化时频率为 ；固定x，只让z变化时频率为 。

求和是对所有频域坐标点 进行。

在频率平面上以原点为中心每隔 （L为海面patch尺寸）取一个点，共NxN个点：

在xz平面上以原点为中心每隔 取一个点，也是共NxN个点：

下图红点为频域点，蓝点为空域点：

对于游戏而言，通常取N=64就够用了，即64^2=4096个频率点，亦即4096个不同频率正弦信号叠加，这个叠加数量比一般gerstner wave多多了，所以细节更丰富。

另外可以验证，在N为偶数的情况下，以上所有频率值k对应的周期T= 均为L的约数，所以 横向和纵向均以L为周期，也就是说海面patch是可以无缝tiling的。

四，菲利普频谱（Phillips spectrum）

那么频谱为何呢？

通常采用的是菲利普频谱，公式如下：

其中， 表示 的共轭复数，k表示 的模。

其中g为重力加速度。

其中 和 是相互独立的随机数，均服从均值为0，标准差为1的正态分布。

其中，w是风向。L和前面那个海面patch尺寸L字母用重了，不是一个东西，L定义见下一条。

其中V为风速。

以上公式源于海洋统计学，我暂时无法给出更多解释。

其中生成服从标准正态分布的随机数对儿 ，方法见：杨超：生成服从标准正态分布的随机数

五，法线

上面用IDFT得到海面高度以后，用差分就可以求法线，但那样求出来的不是很精确。最精确的方法是直接推出法线的解析式。

因为高度是：

其空间梯度为：

（因为 并不包含空间变元，故可从梯度符号内移出）

而其中

所以

又 梯度向量、up向量、法向量 三者具有下图所示几何关系：

所以有（注意，normalize是必要的）：

六，尖浪（Choppy Waves）

在gerstner wave中，为了表现波峰尖角，使用下面公式在xz平面内进行挤压（红框中部分）：

此处IDFT海面，同样需要类似挤压操作，公式为：

不难看出，两者虽然写法不同，含义是一样的：即对sin波进行cos挤压，对cos波进行sin挤压。

当xz平面内挤压过头时，就会出现刺穿（如图所示）。恰好对应浪尖破碎形成泡沫的区域。

当发生刺穿时，局部发生翻转，表现在数学上，即面元有向面积变为负值。

那么如何求面元有向面积呢？同济高数下册里学过：二重积分换元法，雅可比行列式。

（找不到高数书的也可看这个：杨超：二重积分换元法、雅可比行列式）

因为x和z均为以x,z为变元的二元函数，即x=x(x,z),z=z(x,z)，由二重积分换元法知变换后面积元为：

由于dxdz必定为正数，所以dA的正负就取决于雅可比行列式

的正负。

由于 ，所以：

由于我们有 的表达式，所以其实上面各偏导数都是可以求出来的。

如果真的去求，我们会发现 ，亦即 ，验证如下：

七，快速傅里叶变换（FFT）

IDFT如果用暴力求和的方法来计算会卡成翔，需要使用FFT来搞。

FFT演算法本来也算是基础理论，但本文已经比较长，而且今天太晚写不完了，就放到下一篇吧。

下一篇：杨超：fft海面模拟（二）

推荐阅读：

相关文章