針對這兩篇教程：
http://www.keithlantz.net/2011/10/ocean-simulation-part-one-using-the-discrete-fourier-transform/
http://www.keithlantz.net/2011/11/ocean-simulation-part-two-using-the-fast-fourier-transform/
的一些注解：
一，
原文下面一段
框起部分不通順，似乎應該是 "by letting Lx=N and Lz=M. then becomes,"
從這段之后教程中本應是Lx，Lz的地方就都變成M,N了。
由此看來，作者為了方便使用fft計算heightmap，強制使Lx，Lz等於M,N了。
而在教程的非fft的版本中是沒有Lx，Lz必須等於N，M的限制的，甚至Lx,Lz可以不是整數。
即對於非fft版本，ocean的構造函數
cOcean::cOcean(const int N, const float A, const vector2 w, const float length, const bool geometry)
中N和length可以傳不同的值，但對於fft版本，N和length必須傳相同整數值才正確（N須為2的冪）。
補充：坑爹，原來我用的瀏覽器有問題，導致教程中插在句子中的數學符號顯示不出來，剛才用另一個瀏覽器看，發現紅框中那句確實是：
整篇教程中類似情況不勝枚舉，我居然一直都是在缺符號的情況下看的。。。
二，
教程中h_tilde為NxN，vertices為(N+1)x(N+1)。
vertices的初始化代碼為：
for (int m_prime = 0; m_prime < Nplus1; m_prime++) {
for (int n_prime = 0; n_prime < Nplus1; n_prime++) {
index = m_prime * Nplus1 + n_prime;
htilde0 = hTilde_0( n_prime, m_prime);
...
vertices[index].ox = vertices[index].x = (n_prime - N / 2.0f) * length / N;
vertices[index].oy = vertices[index].y = 0.0f;
vertices[index].oz = vertices[index].z = (m_prime - N / 2.0f) * length / N;
...
}
}
從中可以看出h'(x,z,t)的x和z的初始表達式為：
x=(n'-N/2)*L/N
z=(m'-N/2)*L/N
對於fft版本，由於教程中強制令L=N，所以
x=n'-N/2
z=m'-N/2
（注，后面為了實現“卷浪”，x，z在此基礎上還會進行一定的偏移）。
三，
教程中fft版本的總體計算思路是將二維DFT：
拆成兩個一維DFT來進行計算
以N=4為例，計算過程如下：
其中DFT(4)代表(4)式所定義的DFT，DFT(3)代表(3)式所定義的DFT。
即先逐行dft，再逐列dft。
可見對於N=4的情況，完成heightmap的更新需要計算8個DFT，那么對於一般情況則需要計算2N個DFT。
每個DFT都可以使用FFT來計算，對應教程中如下代碼：
for (int m_prime = 0; m_prime < N; m_prime++) {
fft->fft(h_tilde, h_tilde, 1, m_prime * N);
…
}
for (int n_prime = 0; n_prime < N; n_prime++) {
fft->fft(h_tilde, h_tilde, N, n_prime);
…
}
四，
教程代碼中實現無縫循環（tiling）的代碼非常簡單，只是將第0行（列）的數據填充給第N+1行（列）即可：
for (int m_prime = 0; m_prime < N; m_prime++) {
for (int n_prime = 0; n_prime < N; n_prime++) {
index = m_prime * N + n_prime; // index into h_tilde..
index1 = m_prime * Nplus1 + n_prime; // index into vertices
…
// height
vertices[index1].y = h_tilde[index].a;
…
// for tiling
if (n_prime == 0 && m_prime == 0) {
vertices[index1 + N + Nplus1 * N].y = h_tilde[index].a;
…
}
if (n_prime == 0) {
vertices[index1 + N].y = h_tilde[index].a;
…
}
if (m_prime == 0) {
vertices[index1 + Nplus1 * N].y = h_tilde[index].a;
…
}
}
}
言外之意，h'(x,z,t)在x和z方向上均是以L為周期的周期函數，即：
h'(x+L,z,t)=h(x,z,t)
h'(x,z+L,t)=h(x,z,t)
事實確實如此，驗證如下：
因為N為偶數，所以
所以
h'(x+L,z,t)=h(x,z,t)
h'(x,z+L,t)=h(x,z,t)
五，
evaluateWavesFFT函數中如下代碼：
for (int m_prime = 0; m_prime < N; m_prime++) {
fft->fft(h_tilde, h_tilde, 1, m_prime * N);
…
}
for (int n_prime = 0; n_prime < N; n_prime++) {
fft->fft(h_tilde, h_tilde, N, n_prime);
…
}
int sign;
float signs[] = { 1.0f, -1.0f };
vector3 n;
for (int m_prime = 0; m_prime < N; m_prime++) {
for (int n_prime = 0; n_prime < N; n_prime++) {
index = m_prime * N + n_prime; // index into h_tilde..
index1 = m_prime * Nplus1 + n_prime; // index into vertices
sign = signs[(n_prime + m_prime) & 1];//如果n+m為奇數，則sign=signs[1]=-1；如果n+m為偶數，則sign=signs[0]=1
h_tilde[index] = h_tilde[index] * sign;
…
}
}
可見在計算完2N次FFT后，還有一個符號校正（乘以sign）的步驟。
為什么需要符號校正？
因為代碼中的第一個fft並不是完整計算了h''(x,m',t)，而只是計算了下式中紅框內的部分（省略了一個(-1)^x因子）：
同理，代碼中的第二個fft也並不是完整計算了h'(x,z,t)，而只是計算了下式中紅框內的部分（省略了一個(-1)^z因子）：
所以當fft計算全部執行完后，還需要乘以(-1)^(x+z)。
因為在fft版本中x=n'-N/2，z=m'-N/2，所以(-1)^(x+z)=(-1)^(n'+m'-N)，因為N為偶數，所以(-1)^(n'+m'-N)=(-1)^(n'+m')。
六，
蝶形圖及T函數表實現
這一塊發現教程中有些錯誤，我們自己來推導一下N=4的蝶形圖。
在上式中令N=4，得：
因x=n'-N/2，n'=0,1,2,3，所以x=-2,-1,0,1，並適當利用
變形，得：
據此可畫出蝶形圖：
對上面蝶形圖進行等價優化變形，得：
注：以上優化變形的原理是：因為，所以有T4^0=-T4^(-2)，故可將T4^(-2)由支干移至主干，並將T4^0改為-1。同理，可將T4^(-1)由支干移至主干，並將T4^1改為-1。
容易看出變換后的蝶形圖與原蝶形圖等價。此優化思路源於：https://cnx.org/contents/zmcmahhR@7/Decimation-in-time-DIT-Radix-2，其中的Additional Simplification一節。
此即為N=4時的最終蝶形圖，可以驗證此蝶形圖與教程中給出的N=4蝶形圖是不等價的，即教程中蝶形圖有誤。
再看教程中生成T函數表的代碼：
cFFT::cFFT(unsigned int N) : N(N), reversed(0), T(0), pi2(2 * M_PI) {
…
int pow2 = 1;
T = new complex*[log_2_N]; // prep T
for (int i = 0; i < log_2_N; i++) {
T[i] = new complex[pow2];
for (int j = 0; j < pow2; j++) T[i][j] = t(j, pow2 * 2);
pow2 *= 2;
}
…
}
complex cFFT::t(unsigned int x, unsigned int N) {
return complex(cos(pi2 * x / N), sin(pi2 * x / N));
}
顯然也是有錯誤的，其中
for (int j = 0; j < pow2; j++) T[i][j] = t(j, pow2 * 2);
一句，應改為：
for (int j = 0; j < pow2; j++) T[i][j] = t(j-N/2, pow2 * 2);
因為x=n'-N/2，而非x=n'。
此至，教程及代碼中的疑問基本上就都解決了，下面是移植到unity中的結果：
（為了確認無縫，在z軸方向平鋪了一次）
----補充：
一，
蝶形圖的含義：
首先我們知道蝶形圖右邊任何一個輸出都是左邊所有輸入的線性組合，所以假設我們想知道h'''(-1,m',t)的表達式，首先有
然后再根據蝶形圖讀出A,B,C,D的值。
求A的值，看紅色路徑，沒有經過任何值，所以A=1。
求B的值，看綠色路徑，經過T2^(-2)和-1，所以B=-T2^(-2)。
求C的值，看藍色路徑，經過T4^(-1)，所以C=T4^(-1)。
求D的值，看黃色路徑，經過T2^(-2)，-1，T4^(-1)，所以D=-T2^(-2)*T4^(-1)。
所以得：
對比前面（六）中所得：
因為-T2^(-2)=T2^(-1)，所以結果是一樣的。
二，
比特翻轉
根據蝶形圖實現fft，因為蝶形圖左邊N個輸入數據的是亂序（比如對於N=4而言是0,2,1,3；對於N=8而言是0,4,2,6,1,5,3,7），這個次序可由bit reverse來生成，以N=4為例：
0 (00) -- bit reverse --> (00) 0
1 (01) -- bit reverse --> (10) 2
2 (10) -- bit reverse --> (01) 1
3 (11) -- bit reverse --> (11) 3
參考：http://www.dspguide.com/ch12/2.htm
代碼中實現bit reverse的部分為：
cFFT::cFFT(unsigned int N) : N(N), reversed(0), T(0), pi2(2 * M_PI) {
…
log_2_N = log(N)/log(2);
reversed = new unsigned int[N]; // prep bit reversals
for (int i = 0; i < N; i++) reversed[i] = reverse(i);
…
}
unsigned int cFFT::reverse(unsigned int i) {
unsigned int res = 0;
for (int j = 0; j < log_2_N; j++) {
res = (res << 1) + (i & 1);
i >>= 1;
}
return res;
}
----補充2
此教程為cpu實現fft ocean，且fft部分的實現方法用的是最直觀的方法，所以效率是比較不高的方法。
但即使cpu fft ocean再怎么優化，效果肯定也高不到哪兒去，所以此教程僅作為一個起點。接下來移植到gpu是必須的。等移植完我再另寫一篇日志。
2017-6-25更新：
基於gpu的實現已初步試驗成功，見下帖中的2017-6-25：http://www.cnblogs.com/wantnon/p/6985141.html
2017-6-30更新：
gpu fft海面動畫已實現，見下帖中的2017-6-30: http://www.cnblogs.com/wantnon/p/6985141.html
----
另附兩個重要參考:
http://graphics.ucsd.edu/courses/rendering/2005/jdewall/tessendorf.pdf
https://pdfs.semanticscholar.org/0047/8af7044a7f1350d5ec75ffc7c15b40057051.pdf