DirectX11实现FFT海面模拟

FFT_Ocean Simulation with DirectX11

概述：

模拟海洋的方法主要就是两种：

利用Gerstner波（即通过在水平方向上挤压正弦波叠加形成的海洋波谱，来形成形成较尖的浪头和较宽的浪槽），当然本文的重点不在这里，具体实现可查阅以下博客
1. 【DX11学习笔记】GerstnerWave波浪模拟（基于GPU计算着色器的实现） - Neko_YG - 博客园 (cnblogs.com)
2. DirectX11 GerstnerWaves模拟（计算着色器） - YIMG - 博客园 (cnblogs.com)
利用海洋频谱，再通过逆离散傅里叶变化（IDEF）以及我们的核心算法快速傅里叶变化（FFT）来获得由许多不同频率、振幅的波叠加出来的复杂海面。本文将着重介绍如何在使用Direct3D 11 API来实现FFT海面模拟

前置知识：

1. 时域、频域、（逆）傅里叶变化

我们知道，一个复杂的波形可以通过非常多的正弦波叠加来逼近，首先摆出一张很直观的图片

可以看到这里有非常多的坐标轴，所以当我们从不同的角度来观察时，就可以得到不同的含义。

当我们从左下角往右上角观察（把波形图像往频率方向挤压），就可以得到右上角的时域图像。

当我们从右下角往左上角观察（把波形图像往时间轴的负方向挤压），就可以得到左上角的频域图像。

而时域和频域是可以转化的，给我们一个时域的图像（函数），可以通过傅里叶变换转化为频域图像

给我们一个频域的图像（函数），可以通过逆傅里叶变换转化为时域图像

所以，如果我们想要实现一个复杂的波浪，我们只需要生成一个频域图像(也被称为频谱)，然后对其进行逆傅里叶变换就可以得到它的时域图像（也就是我们想要的海洋波形图）

2.离散傅里叶变换

因为本人对数学一窍不通，所以先直接摆出公式来用，具体的推理可以见以下两篇文章（第一篇为上图出处）

【学习笔记】Unity 基于GPU FFT海洋的实现-理论篇 - 知乎 (zhihu.com)
fft海面模拟(一) - 知乎 (zhihu.com)

F(μ)=∑x=0N−1f(x)e−i2πμxNF(\mu)=\sum_{x=0}^{N-1}f(x)e^{-i\frac{2\pi\mu x}{N}} F(μ)=x=0∑N−1f(x)e−iN2πμx

以上就是离散傅里叶变换的公式

F(μ)F(\mu)F(μ)为转换后的频域函数，μ\muμ是频率。可以看上面的图，我们将一个复杂的波形展开成多个正弦波，在将正弦波按照频率大小依次排列，频率和振幅就组成了频域图像。频域图像的横坐标为频率，纵坐标为振幅。f(x)f(x)f(x)是我们时域的函数，e−i2πμxNe^{-i\frac{2\pi\mu x}{N}}e−iN2πμx是一个复数。

其余就不做过多介绍了（也不会…）

而我们实际上使用的是逆离散傅里叶变化，公式其实差不多
f(x)=1N∑μ=0N−1F(μ)ei2πμxNf(x)=\frac{1}{N}\sum_{\mu=0}^{N-1}F(\mu)e^{i\frac{2\pi\mu x}{N}} f(x)=N1μ=0∑N−1F(μ)eiN2πμx
b站上一个介绍FFT的视频里也有一部分介绍DFT和IDFT,以下是视频中的一张截图

快速傅里叶变换(FFT)——有史以来最巧妙的算法？_哔哩哔哩_bilibili

核心内容与具体实现

因为在网上可以搜索到很多理论方面的，写的很详细的教程，这里就不介绍过多的理论知识，直接从代码（实践）的角度介绍,中间会夹杂少量的理论知识和直接摆出来使用的公式

代码部分是在x_jun大佬的dx教程源码的基础上完成的。

教程链接：DirectX11 With Windows SDK–00 目录 - X_Jun - 博客园 (cnblogs.com)

整体思维导图：

（具体内容随后一一介绍）

初始化：

海洋网格和纹理参数初始化

因为我们看到的海洋表面本质上就是一个巨大的海洋的网格平面。其中的每个方形网格由两个三角形构成。我们辛辛苦苦算出偏移纹理也是为了给网格上的顶点用来采样，最终确定顶点的位置（通过网上上顶点在空间中的排布，再贴上表面的海水纹理和泡沫纹理等，最终绘制出海面）。所以在初始化时，我们首先要创建一个海洋的网格

//=========创建海洋网格==========
//
// 顶点行(列)数 - 1 = 网格行(列)数
auto meshData = Geometry::CreateTerrain<VertexPosNormalTex, DWORD>(XMFLOAT2((cols - 1) * spatialStep, (rows - 1) * spatialStep),XMUINT2(cols - 1, rows - 1));//利用教程项目中内部实现的CreateTerrain方法来创建网格m_IndexCount = 6 * (rows - 1) * (cols - 1);//计算顶点（索引）数量

这里的spatialStep（空间步长），就是指网格中相邻顶点之间的距离，可以配合这顶点行数rows和列数cols来设置海洋网格的大小

随后我们再设置好纹理坐标再u,v方向上的单位长度以及纹理的阶次（FFT计算要用到）

海洋频谱参数初始化

FFT海洋公式

如果我们想要得到一个复杂的海面，我们需要生成一个频谱(也就是上面的频域)，然后在通过逆离散傅里叶变换(IDFT)，我们就可以得到由很多不同频率、振幅的波所叠加出来的复杂海面。频谱直接决定了我们的海面形状，这里采用的是一篇论文（Simulating Ocean Water-Jerry Tessendorf）中的公式
h(x→,t)=∑kh~(k~,t)eik~x→h(\overrightarrow{x},t)=\sum_{k}\widetilde{h}(\widetilde{k},t)e^{i\widetilde{k}\overrightarrow{x}} h(x,t)=k∑h(k,t)eikx
海洋公式这块也不知道怎么解释原理，个人理解就是海洋学家们找到了利用风速，风向和一些随机数来模拟海浪的方法（函数），我们只需要知道怎么代入参数计算即可。

（想要深入了解可以参考上面第一篇Unity的文章或是查下论文“Empirical Directional Wave Spectra for Computer Graphics”）

代码实现：

因为这些海洋相关的参数都是要传入到着色器中作为全局变量计算的，所以这里就根据更新的频率将参数打包成了三个数据结构体，并在HLSL中创建对应的结构体

//只需要传递一次的cb
struct {float m_SpatialStep = 0.0f;            // 空间步长int fftSize;   //fft纹理大小DirectX::XMFLOAT2 g_pad1;
} m_CBInitSettings = {};//更新时维护的参数结构体
struct {DirectX::XMFLOAT4 WindAndSeed = { 0.1f, 0.2f, 0.0, 0.0 };//风向和随机种子 xy为风, zw为两个随机种子float Lambda = -1.0;  //控制偏移大小float HeightScale = 1.0;  //高度影响float BubblesScale = 1.0;  //泡沫强度float BubblesThreshold = 1.0;  //泡沫阈值float A = 10;            //海浪高度因子float Time = 0.0f;        // 累积时间DirectX::XMFLOAT2 g_pad2;
} m_CBUpdataSettings = {};// 对应Waves.hlsli的常量缓冲区//专门传递Ns的结构体（用于FFT计算时的阶数）
struct {int ns;DirectX::XMFLOAT3 g_pad3;
}m_CBns = {};

通过海洋公式，我们就可以得到我们的初始高度频谱

（上图为利用Philips频谱和高斯随机数生成初始频谱的过程）

再加入时间的维度后，我们进而得到海洋高度公式以及水平偏移公式。有了这三个公式后，在某一时刻（固定时间变量t），通过代入波矢量k→\overrightarrow{k}k 就可以得到在某一时刻x,y,z三个方向的偏移频谱（也就是有许多不同的波叠加形成的时域图像）

当我们拿到这些频谱后，分别对他们进行IDFT，就会得到水平x,zx,zx,z及高度yyy的偏移图（中间三个灰白的纹理）

将三张xyz三个方向的偏移纹理合在一起变成一张最终用于顶点着色器采样的纹理，也就是我们的位移贴图（彩色的Displacement那张）。有了这样偏移纹理就相当于知道了所有波浪的顶点位置，所以也可以在cs完成法线和泡沫的计算，并也为它们生成两张纹理图。这样在VS和PS阶段就可以直接通过采样来获取法线和泡沫的数据。

代码实现：

填充好描述，创建各种纹理资源，及其相应的无序访问视图（UAV）因为我们生成的纹理是要读写的。因为最终的（上图右边的三个纹理）是要给到VS和PS采样，所以也要为它们创建着色器资源视图（SRV）

关于UAV的介绍：DirectX11 With Windows SDK–26 计算着色器：入门 - X_Jun - 博客园 (cnblogs.com)

//用Comptr管理纹理资源
ComPtr<ID3D11Texture2D> m_pGaussianRandomRT;         //高斯随机数对
ComPtr<ID3D11Texture2D> m_pHeightSpectrumRT;         //高度频谱
ComPtr<ID3D11Texture2D> m_pDisplaceXSpectrumRT;      //x方向偏移频谱
ComPtr<ID3D11Texture2D> m_pDisplaceZSpectrumRT;      //z方向偏移频谱
ComPtr<ID3D11Texture2D> m_pDisplaceRT;               //（总）偏移频谱
ComPtr<ID3D11Texture2D> m_pInputRT;                  //FFT输入纹理
ComPtr<ID3D11Texture2D> m_pOutputRT;                 //FFT输出纹理
ComPtr<ID3D11Texture2D> m_pNormalRT;                 //法线纹理
ComPtr<ID3D11Texture2D> m_pBubblesRT;                //泡沫纹理// 创建缓存计算结果的资源
D3D11_TEXTURE2D_DESC texDesc;
texDesc.Width = cols;
texDesc.Height = rows;
texDesc.MipLevels = 1;
texDesc.ArraySize = 1;
texDesc.Format = DXGI_FORMAT_R16G16B16A16_FLOAT;
texDesc.SampleDesc.Count = 1;
texDesc.SampleDesc.Quality = 0;
texDesc.Usage = D3D11_USAGE_DEFAULT;
texDesc.BindFlags = D3D11_BIND_SHADER_RESOURCE |D3D11_BIND_UNORDERED_ACCESS;
texDesc.CPUAccessFlags = 0;
texDesc.MiscFlags = 0;
hr = device->CreateTexture2D(&texDesc, nullptr, m_pGaussianRandomRT.GetAddressOf());
if (FAILED(hr))return hr;
//...
//为上面每一个纹理都创建资源,此处省略（太长了）
//...// 创建着色器资源视图
D3D11_SHADER_RESOURCE_VIEW_DESC srvDesc;
srvDesc.ViewDimension = D3D11_SRV_DIMENSION_TEXTURE2D;
srvDesc.Format = DXGI_FORMAT_R16G16B16A16_FLOAT;
srvDesc.Texture2D.MipLevels = 1;
srvDesc.Texture2D.MostDetailedMip = 0;
//只需要为三个采样需要用的的纹理资源创建SRV
hr = device->CreateShaderResourceView(m_pDisplaceRT.Get(), &srvDesc, m_pDisplaceSRV.GetAddressOf());
if (FAILED(hr))return hr;
hr = device->CreateShaderResourceView(m_pNormalRT.Get(), &srvDesc, m_pNormalSRV.GetAddressOf());
if (FAILED(hr))return hr;
hr = device->CreateShaderResourceView(m_pBubblesRT.Get(), &srvDesc, m_pBubblesSRV.GetAddressOf());
if (FAILED(hr))return hr;// 创建无序访问视图
D3D11_UNORDERED_ACCESS_VIEW_DESC uavDesc;
uavDesc.ViewDimension = D3D11_UAV_DIMENSION_TEXTURE2D;
uavDesc.Format = DXGI_FORMAT_R16G16B16A16_FLOAT;
uavDesc.Texture2D.MipSlice = 0;
hr = device->CreateUnorderedAccessView(m_pGaussianRandomRT.Get(), &uavDesc, m_pGaussianRandomUAV.GetAddressOf());
if (FAILED(hr))return hr;
//...
//为上面每一个纹理都创建UAV,此处省略（太长了）
//...

生成高斯随机数对

创建好了各种纹理资源后，就可以利用GPU强大的并行计算的能力，在计算着色器CS中生成我们的纹理

首先是生成高斯随机数对。

代码实现：

首先是Waves核心和生成高斯随机数对ComputeGaussianRandom_CS的设计

Waves.hlsli

#define PI 3.14159274f
#define G 9.81f// 用于初始化的cb
cbuffer cbInitSettings : register(b0)
{float OceanLength; //海洋长度int N; //fft纹理大小//uint rngState;float2 g_pad1;
}// 用于更新模拟的cb
cbuffer cbUpdataSettings : register(b1)
{float4 WindAndSeed; //风和随机种子 xy为风, zw为两个随机种子float Lambda; //偏移影响float HeightScale; //高度影响float BubblesScale;  //泡沫强度float BubblesThreshold; //泡沫阈值float A; //phillips谱参数，影响波浪高度float Time; //时间float2 g_pad2;
}//用于FFT计算的阶数
cbuffer cbns : register(b2)
{int Ns;float3 g_pad3;
}//RW纹理
RWTexture2D<float4> GaussianRandomRT : register(u0); //高斯随机数
RWTexture2D<float4> HeightSpectrumRT : register(u1); //高度频谱
RWTexture2D<float4> DisplaceXSpectrumRT : register(u2); //X偏移频谱
RWTexture2D<float4> DisplaceZSpectrumRT : register(u3); //Z偏移频谱
RWTexture2D<float4> DisplaceRT : register(u4); //最后生成的偏移纹理
RWTexture2D<float4> InputRT : register(u5); //输入
RWTexture2D<float4> OutputRT : register(u6); //输出
RWTexture2D<float4> NormalRT : register(u7); //法线纹理
RWTexture2D<float4> BubblesRT : register(u8); //泡沫纹理//=================实用函数声明===================
float DonelanBannerDirectionalSpreading(float2 k);
float PositiveCosineSquaredDirectionalSpreading(float2 k);
float phillips(float2 k);
float dispersion(float2 k);
float2 gaussian(float2 id);
uint wangHash(uint seed);
float rand();
float2 complexMultiply(float2 c1, float2 c2);//=================实用函数定义===================
//计算高斯随机数
//Donelan-Banner方向拓展
float DonelanBannerDirectionalSpreading(float2 k)
{float betaS;float omegap = 0.855f * G / length(WindAndSeed.xy);float ratio = dispersion(k) / omegap;if (ratio < 0.95f){betaS = 2.61f * pow(ratio, 1.3f);}if (ratio >= 0.95f && ratio < 1.6f){betaS = 2.28f * pow(ratio, -1.3f);}if (ratio > 1.6f){float epsilon = -0.4f + 0.8393f * exp(-0.567f * log(ratio * ratio));betaS = pow(10, epsilon);}float theta = atan2(k.y, k.x) - atan2(WindAndSeed.y, WindAndSeed.x);return betaS / max(1e-7f, 2.0f * tanh(betaS * PI) * pow(cosh(betaS * theta), 2));
}
//正余弦平方方向拓展
float PositiveCosineSquaredDirectionalSpreading(float2 k)
{float theta = atan2(k.y, k.x) - atan2(WindAndSeed.y, WindAndSeed.x);if (theta > -PI / 2.0f && theta < PI / 2.0f){return 2.0f / PI * pow(cos(theta), 2);}else{return 0;}
}//计算phillips谱
float phillips(float2 k)
{float kLength = length(k);kLength = max(0.001f, kLength);// kLength = 1;float kLength2 = kLength * kLength;float kLength4 = kLength2 * kLength2;float windLength = length(WindAndSeed.xy);float l = windLength * windLength / G;float l2 = l * l;float damping = 0.001f;float L2 = l2 * damping * damping;//phillips谱return A * exp(-1.0f / (kLength2 * l2)) / kLength4 * exp(-kLength2 * L2);
}
//随机种子
uint wangHash(uint seed)
{seed = (seed ^ 61) ^ (seed >> 16);seed *= 9;seed = seed ^ (seed >> 4);seed *= 0x27d4eb2d;seed = seed ^ (seed >> 15);return seed;
}
//计算均匀分布随机数[0,1)
float rand(uint rng)
{// Xorshift算法rng ^= (rng << 13);rng ^= (rng >> 17);rng ^= (rng << 5);return rng / 4294967296.0f;;
}//计算高斯随机数
float2 gaussian(float2 id)
{//均匀分布随机数uint rngState = wangHash(id.y * N + id.x);float x1 = rand(rngState);float x2 = rand(rngState);x1 = max(1e-6f, x1);x2 = max(1e-6f, x2);//计算两个相互独立的高斯随机数float g1 = sqrt(-2.0f * log(x1)) * cos(2.0f * PI * x2);float g2 = sqrt(-2.0f * log(x1)) * sin(2.0f * PI * x2);return float2(g1, g2);
}
//计算弥散
float dispersion(float2 k)
{return sqrt(G * length(k));
}
//复数相乘
float2 complexMultiply(float2 c1, float2 c2)
{return float2(c1.x * c2.x - c1.y * c2.y,c1.x * c2.y + c1.y * c2.x);
}

ComputeGaussianRandom_CS.hlsl

#include"Waves.hlsli"
//计算高斯随机变量
[numthreads(16, 16, 1)]
void ComputeGaussianRandom(uint3 id : SV_DispatchThreadID)
{float2 g = gaussian(id.xy);GaussianRandomRT[id.xy] = float4(g, 0, 0);
}

然后在c++中创建它（后面的CS的创建方法相同，只在这里贴出来一次）

ComPtr<ID3DBlob> blob;
hr = CreateShaderFromFile(L"HLSL\\ComputeGaussianRandom_CS.cso", L"HLSL\\ComputeGaussianRandom_CS.hlsl", "ComputeGaussianRandom", "cs_5_0", blob.GetAddressOf());
if (FAILED(hr))return hr;
hr = device->CreateComputeShader(blob->GetBufferPointer(), blob->GetBufferSize(), nullptr, m_pCompGaussianRandomCS.GetAddressOf());
if (FAILED(hr))return hr;

因为高斯随机数纹理只需要生成一次，所以我就将它和固定参数的初始化打包成了WaveRender类的一个单独的成员函数。

void WavesRender::InitCS(ID3D11DeviceContext* deviceContext)
{//初始化常量缓冲区//D3D11_MAPPED_SUBRESOURCE mappedData;deviceContext->Map(m_pInitConstantBuffer.Get(), 0, D3D11_MAP_WRITE_DISCARD, 0, &mappedData);memcpy_s(mappedData.pData, sizeof m_CBInitSettings, &m_CBInitSettings, sizeof m_CBInitSettings);deviceContext->Unmap(m_pInitConstantBuffer.Get(),0);//生成高斯随机数////设置CS资源deviceContext->CSSetConstantBuffers(0, 1, m_pInitConstantBuffer.GetAddressOf());deviceContext->CSSetShader(m_pCompGaussianRandomCS.Get(), nullptr, 0);deviceContext->CSSetUnorderedAccessViews(0, 1, m_pGaussianRandomUAV.GetAddressOf(), nullptr);//调用线程运算deviceContext->Dispatch(m_CBInitSettings.fftSize / 16, m_CBInitSettings.fftSize / 16, 1);//解除绑定ID3D11UnorderedAccessView* nullpUAV[1] = { nullptr };deviceContext->CSSetUnorderedAccessViews(0, 1, nullpUAV, nullptr);
}

在绑定纹理视图时需要注意的一点，在使用同时含有UAV和SRV的纹理时，着色器一次只能绑定该纹理的其中一个视图，当其中一种视图使用完后，要记得解除绑定，否则就会出现冲突。

至此，初始化的部分基本就已经完成了（还包括一些着色器的初始化，顶点、索引…缓冲区的初始化，还未出场的CS的初始化，海洋表面纹理的读取等就不介绍了）

更新：

接下来就到了重头戏，也就是波浪的更新了。这里要完成每次更新波浪所用的各种纹理的计算、生成

生成高度频谱和偏移频谱

首先要更新风向（也可以初始化给定值后不更新）。随后就可以结合我们在初始化时得到的高斯随机数谱生成高度频谱，进而得到偏移频谱了

代码实现：

这里先判断累积时间大于时间步长m_TimeStep（自定义的更新时间间隔）才更新，随后更新风向和种子，再利用Map更新HLSL中的数据结构体

m_AccumulateTime += dt;if (m_AccumulateTime < m_TimeStep)return;m_CBUpdataSettings.Time += m_AccumulateTime;//更新时间
//更新海洋相关参数
//
//生成随机数,更新风向
std::uniform_real_distribution<float> u(1.0, 10.0);
m_CBUpdataSettings.WindAndSeed.z = u(e);
m_CBUpdataSettings.WindAndSeed.w = u(e);
XMFLOAT2 wind = { m_CBUpdataSettings.WindAndSeed.x, m_CBUpdataSettings.WindAndSeed.y };
XMVECTOR tmp = XMLoadFloat2(&wind);
tmp = XMVector2Normalize(tmp);
XMStoreFloat2(&wind, tmp);
m_CBUpdataSettings.WindAndSeed.x = wind.x * WindScale;
m_CBUpdataSettings.WindAndSeed.y = wind.y * WindScale;//更新常量缓冲区
//
D3D11_MAPPED_SUBRESOURCE mappedData;
deviceContext->Map(m_pUpdataConstantBuffer.Get(), 0, D3D11_MAP_WRITE_DISCARD, 0, &mappedData);
memcpy_s(mappedData.pData, sizeof m_CBUpdataSettings, &m_CBUpdataSettings, sizeof m_CBUpdataSettings);
deviceContext->Unmap(m_pUpdataConstantBuffer.Get(), 0);

接下来生成高度频谱，并进而生成偏移频谱

//生成高度频谱
//
//设置CS资源
deviceContext->CSSetConstantBuffers(1, 1, m_pUpdataConstantBuffer.GetAddressOf());
deviceContext->CSSetShader(m_pCreateHeightSpectrumCS.Get(), nullptr, 0);
ID3D11UnorderedAccessView* pUAVs2[2] = { m_pGaussianRandomUAV.Get(),m_pHeightSpectrumUAV.Get()};
deviceContext->CSSetUnorderedAccessViews(0, 2, pUAVs2, nullptr);
//调用线程运算
deviceContext->Dispatch(m_CBInitSettings.fftSize / 16, m_CBInitSettings.fftSize / 16, 1);
//解除绑定
pUAVs2[0] = pUAVs2[1] = nullptr;
deviceContext->CSSetUnorderedAccessViews(0, 2, pUAVs2, nullptr);//生成偏移频谱
//
//设置CS资源
deviceContext->CSSetShader(m_pCreateDisplaceSpectrumCS.Get(), nullptr, 0);
ID3D11UnorderedAccessView* pUAVs3[3] = { m_pHeightSpectrumUAV.Get(),m_pDisplaceXSpectrumUAV.Get(),m_pDisplaceZSpectrumUAV.Get() };
deviceContext->CSSetUnorderedAccessViews(1, 3, pUAVs3, nullptr);
//调用线程运算
deviceContext->Dispatch(m_CBInitSettings.fftSize / 16, m_CBInitSettings.fftSize / 16, 1);
//解除绑定
pUAVs3[0] = pUAVs3[1] = pUAVs3[2] = nullptr;
deviceContext->CSSetUnorderedAccessViews(1, 3, pUAVs3, nullptr);

下面是生成高度频谱CreateHeightSpectrum_CS和CreateDisplaceSpectrum_CS的设计

#include"Waves.hlsli"//生成高度频谱
[numthreads(16, 16, 1)]
void CreateHeightSpectrum(uint3 id : SV_DispatchThreadID)
{float2 k = float2(2.0f * PI * id.x / N - PI, 2.0f * PI * id.y / N - PI);float2 gaussian = GaussianRandomRT[id.xy].xy;float2 hTilde0 = gaussian * sqrt(abs(phillips(k) * DonelanBannerDirectionalSpreading(k)) / 2.0f);float2 hTilde0Conj = gaussian * sqrt(abs(phillips(-k) * DonelanBannerDirectionalSpreading(-k)) / 2.0f);hTilde0Conj.y *= -1.0f;float omegat = dispersion(k) * Time;float c = cos(omegat);float s = sin(omegat);float2 h1 = complexMultiply(hTilde0, float2(c, s));float2 h2 = complexMultiply(hTilde0Conj, float2(c, -s));float2 HTilde = h1 + h2;HeightSpectrumRT[id.xy] = float4(HTilde, 0, 0);
}

#include"Waves.hlsli"
//生成偏移频谱
[numthreads(16, 16, 1)]
void CreateDisplaceSpectrum(uint3 id : SV_DispatchThreadID)
{float2 k = float2(2 * PI * id.x / N - PI, 2 * PI * id.y / N - PI);k /= max(0.001f, length(k));float2 HTilde = HeightSpectrumRT[id.xy].xy;float2 KxHTilde = complexMultiply(float2(0, -k.x), HTilde);float2 kzHTilde = complexMultiply(float2(0, -k.y), HTilde);DisplaceXSpectrumRT[id.xy] = float4(KxHTilde, 0, 0);DisplaceZSpectrumRT[id.xy] = float4(kzHTilde, 0, 0);
}

FFT计算xyz方向的偏移纹理

现在回到我们的海洋高度公式，h(x→,t)=∑kh~(k~,t)eik~x→h(\overrightarrow{x},t)=\sum_{k}\widetilde{h}(\widetilde{k},t)e^{i\widetilde{k}\overrightarrow{x}}h(x,t)=∑kh(k,t)eikx，将其展开成标量形式
h(x,z,t)=∑m=−M2M2−1∑n=−N2N2−1h~(2πnLx,2πmLz,t)ei(2πnLxx+2πmLzz)h(x,z,t)=\sum_{m=-\frac{M}{2}}^{\frac{M}{2}-1}\sum_{n=-\frac{N}{2}}^{\frac{N}{2}-1}\widetilde{h}(\frac{2\pi n}{L_x},\frac{2\pi m}{L_z},t)e^{i(\frac{2\pi n}{L_x}x+\frac{2\pi m}{L_z}z)} h(x,z,t)=m=−2M∑2M−1n=−2N∑2N−1h(Lx2πn,Lz2πm,t)ei(Lx2πnx+Lz2πmz)
不难看出，这里我们要进行N2M2N^2M^2N2M2次乘法运算，时间复杂度是平方级的，而快速傅里叶变化（FFT）算法就是用于快速计算离散傅里叶变换的算法（注：FFT只用于计算N为2的幂的DFT)，算法大概的思路就是。将要累加的N个乘法结果按序号分成奇偶两组。像下图那样

（图片来源：fft海面模拟（二） - 知乎 (zhihu.com)）

则有
G(k)=∑n=0N2−1g(n)e−i2πkN2n=∑n=0N2−1x(2n)e−i2πkN2n,k∈{0,1,...,N2−1}G(k)=\sum_{n=0}^{\frac{N}{2}-1}g(n)e^{-i\frac{2\pi k}{\frac{N}{2}}n}=\sum_{n=0}^{\frac{N}{2}-1}x(2n)e^{-i\frac{2\pi k}{\frac{N}{2}}n},k\in\{0,1,...,\frac{N}{2}-1\} G(k)=n=0∑2N−1g(n)e−i2N2πkn=n=0∑2N−1x(2n)e−i2N2πkn,k∈{0,1,...,2N−1}

H(k)=∑n=0N2−1h(n)e−i2πkN2n=∑n=0N2−1x(2n+1)e−i2πkN2n,k∈{0,1,...,N2−1}H(k)=\sum_{n=0}^{\frac{N}{2}-1}h(n)e^{-i\frac{2\pi k}{\frac{N}{2}}n}=\sum_{n=0}^{\frac{N}{2}-1}x(2n+1)e^{-i\frac{2\pi k}{\frac{N}{2}}n},k\in\{0,1,...,\frac{N}{2}-1\} H(k)=n=0∑2N−1h(n)e−i2N2πkn=n=0∑2N−1x(2n+1)e−i2N2πkn,k∈{0,1,...,2N−1}

而通过复杂的推导（详细推导见图片下方链接）后可以得到用G(k)和H(k)得到X(k)的方式

令WNk=e−i2πkNW^k_N=e^{-i\frac{2\pi k}{N}}WNk=e−iN2πk,则
X(k)=G(k)+WNkH(k),k∈{0,1,...,N2−1}X(k)=G(k−N2)+WNkH(k−N2),k∈{N2,N2+1,...,N−1}X(k)=G(k)+W^k_NH(k) ,k\in\{0,1,...,\frac{N}{2}-1\}\\X(k)=G(k-\frac{N}{2})+W^k_NH(k-\frac{N}{2}) ,k\in\{\frac{N}{2},\frac{N}{2}+1,...,N-1\} X(k)=G(k)+WNkH(k),k∈{0,1,...,2N−1}X(k)=G(k−2N)+WNkH(k−2N),k∈{2N,2N+1,...,N−1}

补充完整就是这样

所以这样，我们就实现了用两个分别为N/2个乘法结果的式子来表示N个乘法结果的方法。而每个N/2个点的计算又可以用同样的方式处理，用两个N/4个点的结果来表示。这样嵌套下去，最终，我们计算这个DFT的算法复杂度就可以降低到nlogN

回到算法的代码实现，因为GPU不适合用递归嵌套的方法来实现FFT,所以这里就要用到蝶形网格的方法将递归展开，具体原理可见第一篇Unity的文章

代码实现：

因为海面的高度公式是一个二维的DFT,所以需要两轮的计算，先横向FFT计算一遍,再纵向FFT一遍。利用蝶形网格计算FFT是按一个一个阶段往前（下）推进（细分）的，而一个N个点的FFT共又m=log2N个阶段。从大佬们推导出的公式可以看出，最后一个阶段的计算和前面m-1个阶段略有不同，所以也单独写成一个cs,一共就是两个用于m-1阶FFT的CS和两个专门计算最后一阶段的CS

FFTHorizontal_CS.hlsl

#include"Waves.hlsli"
//横向FFT计算,只针对第m-1阶段
[numthreads(16, 16, 1)]
void FFTHorizontal(uint3 id : SV_DispatchThreadID)
{int2 idxs = id.xy;idxs.x = floor(id.x / (Ns * 2.0f)) * Ns + id.x % Ns;float angle = 2.0f * PI * (id.x / (Ns * 2.0f));float2 w = float2(cos(angle), sin(angle));float2 x0 = InputRT[idxs].xy;float2 x1 = InputRT[int2(idxs.x + N * 0.5f, idxs.y)].xy;float2 output = x0 + float2(w.x * x1.x - w.y * x1.y, w.x * x1.y + w.y * x1.x);OutputRT[id.xy] = float4(output, 0, 0);
}

FFTHorizontalEnd_CS.hlsl

#include"Waves.hlsli"
//横向FFT最后阶段计算,需要进行特别处理
[numthreads(16, 16, 1)]
void FFTHorizontalEnd(uint3 id : SV_DispatchThreadID)
{int2 idxs = id.xy;idxs.x = floor(id.x / (Ns * 2.0f)) * Ns + id.x % Ns;float angle = 2.0f * PI * (id.x / (Ns * 2.0f));float2 w = float2(cos(angle), sin(angle));/*********不同点***********/w *= -1;/***************************/float2 x0 = InputRT[idxs].xy;float2 x1 = InputRT[int2(idxs.x + N * 0.5f, idxs.y)].xy;float2 output = x0 + float2(w.x * x1.x - w.y * x1.y, w.x * x1.y + w.y * x1.x);/*********不同点***********/int x = id.x - N * 0.5f;output *= ((x + 1) % 2.0f) * 1 + (x % 2.0f) * (-1);/***************************/OutputRT[id.xy] = float4(output, 0, 0);
}

水平方向的两个CS类似，就不贴代码了

然后就是到了C++端的逻辑代码。FFT的计算必须计算2的N次幂，这里的N就是前面初始化时的“纹理阶次”FFTPow，知道了阶次，也就相当于知道了FFT在蝶形网格中的阶段数。我们利用两个循环来一一处理每个阶段的计算。

在每个阶段中，需要计算高度，和水平xz共三个DFT。并根据循环增量是否等于FFTPow来选择计算FFT的CS以区分最后一个阶段的计算。

 //计算FFT////ffttype 1.Horizontal  2.HorizontalEnd  3.Vertical  4.VerticalEndint FFTType=0;//横向FFTfor (int m = 1; m <= FFTPow; m++){//更新Nsm_CBns.ns = (int)pow(2, m - 1);D3D11_MAPPED_SUBRESOURCE mappedData;deviceContext->Map(m_pNsConstantBuffer.Get(), 0, D3D11_MAP_WRITE_DISCARD, 0, &mappedData);memcpy_s(mappedData.pData, sizeof m_CBns, &m_CBns, sizeof m_CBns);deviceContext->Unmap(m_pNsConstantBuffer.Get(), 0);FFTType = m == FFTPow ? 2 : 1;//设置fft类型ComputeFFT(deviceContext, m_pHeightSpectrumRT.GetAddressOf(), FFTType);ComputeFFT(deviceContext, m_pDisplaceXSpectrumRT.GetAddressOf(), FFTType);ComputeFFT(deviceContext, m_pDisplaceZSpectrumRT.GetAddressOf(), FFTType);}//纵向FFTfor (int m = 1; m <= FFTPow; m++){//更新Nsm_CBns.ns = (int)pow(2, m - 1);D3D11_MAPPED_SUBRESOURCE mappedData;deviceContext->Map(m_pNsConstantBuffer.Get(), 0, D3D11_MAP_WRITE_DISCARD, 0, &mappedData);memcpy_s(mappedData.pData, sizeof m_CBns, &m_CBns, sizeof m_CBns);deviceContext->Unmap(m_pNsConstantBuffer.Get(), 0);FFTType = m == FFTPow ? 4 : 3;//设置fft类型ComputeFFT(deviceContext, m_pHeightSpectrumRT.GetAddressOf(), FFTType);ComputeFFT(deviceContext, m_pDisplaceXSpectrumRT.GetAddressOf(), FFTType);ComputeFFT(deviceContext, m_pDisplaceZSpectrumRT.GetAddressOf(), FFTType);}

其中，这里要更新每个阶段FFT计算时HLSL中要用到的参数Ns

然后是具体每个FFT的计算。

这里是把需要计算的纹理作为参数传入，作为Input给到计算着色器来访问，当线程组完成运算得到一个小阶段的输出并存放在OutputRT后，再返回给传入的纹理，从而达到不断更新三个纹理的目的。

void WavesRender::ComputeFFT(ID3D11DeviceContext* deviceContext, ID3D11Texture2D** ppInputTex, int fftType)
{m_pInputRT = *ppInputTex;//设置CS资源//type 1.Horizontal  2.HorizontalEnd  3.Vertical  4.VerticalEndswitch (fftType){case 1:deviceContext->CSSetShader(m_pFFTHorizontalCS.Get(), nullptr, 0);break;case 2:deviceContext->CSSetShader(m_pFFTHorizontalEndCS.Get(), nullptr, 0);break;case 3:deviceContext->CSSetShader(m_pFFTVerticalCS.Get(), nullptr, 0);break;case 4:deviceContext->CSSetShader(m_pFFTVerticalEndCS.Get(), nullptr, 0);break;default:break;}ID3D11UnorderedAccessView* pUAVs[2] = { m_pInputUAV.Get(),m_pOutputUAV.Get() };deviceContext->CSSetUnorderedAccessViews(5, 2, pUAVs, nullptr);//调用线程运算deviceContext->Dispatch(m_CBInitSettings.fftSize / 16, m_CBInitSettings.fftSize / 16, 1);//解除绑定pUAVs[0] = pUAVs[1] = nullptr;deviceContext->CSSetUnorderedAccessViews(5, 2, pUAVs, nullptr);//交换输入输出纹理ID3D11Texture2D* ptmpRT = *ppInputTex;*ppInputTex = m_pOutputRT.Get();m_pOutputRT = ptmpRT;
}

完成了xyz三个纹理横纵两个方向的FFT后，再调用生成偏移纹理的CS

TextureGenerationDisplace_CS.hlsl

#include"Waves.hlsli"
//生成偏移纹理
[numthreads(16, 16, 1)]
void TextureGenerationDisplace(uint3 id : SV_DispatchThreadID)
{float y = length(HeightSpectrumRT[id.xy].xy) / (N * N) * HeightScale; //高度float x = length(DisplaceXSpectrumRT[id.xy].xy) / (N * N) * Lambda; //x轴偏移float z = length(DisplaceZSpectrumRT[id.xy].xy) / (N * N) * Lambda; //z轴偏移HeightSpectrumRT[id.xy] = float4(y, y, y, 0);DisplaceXSpectrumRT[id.xy] = float4(x, x, x, 0);DisplaceZSpectrumRT[id.xy] = float4(z, z, z, 0);DisplaceRT[id.xy] = float4(x, y, z, 0);
}

//计算纹理偏移
//
//设置CS资源
deviceContext->CSSetShader(m_pTextureGenerationDisplaceCS.Get(), nullptr, 0);
ID3D11UnorderedAccessView* pUAVs4[4] = { m_pHeightSpectrumUAV.Get(),m_pDisplaceXSpectrumUAV.Get(),m_pDisplaceZSpectrumUAV.Get(),m_pDisplaceUAV.Get() };
deviceContext->CSSetUnorderedAccessViews(1, 4, pUAVs4, nullptr);
//调用线程运算
deviceContext->Dispatch(m_CBInitSettings.fftSize / 16, m_CBInitSettings.fftSize / 16, 1);
//解除绑定
pUAVs4[0] = pUAVs4[1] = pUAVs4[2] = pUAVs4[3]= nullptr;
deviceContext->CSSetUnorderedAccessViews(1, 4, pUAVs4, nullptr);

我们就得到了xyz方向的偏移纹理（也就是前面很多纹理的那张图中中间灰色的三张）。以及最终的一张偏移纹理（包含xyz方向）

生成法线和泡沫纹理

随后我们可以通过采样当前位置点周围的四个点的，利用差分的方法计算出法线的方向

利用差分计算较为简单，但准确度一般。最好的方式是再次利用DFT计算出该点的空间梯度，再借助几何关系得到准确的法向量

同时还可以使用雅可比行列式计算顶点重合产生的泡沫，具体原理见前面几篇文章。

代码实现：

TextureGenerationNormalBubbles_CS.hlsl

#include"Waves.hlsli"
//生成法线和泡沫纹理
[numthreads(16, 16, 1)]
void TextureGenerationNormalBubbles(uint3 id : SV_DispatchThreadID)
{//计算法线float uintLength = OceanLength / (N - 1.0f); //两点间单位长度//获取当前点，周围4个点的uv坐标uint2 uvX1 = uint2((id.x - 1.0f + N) % N, id.y);uint2 uvX2 = uint2((id.x + 1.0f + N) % N, id.y);uint2 uvZ1 = uint2(id.x, (id.y - 1.0f + N) % N);uint2 uvZ2 = uint2(id.x, (id.y + 1.0f + N) % N);//以当前点为中心，获取周围4个点的偏移值float3 x1D = DisplaceRT[uvX1].xyz; //在x轴 第一个点的偏移值float3 x2D = DisplaceRT[uvX2].xyz; //在x轴 第二个点的偏移值float3 z1D = DisplaceRT[uvZ1].xyz; //在z轴 第一个点的偏移值float3 z2D = DisplaceRT[uvZ2].xyz; //在z轴 第二个点的偏移值//以当前点为原点，构建周围4个点的坐标float3 x1 = float3(x1D.x - uintLength, x1D.yz); //在x轴 第一个点的坐标float3 x2 = float3(x2D.x + uintLength, x2D.yz); //在x轴 第二个点的坐标float3 z1 = float3(z1D.xy, z1D.z - uintLength); //在z轴 第一个点的坐标float3 z2 = float3(z1D.xy, z1D.z + uintLength); //在z轴 第二个点的坐标//计算两个切向量float3 tangentX = x2 - x1;float3 tangentZ = z2 - z1;//计算法线float3 normal = normalize(cross(tangentZ, tangentX));//计算泡沫float3 ddx = x2D - x1D;float3 ddz = z2D - z1D;//雅可比行列式float jacobian = (1.0f + ddx.x) * (1.0f + ddz.z) - ddx.z * ddz.x;jacobian = saturate(max(0, BubblesThreshold - saturate(jacobian)) * BubblesScale);NormalRT[id.xy] = float4(normal, 0);BubblesRT[id.xy] = float4(jacobian, jacobian, jacobian, 0);
}

//生成法线和泡沫纹理
//设置CS资源
deviceContext->CSSetShader(m_pTextureGenerationNormalBubblesCS.Get(), nullptr, 0);
deviceContext->CSSetUnorderedAccessViews(4, 1, m_pDisplaceUAV.GetAddressOf(), nullptr);
deviceContext->CSSetUnorderedAccessViews(7, 1, m_pNormalUAV.GetAddressOf(), nullptr);
deviceContext->CSSetUnorderedAccessViews(8, 1, m_pBubblesUAV.GetAddressOf(), nullptr);
//调用线程运算
deviceContext->Dispatch(m_CBInitSettings.fftSize / 16, m_CBInitSettings.fftSize / 16, 1);
//解除绑定
ID3D11UnorderedAccessView* nullpUAV[1] = { nullptr };
deviceContext->CSSetUnorderedAccessViews(4, 1, nullpUAV, nullptr);
deviceContext->CSSetUnorderedAccessViews(7, 1, nullpUAV, nullptr);
deviceContext->CSSetUnorderedAccessViews(8, 1, nullpUAV, nullptr);m_AccumulateTime = 0.0f;      // 重置累积时间

至此，一次完整的海浪偏移纹理的更新也就完成了，最后记得重置累积时间。

绘制：

绘制就比较简单了。IA阶段设置好顶点，索引缓冲区。再利用xjun项目框架中的特效管理类来设置好管线的资源和状态，Apply应用更改后就可以按索引绘制出波浪了。

void WavesRender::Draw(ID3D11DeviceContext* deviceContext, BasicEffect& effect)
{UINT strides[1] = { sizeof(VertexPosNormalTex) };UINT offsets[1] = { 0 };deviceContext->IASetVertexBuffers(0, 1, m_pVertexBuffer.GetAddressOf(), strides, offsets);deviceContext->IASetIndexBuffer(m_pIndexBuffer.Get(), DXGI_FORMAT_R32_UINT, 0);effect.SetWavesStates(true, 1.0f / m_NumCols, 1.0f / m_NumCols, m_CBInitSettings.m_SpatialStep); // 开启波浪绘制effect.SetMaterial(m_Material);//设置材质effect.SetTextureDiffuse(m_pTextureDiffuse.Get());    //设置波浪纹理effect.SetTextureDisplace(m_pDisplaceSRV.Get());    // 设置位移贴图effect.SetTextureNormal(m_pNormalSRV.Get());   // 设置法线effect.SetTextureBubbles(m_pBubblesSRV.Get());   // 设置泡沫effect.SetWorldMatrix(m_Transform.GetLocalToWorldMatrixXM());effect.Apply(deviceContext);deviceContext->DrawIndexed(m_IndexCount, 0, 0);  //绘制波浪// 解除占用effect.SetTextureDisplace(nullptr);    effect.SetTextureNormal(nullptr);   effect.SetTextureBubbles(nullptr);  effect.SetWavesStates(false);   // 关闭波浪绘制effect.Apply(deviceContext);
}

同样不要忘记解除绑定

以上就是这篇介绍用Direct3D 11实现FFT海面模拟的全部内容了。因为本人数学不好，很多公式的原理和算法理解的还不是很透彻。所以基本上都是直接搬来套用。同时对于Dx的知识还是由很多不熟悉和不理解的。所以写出来的东西也比较浅薄，基本都是自己的一种通俗的理解。

后续完善项目后应该会翻新（或重新）写一篇思路更清晰的介绍。

欢迎指出里面的（一堆）错误，谢谢！