FT 在图形渲染中的应用：基于 FFT 的海浪模拟

接上文：FT 在图像处理中的应用

五、一个大型案例：基于 FFT 的海浪模拟

前置：​​​​​

• 傅里叶级数与傅里叶变换
• 离散傅里叶变换(DFT)
• ​​​​​​​​FT 在图像处理中的应用​

5.1 FFT 海洋公式：二维 IDFT

https://tore.tuhh.de/bitstream/11420/1439/1/GPGPU_FFT_Ocean_Simulation.pdf

https://www.youtube.com/watch?v=ClW3fo94KR4

既然任意连续且收敛的函数  都可以分解为无数个不同频率、不同幅值的正、余弦信号，且这个性质可以扩展到二维，即任意一个二维图像都可以分解为无数个复平面波  的叠加

那么海浪作为“波”的典型，是否也可以将其拆成任意方向不同频率强度的基波的叠加？当然没有问题：令  为  时刻，位置  上海面的高度，其构建公式就为

可以看出，这正是二维离散傅里叶变换的逆变换(IDFT)版本，只不过多了一个量度即时间，考虑到 DFT 形式下任意  将对应所有可能的复平面波 ，就能顺理成章的推出波矢量

其中常量  为海平面大小，采样点  满足 ，N 为采样点数量，当然也可以做个变换，即将范围映射到 ，此时

 本质不变

代入上面的海洋方程，并将所有向量全部展开，为方便可以直接任  就有

当然也可以任 ，这在计算方式上是相同的

5.1.1 二维 IFFT

前面已经详细介绍过了一维 FFT 的原理和计算过程，这里同样要扩展二维以得到最终的 ，其实看这个二维 IDFT 式子，你就会发现结论很明显了：

它其实是可以拆成两个一维的 IDFT 的：

• 
• 

那么最终我们计算海洋公式的步骤就应如下：

1. 计算  生成频谱，如果是交给 GPU 计算，可以将结果存储在一张 NxN 的纹理中（关于  的具体内容及计算方式可以参考下一节 5.2 的内容），t 可以理解成第几帧，必然每一帧这个图像都会不同

2. 执行 N 次水平 FFT：，上一步的纹理作为输入，输出结果作为下一步的输入

3. 执行 N 次纵向 FFT：，得到的结果再进行最后的符号计算：即乘上，搞定

到此，我们就拿到了一张高度图，可以用于后续的顶点动画及渲染着色了

5.2 海浪统计模型

事实上，这种基于 FFT 或者说基于频谱逆向推出海洋水面高度的方法，本质上是一种统计模型方法，这里面涉及不少其他领域内的知识，例如海洋物理学，又或是对现实海洋的实时观测统计与实验，我们当然不用深入了解这部分的内容，只要拿对应的统计结果来用就可以了

上面的内容还有一个很重要的东西还未知，就是用于表达海洋波的频率的集合、IDFT 中的频域函数 ，搞定了这个之后，我们才能通过 IDFT 算法将频谱转化为时域的值，从而得到海洋随时间变化的高度图

公式都给你写好了，照着抄就行

5.2.1 菲利浦频谱(Phillips spectrum)

菲利浦频谱来自海洋动力学的成果，是一个反馈了真实海洋的经验模型，也是这里一切的核心，当然这个模型并不唯一

对于菲利浦频谱 ，其中

• ， 是风速， 为风向， 为引力常数 
•  为  的幅值，即频率大小

从论文得到的信息是：菲利浦频谱主要由两个单独部分组成，分别为

1. 非定向波谱 ，由海洋统计分析得出
2. 方向扩展函数 ，体现的是风向对波最终幅值的贡献

5.2.2 海洋初始状态

回到前面的 ，可以先考虑它的初始状态，也就是 t = 0 的时刻的振幅  有：

 和  为两个相互独立的服从均值为0，标准差为1的高斯随机数，实部表初始幅度，虚部表初始相位，而  正是前面的方向波谱

这里再介绍另一个案例，其实后面的 FFT 实现也正是借鉴的这篇视频，这个案例中额外考虑了频谱到波的离散转化：

5.2.3 海洋公式最终形式

最后就是把时间参数加上：

，这里又有

• 角频率与波长的频散关系  满足 ，其考虑了引力风波和重力场的关系，当然如果额外考虑水深 h，其频散关系 
•  是  的共轭复数

到此为止，整个 IDFT 公式就明朗了，然后就是套用上面 IFFT 的步骤

5.3 基于 GPGPU 的 Ocean IFFT 计算

通过上面的内容，可以预见的是生成高度图会有大量的数学计算，尽管这已经是算法加速后的结果了，因此大多数情况下都会把这部分计算任务交给更擅长并行计算的 GPU 而非 CPU

在 Unity 中，往往使用 ComputeShader 来实现

5.3.1 初始频谱计算

先是计算菲利浦频谱  以生成对应纹理，看上去最陌生的公式其实直接抄作业就行了，其中  作为纹理坐标，其余除了风向  都为常量。很明显，游戏中我们不太可能每帧去改变风向，因此风向  也可以视作常量（由美术或策划去配置），在这种情况下 的计算理论可以离线或者只在风向改变时做一次，不需要每帧去重复计算

然后就是高斯随机数  生成，这个只是为了得到不同的初始状态，因此同样可以离线计算：

//如果本地已经存在 GuassTexture，则不需要生成了，直接读取，否则在 Awake 时生成一份
private void Awake()
{
    gaussianNoise = GetNoiseTexture(size);
}

Texture2D GetNoiseTexture(int size)
{
    string filename = "GaussianNoiseTexture" + size.ToString() + "x" + size.ToString();
    Texture2D noise = Resources.Load<Texture2D>("GaussianNoiseTextures/" + filename);
    return noise ? noise : GenerateNoiseTexture(size, true);
}

float NormalRandom()
{
    return Mathf.Cos(2 * Mathf.PI * Random.value) * Mathf.Sqrt(-2 * Mathf.Log(Random.value));
}
    
Texture2D GenerateNoiseTexture(int size, bool saveIntoAssetFile)
{
    Texture2D noise = new Texture2D(size, size, TextureFormat.RGFloat, false, true);
    noise.filterMode = FilterMode.Point;
    for (int i = 0; i < size; i++)
    {
        for (int j = 0; j < size; j++)
        {
            noise.SetPixel(i, j, new Vector4(NormalRandom(), NormalRandom()));
        }
    }
    noise.Apply();

    //此处尝试将 GuassTexture 保存到本地，代码略
    return noise;
}

GPU Gems3:Chapter 37 中提到过如何生成满足服从均值为0，标准差为1的高斯随机数，公式如下：

• 
• 

其中  为均匀分布的随机数

然后两者相组合，就能得到 ：

public void CalculateInitials(WavesSettings wavesSettings, float lengthScale)
{
    //设置计算频谱的各项参数，例如风速、风向、重力加速度常量等
    wavesSettings.SetParametersToShader(initialSpectrumShader, KERNEL_INITIAL_SPECTRUM, paramsBuffer);

    //将计算的 H0k 保存到一张纹理中
    initialSpectrumShader.SetTexture(KERNEL_INITIAL_SPECTRUM, H0K_PROP, buffer);
    //计算角频率与波长的频散关系 w 保存到另一张纹理中
    initialSpectrumShader.SetTexture(KERNEL_INITIAL_SPECTRUM, PRECOMPUTED_DATA_PROP, precomputedData);
    //传入高斯随机数纹理，用于计算 H0K
    initialSpectrumShader.SetTexture(KERNEL_INITIAL_SPECTRUM, NOISE_PROP, gaussianNoise);
    //交给 GPGPU 开始计算
    initialSpectrumShader.Dispatch(KERNEL_INITIAL_SPECTRUM, size / LOCAL_WORK_GROUPS_X, size / LOCAL_WORK_GROUPS_Y, 1);

    //计算得到 H0k 后，把共轭结果 H*0K 也存下来
    initialSpectrumShader.SetTexture(KERNEL_CONJUGATE_SPECTRUM, H0_PROP, initialSpectrum);
    initialSpectrumShader.SetTexture(KERNEL_CONJUGATE_SPECTRUM, H0K_PROP, buffer);
    initialSpectrumShader.Dispatch(KERNEL_CONJUGATE_SPECTRUM, size / LOCAL_WORK_GROUPS_X, size / LOCAL_WORK_GROUPS_Y, 1);
}

中间用到了一个专门用于计算频谱函数的 ComputeShader，里面就是纯计算，并且由于计算模型并不唯一，具体代码就不贴了

5.3.2 任意时刻频谱计算

关于  的代码实现更简单，毕竟前面已经得到了  以及 （该图片来源于 youtube《Ocean waves simulation with Fast Fourier transform》，后同）：

float4 wave = WavesData[id.xy];
float phase = wave.w * Time;
float2 exponent = float2(cos(phase), sin(phase));
float2 h = ComplexMult(H0[id.xy].xy, exponent)
    + ComplexMult(H0[id.xy].zw, float2(exponent.x, -exponent.y));
float2 ih = float2(-h.y, h.x);

不过考虑到海面并不只有高度变化，还有水平方向的流动，后续只有一张 y 轴的高度图还是不够的，还需要 x 轴和 z 轴的偏移图，这里可以先根据  获取水平方向的频谱

• 
• 

Gerstner 波为了展现波峰尖角，也同样对 xz 平面做了变换

前面代码将  一起塞到了一张4通道纹理中，这里继续拿来用：

float2 displacementX = ih * wave.x * wave.y;
float2 displacementY = h;
float2 displacementZ = ih * wave.z * wave.y;

5.3.3 海洋法线

其实到这这一步应该就可以准备 IFFT 计算了，不过中间还有一个东西没有考虑，那就是法线图的生成，因为你要计算光照，法线图是必须的，实时的海浪必然意味着法线图也需要实时计算

这里关于法线计算有两种方案

1. 得到高度图后，可以通过差分方法来求法线，即对变换后的新顶点 ，求出与其相邻的两组顶点   和  ，交叉连线后得到两个向量，叉乘得出法线，

2. 通过偏导计算求出曲面梯度，然后再根据梯度得到法线

方案①逻辑非常简单，可以直接在着色器中实时计算得到法线，缺点也很明显，那就是仅能得到近似法线，这个结果是不会准确的，所以后续若想要得到一个比较好的光照效果，就只能采取方案②

梯度的本意是一个向量，它的方向与取得最大方向导数的方向一致，而它的模为 方向导数的最大值，也可以理解为函数在该点处沿着该方向变化最快，变化率最大

设函数  在平面区域 D 内具有一阶连续偏导数，则对于每一点 ，都可定出一个向量 ，该向量称为函数  在点  的梯度，记作

可以证明：

对于高度图 ，其梯度满足

，其中偏导

• 
• 

float2 displacementY_dx = ih * wave.x;
float2 displacementY_dz = ih * wave.z;

不过还不够，我们还有两张水平方向的偏移图：

• 
• 

此时不止高度会发生变化，延 x 方向和 z 方向也有一段偏移，这段偏移是

 和 ，对应的偏导

• 
• 

float2 displacementX_dx = -h * wave.x * wave.x * wave.y;
float2 displacementZ_dz = -h * wave.z * wave.z * wave.y;

因此，考虑了水平偏移后，真正的海浪梯度 ​​​​​​​ 就为：

看上去这个公式很复杂，其实本质上就是 (高度图对 x 偏导 / 扭曲挤压后的新海面对 x 偏导, 高度图对 z 偏导 / 扭曲挤压后的新海面对 z 偏导)

---

高度图很好理解，但是平面延 x 和 z 方向的偏移过程可能不是特别好理解：这里可以沿用一张 wiki 上的图片来做参考：

可以看到，在应用水平偏移后，整个海平面就不再是一个线性的平面了，其坐标会被扭曲，甚至会出现挤压（有向面元在变换后面积为负数），

这样再看上面那个极其复杂的海浪梯度向量，其实分子部分就是

分母部分即是对上图右部分的新曲面求偏导，即 

也可以以实际应用为例子，直接看最终海洋的顶点变换表现，当你没有应用水平偏移时，最终顶点动画的效果如左：

可以看出顶点网格水平方向上其实是静止的，只有高度在发生变化，而如果应用了水平方向的偏移，最终顶点动画效果如右

差别明显，如果我们把视角调成垂直往下，就可以直接观察到变换后海面的顶点水平分布情况：

搞定梯度后，想要计算法线就很简单了，只需要拿上向量  减去梯度向量，得到的就是法向量，这样下来对于空间中海面上的一点

，即

其精确法线（未考虑归一化）就为

搞定

5.3.4 IFFT 蝶形图生成

第三章已经对 FFT 蝶形算法做过了详细的解析，一个 N = 8 的蝶形算法如下图：

对于 IDFT ，，令 ，其满足性质：

•  ，
• ，

其中 ，

public void IFFT2D(RenderTexture input, RenderTexture buffer, bool outputToInput = false, bool scale = true, bool permute = false)
{
    int logSize = (int)Mathf.Log(size, 2);
    bool pingPong = false;
    fftShader.SetTexture(KERNEL_HORIZONTAL_STEP_IFFT, PROP_ID_PRECOMPUTED_DATA, precomputedData);
    fftShader.SetTexture(KERNEL_HORIZONTAL_STEP_IFFT, PROP_ID_BUFFER0, input);
    fftShader.SetTexture(KERNEL_HORIZONTAL_STEP_IFFT, PROP_ID_BUFFER1, buffer);
    for (int i = 0; i < logSize; i++)
    {
        pingPong = !pingPong;
        fftShader.SetInt(PROP_ID_STEP, i);
        fftShader.SetBool(PROP_ID_PINGPONG, pingPong);
        fftShader.Dispatch(KERNEL_HORIZONTAL_STEP_IFFT, size / LOCAL_WORK_GROUPS_X, size / LOCAL_WORK_GROUPS_Y, 1);
    }
    fftShader.SetTexture(KERNEL_VERTICAL_STEP_IFFT, PROP_ID_PRECOMPUTED_DATA, precomputedData);
    fftShader.SetTexture(KERNEL_VERTICAL_STEP_IFFT, PROP_ID_BUFFER0, input);
    fftShader.SetTexture(KERNEL_VERTICAL_STEP_IFFT, PROP_ID_BUFFER1, buffer);
    for (int i = 0; i < logSize; i++)
    {
        pingPong = !pingPong;
        fftShader.SetInt(PROP_ID_STEP, i);
        fftShader.SetBool(PROP_ID_PINGPONG, pingPong);
        fftShader.Dispatch(KERNEL_VERTICAL_STEP_IFFT, size / LOCAL_WORK_GROUPS_X, size / LOCAL_WORK_GROUPS_Y, 1);
    }
    
    if (pingPong)
    {
        Graphics.Blit(buffer, input);
    }
    fftShader.SetInt(PROP_ID_SIZE, size);
    fftShader.SetTexture(KERNEL_PERMUTE, PROP_ID_BUFFER0, outputToInput ? input : buffer);
    fftShader.Dispatch(KERNEL_PERMUTE, size / LOCAL_WORK_GROUPS_X, size / LOCAL_WORK_GROUPS_Y, 1);
}

如果递归计算上面按所有的 X, G, H，总体复杂度为 O(NlogN)，例如 N = 8 时，从上图可以看出其递归深度为3，每次计算量都为8，即总共 3x8=24 个蝶形计算单元：

一次蝶形计算单元如上，由于我们可能需要进行多次 IFFT，但事实上由于每次 IFFT 的大小（采样次数）N 是固定的，因此可以预计算所有的 ，并将其复数结果放入一张 logN x N 大小的纹理的 R 和 G 通道中（纹理坐标 x, y 取值从0开始），令 ，则对于纹理  位置的值  满足：

还有两个通道 B 和 A，用于放入图中的上一层 IFFT 的索引 a 和 b，a 和 b 的计算公式如下：

• 
• 

对于 N = 8 的情况下，a, b, k 和 d 的取值即如下，纵轴对应 y = 0~8，横轴对应 x = 0~2

• x = 0：a = 0 1 2 3 0 1 2 3  b = 4 5 6 7 4 5 6 7  k = 0 0 0 0 4 4 4 4  d = 4
• x = 1：a = 0 1 4 5 0 1 4 5  b = 2 3 6 7 2 3 6 7  k = 0 0 2 2 4 4 8 8  d = 2
• x = 2：a = 0 2 4 6 0 2 4 6  b = 1 3 5 7 1 3 5 7  k = 0 1 2 3 4 5 6 7  d = 1

对于实际应用的例子中 N = 256 的情况，生成的蝶形纹理如下：

void PrecomputeTwiddleFactorsAndInputIndices(uint3 id : SV_DispatchThreadID)
{
    uint b = Size >> (id.x + 1);
    float2 mult = 2 * PI * float2(0, 1) / Size;
    uint i = (2 * b * (id.y / b) + id.y % b) % Size;
    float2 twiddle = ComplexExp(-mult * ((id.y / b) * b));
    PrecomputeBuffer[id.xy] = float4(twiddle.x, twiddle.y, i, i + b);
    PrecomputeBuffer[uint2(id.x, id.y + Size / 2)] = float4(-twiddle.x, -twiddle.y, i, i + b);
}

这个无论是 k 的计算，还是最后索引 a 的计算，可能单看数字不是特别明白，但是没有关系，只要代入到后面的计算场景就会立刻清晰：

如果采样递归的手段计算 FFT，那么逻辑写起来会相对简单一点，因为它是顺着思路的，之前就有举过一道多项式乘法的例子，并且给出了代码，然而由于我们想要 GPU 帮我们并行计算，就不好递归去做，这个时候只能逆过来算，也就是按照蝶形图从左至右（stage 1~3 的顺序）的计算每一个单元：这个时候再看我们关于索引 a, b 以及 k 的计算就好理解多了

当然每一层计算完毕后，索引都会按照该层的计算顺序更新排列

不仅如此，由于我们计算每一个蝶形单元时，只依赖上一个层(stage)的结果，因此每一层的运算都可以并行处理，效率极高：

void HorizontalStepInverseFFT(uint3 id : SV_DispatchThreadID)
{
    float4 data = PrecomputedData[uint2(Step, id.x)];
    uint2 inputsIndices = (uint2)data.ba;
    if (PingPong)
    {
        Buffer1[id.xy] = Buffer0[uint2(inputsIndices.x, id.y)]
            + ComplexMult(float2(data.r, -data.g), Buffer0[uint2(inputsIndices.y, id.y)]);
    }
    else
    {
        Buffer0[id.xy] = Buffer1[uint2(inputsIndices.x, id.y)]
            + ComplexMult(float2(data.r, -data.g), Buffer1[uint2(inputsIndices.y, id.y)]);
    }
}

对所有行来一次横向 FFT 后，就是所有列的纵向 FFT，由于对应的计算方式几乎完全一致，就不再做详细介绍了

5.3.5 收尾

到此整个 FFT Ocean 的物理计算流程就要到尾声了，经过 IFFT 计算后，我们就能通过后续纹理合并，以及向量计算拿到：

1. 一张海浪高度图
2. 水平偏移图
3. 法线图

其中①和②可以合到一张纹理中，然后运用到顶点变换：

float3 displacement = 0;
displacement += tex2Dlod(_Displacement_c0, worldUV / LengthScale0) * lod_c0;
v.vertex.xyz += mul(unity_WorldToObject, displacement);

法线图的计算则如下

float4 derivatives = tex2D(_Derivatives_c0, IN.worldUV / LengthScale0);
float2 slope = float2(derivatives.x / (1 + derivatives.z),
    derivatives.y / (1 + derivatives.w));
float3 worldNormal = normalize(float3(-slope.x, 1, -slope.y));

搞定，效果放在了文章最前面，代码的算法部分实际上是直接借鉴的 gasgiant/FFT-Ocean，当然为了能让高配手机上也能运行，以及要做非真实感的海洋，还是要做不少工作的，例如着色，水体交互、海浪泡沫、性能优化等等。考虑到本篇的原意只是介绍 FT 及其在游戏开发领域的应用，就先在这画上句号吧