接上文：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

既然任意连续且收敛的函数 $f(x)$ 都可以分解为无数个不同频率、不同幅值的正、余弦信号，且这个性质可以扩展到二维，即任意一个二维图像都可以分解为无数个复平面波 $e^{-i 2 \pi(u x+v y)}$ 的叠加

$F(u, v)=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x, y) e^{-i 2 \pi(u x+v y)} d x d y$

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

$h(\vec{x}, t)=\sum_k \tilde{h}(\vec{k}, t) e^{i \vec{k} \cdot \vec{x}}$

可以看出，这正是二维离散傅里叶变换的逆变换(IDFT)版本，只不过多了一个量度即时间，考虑到 DFT 形式下任意 $\vec{k}$ 将对应所有可能的复平面波 $e^{-i 2 \pi(n x+m z)}$ ，就能顺理成章的推出波矢量

$\vec{k}=\left(k, l\right)=\left(\frac{2 \pi k}{L}, \frac{2 \pi l}{L}\right)$

其中常量 $L$ 为海平面大小，采样点 $k, l$ 满足 $-N / 2 \leq k, l<N / 2$ ，N 为采样点数量，当然也可以做个变换，即将范围映射到 $0<k, l<N-1$ ，此时

$\vec{k}=\left(\frac{2 \pi k-\pi N}{L}, \frac{2 \pi l-\pi N}{L}\right)$ 本质不变

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

$h(x, z, t)=\frac{1}{N^2} \sum_{k=0}^{N-1} \sum_{l=0}^{N-1} \tilde{h}(k, l, t)e^{i \frac{2 \pi k x-\pi N x}{N}}e^{ i \frac{2 \pi l z-\pi N z}{N}} \\ =\frac{1}{N^2} \sum_{k=0}^{N-1} \sum_{l=0}^{N-1} \tilde{h}(k, l, t)(-1)^x e^{ i \frac{2 \pi kx}{N}}(-1)^z e^{i \frac{2 \pi lz}{N}} \\ \frac{1}{N^2}(-1)^x \sum_{k=0}^{N-1}\left[(-1)^z \sum_{l=0}^{N-1} \tilde{h}(k, l, t) e^{i \frac{2 \pi lz}{N}}\right] e^{i \frac{2 \pi kx}{N}}$

当然也可以任 $N\neq L$ ，这在计算方式上是相同的

5.1.1 二维 IFFT

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

$h(x, z, t)=\frac{1}{N^2}(-1)^x \sum_{k=0}^{N-1}\left[(-1)^z \sum_{l=0}^{N-1} \tilde{h}(k, l, t) e^{i \frac{2 \pi lz}{N}}\right] e^{i \frac{2 \pi kx}{N}}$

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

$h'(k,z,t)=\frac{(-1)^z}{N} \sum_{l=0}^{N-1} \tilde{h}(k, l, t) e^{i \frac{2 \pi lz}{N}}$
$h(x, z, t)=\frac{(-1)^x}{N} \sum_{k=0}^{N-1}h'(k,z,t)e^{i \frac{2 \pi kx}{N}}$

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

计算 $\tilde{h}(k, l, t)$ 生成频谱，如果是交给 GPU 计算，可以将结果存储在一张 NxN 的纹理中（关于 $\tilde{h}(k, l, t)$ 的具体内容及计算方式可以参考下一节 5.2 的内容），t 可以理解成第几帧，必然每一帧这个图像都会不同
执行 N 次水平 FFT： $\sum_{l=0}^{N-1} \tilde{h}(k, l, t) e^{i \frac{2 \pi lz}{N}}$ ，上一步的纹理作为输入，输出结果作为下一步的输入
执行 N 次纵向 FFT： $\sum_{k=0}^{N-1}h'(k,z,t)e^{i \frac{2 \pi kx}{N}}$ ，得到的结果再进行最后的符号计算：即乘上 $\frac{(-1)^x(-1)^y}{N}$ ，搞定

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

5.2 海浪统计模型

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

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

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

5.2.1 菲利浦频谱(Phillips spectrum)

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

对于菲利浦频谱 $P_h(\vec{k})=S(\omega)D(\omega, \theta)=A \frac{e^{-1 /(k L)^2}}{k^4}|\vec{k} \cdot \vec{\omega}|^2$ ，其中

$L=\frac{V^2}{ g}$ ， $V$ 是风速， $\vec{\omega}$ 为风向， $g$ 为引力常数 $9.8 \mathrm{~m} / \mathrm{sec}^2$
$k$ 为 $\vec{k}$ 的幅值，即频率大小

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

非定向波谱 $S(\omega)=A\frac{e^{-1 /(k L)^2}}{k^4}$ ，由海洋统计分析得出
方向扩展函数 $D(\omega, \theta)=A \frac{e^{-1 /(k L)^2}}{k^4}$ ，体现的是风向对波最终幅值的贡献

5.2.2 海洋初始状态

回到前面的 $\tilde{h}(\vec{k}, t)$ ，可以先考虑它的初始状态，也就是 t = 0 的时刻的振幅 $\tilde{h_0}(\vec{k_0})$ 有：

$\tilde{h}_0(\vec{k})=\frac{1}{\sqrt{2}}\left(\xi_r+i \xi_i\right) \sqrt{P_h(\vec{k})}$

$\xi_r$ 和 $\xi_i$ 为两个相互独立的服从均值为0，标准差为1的高斯随机数，实部表初始幅度，虚部表初始相位，而 $P_h(\vec{k})$ 正是前面的方向波谱

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

$h_0=\frac{1}{\sqrt{2}}\left(\xi_1+i \xi_2\right) \sqrt{2 S(\omega) D(\theta, \omega) \frac{d \omega(k)}{d k} \frac{1}{k} \Delta k_x \Delta k_z}$

5.2.3 海洋公式最终形式

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

$\tilde{h}(\vec{k}, t)=\tilde{h}_0(\vec{k}) e^{i \omega(k) t}+\tilde{h}_0^*(-\vec{k}) e^{-i \omega(k) t}$ ，这里又有

角频率与波长的频散关系 $\omega(k)$ 满足 $\omega(k)=\sqrt{g k}$ ，其考虑了引力风波和重力场的关系，当然如果额外考虑水深 h，其频散关系 $\omega(k)=\sqrt{g k \tanh (k h)}$
$\tilde{h}_0$ 是 $\tilde{h}_0^*$ 的共轭复数

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

5.3 基于 GPGPU 的 Ocean IFFT 计算

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

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

5.3.1 初始频谱计算

先是计算菲利浦频谱 $P_h(\vec{k})=A \frac{e^{-1 /(k L)^2}}{k^4}|\vec{k} \cdot \vec{\omega}|^2$ 以生成对应纹理，看上去最陌生的公式其实直接抄作业就行了，其中 $\vec{k}=\left(k, l\right)$ 作为纹理坐标，其余除了风向 $\vec{\omega}$ 都为常量。很明显，游戏中我们不太可能每帧去改变风向，因此风向 $\vec{\omega}$ 也可以视作常量（由美术或策划去配置），在这种情况下 $P_h(\vec{k})$ 的计算理论可以离线或者只在风向改变时做一次，不需要每帧去重复计算

然后就是高斯随机数 $\left(\xi_r+i \xi_i\right)$ 生成，这个只是为了得到不同的初始状态，因此同样可以离线计算：

//如果本地已经存在 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的高斯随机数，公式如下：

$\xi_r =\sin \left(2 \pi u\right) \sqrt{-2 \log \left(v\right)}$
$\xi_i =\cos \left(2 \pi u\right) \sqrt{-2 \log \left(v\right)}$

其中 $u,v$ 为均匀分布的随机数

然后两者相组合，就能得到 $\tilde{h}_0(\vec{k})$ ：

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 任意时刻频谱计算

关于 $\tilde{h}(\vec{k}, t)=\tilde{h}_0(\vec{k}) e^{i \omega(k) t}+\tilde{h}_0^*(-\vec{k}) e^{-i \omega(k) t}$ 的代码实现更简单，毕竟前面已经得到了 $\tilde{h}_0(\vec{k})$ 以及 $\omega(k)$ （该图片来源于 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 轴的偏移图，这里可以先根据 $\tilde{h}(\vec{k}, t)$ 获取水平方向的频谱

$\tilde{d_x}(\vec{k}, t)=i \frac{k_x}{|k|} \tilde{h}(\vec{k}, t)$
$\tilde{d_z}(\vec{k}, t)=i \frac{k_z}{|k|} \tilde{h}(\vec{k}, t)$

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

前面代码将 $(k_x, \frac{1}{|k|}, k_z, w(k))$ 一起塞到了一张4通道纹理中，这里继续拿来用：

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

5.3.3 海洋法线

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

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

得到高度图后，可以通过差分方法来求法线，即对变换后的新顶点 $M_0(x,y)$ ，求出与其相邻的两组顶点 $M_1(x+1,y)$ $M_1'(x-1,y)$ 和 $M_2(x,y-1)$ $M_2'(x,y+1)$ ，交叉连线后得到两个向量，叉乘得出法线，
通过偏导计算求出曲面梯度，然后再根据梯度得到法线

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

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

设函数 $z=f(x, y)$ 在平面区域 D 内具有一阶连续偏导数，则对于每一点 $(x, y)\in D$ ，都可定出一个向量 $\frac{\partial f}{\partial x} i+\frac{\partial f}{\partial y} j$ ，该向量称为函数 $z=f(x, y)$ 在点 $P(x, y)$ 的梯度，记作 $\nabla f(x, y)$

可以证明： $\nabla f(x, y)=\left\{\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}\right\}=f_x(x, y) \bar{i}+f_y(x, y) \bar{j}$

对于高度图 $h(\vec{x}, t)=\sum_k \tilde{h}(\vec{k}, t) e^{i \vec{k} \cdot \vec{x}}$ ，其梯度满足

$\nabla h(\vec{x}, t)=\sum_{{k}} \tilde{h}(\vec{k}, t) \nabla e^{i \vec{k} \cdot \vec{x}} = \sum_{{k}} i\vec{k}\tilde{h}(\vec{k}, t) e^{i \vec{k} \cdot \vec{x}}$ ，其中偏导

$\frac{\partial h(\vec{x}, t)}{\partial x}=\sum_{{k}} ik_x\tilde{h}(\vec{k}, t) e^{i \vec{k} \cdot \vec{x}}$
$\frac{\partial h(\vec{x}, t)}{\partial z}=\sum_{{k}} ik_z\tilde{h}(\vec{k}, t) e^{i \vec{k} \cdot \vec{x}}$

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

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

$D_x(\vec{x}, t)=\sum_ki \frac{k_x}{k} \tilde{h}(\vec{k}, t) e^{i \vec{k} \cdot \vec{x}}$
$D_z(\vec{x}, t)=\sum_ki \frac{k_z}{k} \tilde{h}(\vec{k}, t) e^{i \vec{k} \cdot \vec{x}}$

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

$D_x'(\vec{x}, t)=x+ \vec{D_x}(\vec{x}, t)$ 和 $D_z'(\vec{x}, t)=z+ \vec{D_z}(\vec{z}, t)$ ，对应的偏导

$\frac{\partial D_x'(\vec{x}, t)}{\partial x}=1+ \sum_k- \frac{k_x^2}{k} \tilde{h}(\vec{k}, t) e^{i \vec{k} \cdot \vec{x}}$
$\frac{\partial D_z'(\vec{x}, t)}{\partial z}=1+ \sum_k- \frac{k_z^2}{k} \tilde{h}(\vec{k}, t) e^{i \vec{k} \cdot \vec{x}}$

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

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

$(\frac{\sum_{{k}} ik_x\tilde{h}(\vec{k}, t) e^{i \vec{k} \cdot \vec{x}}}{1+ \sum_k- \frac{k_x^2}{k} \tilde{h}(\vec{k}, t) e^{i \vec{k} \cdot \vec{x}}} ,\frac{\sum_{{k}} ik_z\tilde{h}(\vec{k}, t) e^{i \vec{k} \cdot \vec{x}} }{1+ \sum_k- \frac{k_z^2}{k} \tilde{h}(\vec{k}, t) e^{i \vec{k} \cdot \vec{x}}})$

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

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

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

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

$\nabla h(\vec{x}, t)=(\frac{\partial h(\vec{x}, t)}{\partial x},\frac{\partial h(\vec{x}, t)}{\partial z})$

分母部分即是对上图右部分的新曲面求偏导，即 $(\frac{\partial D_x'(\vec{x}, t)}{\partial x},\frac{\partial D_z'(\vec{x}, t)}{\partial z})$

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

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

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

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

$P(x+ \vec{D_x}(\vec{x}, t), h(\vec{x}, t), z+ \vec{D_z}(\vec{z}, t))$ ，即

$P(x+\sum_ki \frac{k_x}{k} \tilde{h}(\vec{k}, t) e^{i \vec{k} \cdot \vec{x}}, \sum_k \tilde{h}(\vec{k}, t) e^{i \vec{k} \cdot \vec{x}}, z+\sum_ki \frac{k_z}{k} \tilde{h}(\vec{k}, t) e^{i \vec{k} \cdot \vec{x}})$

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

$(0,1,0)-(\frac{\frac{\partial h(\vec{x}, t)}{\partial x}}{1 + \frac{\partial D_x(\vec{x}, t)}{\partial x}}, 0, \frac{\frac{\partial h(\vec{x}, t)}{\partial z}}{1 + \frac{\partial D_z(\vec{x}, t)}{\partial z}}) \\ = P_{N}(-\frac{\sum_{{k}} ik_x\tilde{h}(\vec{k}, t) e^{i \vec{k} \cdot \vec{x}}}{1+ \sum_k- \frac{k_x^2}{k} \tilde{h}(\vec{k}, t) e^{i \vec{k} \cdot \vec{x}}}, 1, -\frac{\sum_{{k}} ik_z\tilde{h}(\vec{k}, t) e^{i \vec{k} \cdot \vec{x}} }{1+ \sum_k- \frac{k_z^2}{k} \tilde{h}(\vec{k}, t) e^{i \vec{k} \cdot \vec{x}}})$

搞定

5.3.4 IFFT 蝶形图生成

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

对于 IDFT $x[n]=\frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{i \frac{2 \pi}{N} kn}$ ， $n(0{\leq}n{<}N)$ ，令 $W_N=e^{-i\frac{2{\pi}}{N}}$ ，其满足性质：

$X(n)=G(n)+W_N^{-n} H(n)$ ， $0{\leq}n{<}\frac{N}{2}$
$X(n)=G(n-\frac{N}{2})-W_N^{-(n-\frac{N}{2})} H(n-\frac{N}{2})$ ， $\frac{N}{2}{\leq}n{<}N$

其中 $G(n)=\sum_{r=0}^{\frac{N}{ 2}-1} x(2 r) W_{\frac{N}{ 2}}^{-r n}$ ， $H(n)=\sum_{r=0}^{\frac{N}{ 2}-1} x(2 r+1) W_{\frac{N}{ 2}}^{-r n}$

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 是固定的，因此可以预计算所有的 $W_N^k$ ，并将其复数结果放入一张 logN x N 大小的纹理的 R 和 G 通道中（纹理坐标 x, y 取值从0开始），令 $d = \frac{N}{2^{x+1}}$ ，则对于纹理 $(x,y)$ 位置的值 $W_N^k$ 满足：

$k=\lfloor\frac{y }{d} \rfloor\cdot d$

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

$a=(2k+y\bmod d)\bmod{N}$
$b=a+d$

对于 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 计算后，我们就能通过后续纹理合并，以及向量计算拿到：

一张海浪高度图
水平偏移图
法线图

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

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 及其在游戏开发领域的应用，就先在这画上句号吧