fft ocean注解
針對這兩篇教程：
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/
的一些注解：
一，
原文下面一段
框起部分不通順，似乎應(yīng)該是 "by letting Lx=N and Lz=M. then becomes,"
從這段之后教程中本應(yīng)是Lx，Lz的地方就都變成M,N了。
由此看來，作者為了方便使用fft計算heightmap，強(qiáng)制使Lx，Lz等于M,N了。
而在教程的非fft的版本中是沒有Lx，Lz必須等于N，M的限制的，甚至Lx,Lz可以不是整數(shù)。
即對于非fft版本，ocean的構(gòu)造函數(shù)
cOcean::cOcean(const int N, const float A, const vector2 w, const float length, const bool geometry)
中N和length可以傳不同的值，但對于fft版本，N和length必須傳相同整數(shù)值才正確（N須為2的冪）。
補(bǔ)充：坑爹，原來我用的瀏覽器有問題，導(dǎo)致教程中插在句子中的數(shù)學(xué)符號顯示不出來，剛才用另一個瀏覽器看，發(fā)現(xiàn)紅框中那句確實(shí)是：
整篇教程中類似情況不勝枚舉，我居然一直都是在缺符號的情況下看的。。。
?
二，
教程中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的初始表達(dá)式為：
x=(n'-N/2)*L/N
z=(m'-N/2)*L/N
對于fft版本，由于教程中強(qiáng)制令L=N，所以
x=n'-N/2
z=m'-N/2
（注，后面為了實(shí)現(xiàn)“卷浪”，x，z在此基礎(chǔ)上還會進(jìn)行一定的偏移）。
三，
教程中fft版本的總體計算思路是將二維DFT：
拆成兩個一維DFT來進(jìn)行計算
以N=4為例，計算過程如下：
其中DFT(4)代表(4)式所定義的DFT，DFT(3)代表(3)式所定義的DFT。
即先逐行dft，再逐列dft。
可見對于N=4的情況，完成heightmap的更新需要計算8個DFT，那么對于一般情況則需要計算2N個DFT。
每個DFT都可以使用FFT來計算，對應(yīng)教程中如下代碼：
??? 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);
?? ??? ?…
?? ?}
四，
教程代碼中實(shí)現(xiàn)無縫循環(huán)（tiling）的代碼非常簡單，只是將第0行（列）的數(shù)據(jù)填充給第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為周期的周期函數(shù)，即：
h'(x+L,z,t)=h(x,z,t)
h'(x,z+L,t)=h(x,z,t)
事實(shí)確實(shí)如此，驗證如下：
因為N為偶數(shù)，所以
所以
h'(x+L,z,t)=h(x,z,t)
h'(x,z+L,t)=h(x,z,t)
五，
evaluateWavesFFT函數(shù)中如下代碼：
??? 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為奇數(shù)，則sign=signs[1]=-1；如果n+m為偶數(shù)，則sign=signs[0]=1
?? ??? ??? ?h_tilde[index]???? = h_tilde[index] * sign;
?? ??? ??? ?…
?? ??? ?}
?? ?}
可見在計算完2N次FFT后，還有一個符號校正（乘以sign）的步驟。
為什么需要符號校正？
因為代碼中的第一個fft并不是完整計算了h''(x,m',t)，而只是計算了下式中紅框內(nèi)的部分（省略了一個(-1)^x因子）：
同理，代碼中的第二個fft也并不是完整計算了h'(x,z,t)，而只是計算了下式中紅框內(nèi)的部分（省略了一個(-1)^z因子）：
所以當(dāng)fft計算全部執(zhí)行完后，還需要乘以(-1)^(x+z)。
因為在fft版本中x=n'-N/2，z=m'-N/2，所以(-1)^(x+z)=(-1)^(n'+m'-N)，因為N為偶數(shù)，所以(-1)^(n'+m'-N)=(-1)^(n'+m')。
六，
蝶形圖及T函數(shù)表實(shí)現(xiàn)
這一塊發(fā)現(xiàn)教程中有些錯誤，我們自己來推導(dǎo)一下N=4的蝶形圖。
在上式中令N=4，得：
因x=n'-N/2，n'=0,1,2,3，所以x=-2,-1,0,1，并適當(dāng)利用變形，得：
據(jù)此可畫出蝶形圖：
對上面蝶形圖進(jìn)行等價優(yōu)化變形，得：
注：以上優(yōu)化變形的原理是：因為，所以有T4^0=-T4^(-2)，故可將T4^(-2)由支干移至主干，并將T4^0改為-1。同理，可將T4^(-1)由支干移至主干，并將T4^1改為-1。
容易看出變換后的蝶形圖與原蝶形圖等價。此優(yōu)化思路源于：https://cnx.org/contents/zmcmahhR@7/Decimation-in-time-DIT-Radix-2，其中的Additional Simplification一節(jié)。
?
此即為N=4時的最終蝶形圖，可以驗證此蝶形圖與教程中給出的N=4蝶形圖是不等價的，即教程中蝶形圖有誤。
再看教程中生成T函數(shù)表的代碼：
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);
一句，應(yīng)改為：
for (int j = 0; j < pow2; j++) T[i][j] = t(j-N/2, pow2 * 2);
因為x=n'-N/2，而非x=n'。
此至，教程及代碼中的疑問基本上就都解決了，下面是移植到unity中的結(jié)果：
（為了確認(rèn)無縫，在z軸方向平鋪了一次）
----補(bǔ)充：
一，
蝶形圖的含義：
首先我們知道蝶形圖右邊任何一個輸出都是左邊所有輸入的線性組合，所以假設(shè)我們想知道h'''(-1,m',t)的表達(dá)式，首先有
然后再根據(jù)蝶形圖讀出A,B,C,D的值。
求A的值，看紅色路徑，沒有經(jīng)過任何值，所以A=1。
求B的值，看綠色路徑，經(jīng)過T2^(-2)和-1，所以B=-T2^(-2)。
求C的值，看藍(lán)色路徑，經(jīng)過T4^(-1)，所以C=T4^(-1)。
求D的值，看黃色路徑，經(jīng)過T2^(-2)，-1，T4^(-1)，所以D=-T2^(-2)*T4^(-1)。
所以得：
對比前面（六）中所得：
因為-T2^(-2)=T2^(-1)，所以結(jié)果是一樣的。
二，
比特翻轉(zhuǎn)
根據(jù)蝶形圖實(shí)現(xiàn)fft，因為蝶形圖左邊N個輸入數(shù)據(jù)的是亂序（比如對于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
代碼中實(shí)現(xiàn)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;
}
?
----補(bǔ)充2
此教程為cpu實(shí)現(xiàn)fft ocean，且fft部分的實(shí)現(xiàn)方法用的是最直觀的方法，所以效率是比較不高的方法。
但即使cpu fft ocean再怎么優(yōu)化，效果肯定也高不到哪兒去，所以此教程僅作為一個起點(diǎn)。接下來移植到gpu是必須的。等移植完我再另寫一篇日志。
2017-6-25更新：
基于gpu的實(shí)現(xiàn)已初步試驗成功，見下帖中的2017-6-25：http://www.cnblogs.com/wantnon/p/6985141.html
2017-6-30更新：
gpu fft海面動畫已實(shí)現(xiàn)，見下帖中的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
?
總結(jié)
以上是生活随笔為你收集整理的fft ocean注解的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： spring Boot打可执行的jar包
- 下一篇： RabbitMQ安装与初始配置