基於快速傅立葉轉換的水波生成(4)

• 前言

終於到了最後一篇，為了寫出這長達四篇的文章，花了我將近1個月的時間...

學習怎麼使用blog，怎麼嵌入程式碼，怎麼把數學符號放上去，讓我數度感到挫折。

廢話就到此為止，開始說明本篇內容。

本篇仍參考 Keith Lantz 所寫的 Ocean simulation part two: using the fast Fourier transform

(以下簡稱 Keith Lantz)文章進行實作。

該文章演示如何利用CPU透過快速傅立葉轉換(以下簡稱FFT)產生水面波浪。

我的目標是修改內容，利用GPU來實現，並在手機上可以有60fps的效能。

• 實作原理分析

回頭觀察水面波浪的方程式和離散傅立葉轉換

可以發現彼此架構一樣，所以可以用快速傅立葉轉換去進行運算。

在 Keith Lantz 裡面，有花了很大的篇幅去證明如何將波浪方程式分解成奇數項和偶數項，

在此我不詳述證明原理，列出幾個要點就好

1.Keith Lantz 將2維的水面波浪方程式展開成兩個1維的水面波浪方程式

也就是將 X 軸 和 Z 軸分開計算，
接下來只要將其中一軸分解出奇數項和偶數項，另外一軸同理可得。

2.為了方便計算，Keith Lantz 將波形長度和取樣點統一



之後分解成奇數項和偶數項



其中，要注意的是 -1^x 這個值。
由於他不再快速傅立葉轉換的公式裡面，因此經過快速傅立葉轉換後，要根據 X 和 Z 的索引值，決定乘上 1 或 -1。

• 程式碼實作

在實作上，為了在手機上能有60fps的效能，我設定了一個前提是假設若風速不會太大，那海浪的變化也不會太大，所以更新也不用太頻繁。

因此，我將波浪的計算分成4個步驟

1. 更新時域頻譜
2. 對X軸座標做FFT運算
3. 對Z軸座標做FFT運算
4. 將計算出的波浪位移和法線值儲存在兩張貼圖上，並Tiling化(或稱seamless)

把這4個步驟分散在4個Frame做批次運算，避免在同一個Frame計算造成計算峰值過高，導致遊戲會有突然會有卡頓的現象發生。

在海浪完成更新前，海浪的波形仍然需要被渲染，因此需要將之前計算出來的波形位移和法線值儲存，作為下一次更新前的海浪渲染資料來源。

由於海浪波形是個週期波，所以可以利用Tiling的方式來去對大範圍的海面做波浪計算。

原本的演算法是只對取樣點 0 到 N-1 之間做計算，並將結果儲存到貼圖內；但貼圖內沒有 N - 1 到 N (也就是 0)之間的波浪資料。因此，若直接拿來做Tiling，會出現一個明顯的斷層

為了解決這個問題，勢必要將取樣點 N-1 到 N (也就是0)的資料一併寫入貼圖內。

假設 N = 128，實作的流程是：

首先，產生1個 129 x 129 個節點的矩形Mesh，然後將X、Z軸的取樣點索引寫入UV值內。

接著將步驟3完成的波浪位移和法向量資料，根據UV的取樣索引點寫入到兩張貼圖儲存。

為了有足夠的空間能夠儲存所有資料，貼圖大小必須上升一個冪次(也就是256 x 256)。

如此，就可以得到一個Tiling化(seamless)的波浪位移與法向量貼圖。

但這個做法存在著一個問題點，Unity Mesh的節點數量上限是 65535。

當取樣超過128時，Mesh 的節點數會超過 65535。

雖然可以用多個Mesh去解決這個問題，

但考量在手機上運作，128 x 128個取樣點已經相當足夠了。

因此我設定取樣點上限為128 x 128來避免這個問題。

• 測試結果

我的測試環境是在 Sony Xperia L3 上運作，並使用 ApowerREC 錄製。

FPS測定是使用 SRDebugger 。

在同時錄製下，仍可維持 fps 60，算是達到我的目標。

但在計算白沫(whitecap)時，我使用 ddx 和 ddy 來計算 jacobian determinant，

似乎在手機環境下，ddx 和 ddy 回傳的數值和 PC 環境比較起來明顯的粗糙。

所以白沫部分會顯得像是馬賽克一樣。

手機上的白沫計算可能還是採用 GPU Gems 2 提到的利用高度門檻

來進行渲染的方式會比較適合。

• 總結

網路上關於使用FFT做波浪生成的相關文章其實並不少，但大多是國外(簡體也有)。

因此作為自己的第一篇文章給有志之士批評指教。

在這裡介紹兩個不錯的Unity專案，同樣是使用FFT生成海浪來做為參考與學習

(我的專案也是參考這兩個專案做成的)

• Phillips-Ocean

也是根據 Keith Lantz 文章的內容所實作的海浪生成，

和我的專案不同的地方是 Phillips-Ocean 是在 CPU 運算(使用 Thread)，而我是在GPU運算。

對於不熟悉 Unity Shader 的工程師，是個很好的學習教材

• Brunetons-Ocean 

使用 JONSWAP Spectrum 作為海浪頻譜，並示範如何混合兩個頻譜輸出，

海洋渲染實現了  Real-time Realistic Ocean Lighting using Seamless Transitions from Geometry to BRDF [1]，提及的方式。

如果要在專案中做出海洋效果，直接參考這個專案比較好。

最後附上我的專案 FastFourierOcean

參考

1.Bruneton, E., Neyret, F., Holzschuch, N: Real-time Realistic Ocean Lighting using Seamless Transitions from Geometry to BRDF. Comput. Graph. Forum 29(2), 487-496 (2010),Blackwell Publishing Ltd.

基於快速傅立葉轉換的水波生成(3)

• 前言

前兩篇文章簡單敘述如何用離散傅立葉轉換做水波生成，但是存在著效能問題。

在我的桌上型主機上，長、寬各取64個取樣點時，勉強還能運行；

但當取樣點來到128個時候，Unity就沒有反應，只能出動檔案管理員來關閉執行

就算做出來，也無法實際運用在遊戲內。

因此，需要使用快速傅立葉轉換來改善效能。

在實際實作之前，本篇文章會講述快速傅立葉的原理，以方便加深對實作程式碼的理解。

• 快速傅立葉轉換

先回到傅立葉轉換的公式

為了方便說明，我假設所有的訊號的是週期波，範圍在 0 ~ 2PI 內。

然後將公式離散化得到

然後展開數學式，根據 n 的值，分成奇數項(n = 1,3,5,...)的和與偶數項(n = 2,4,6,...)的和



接著，移動相位前進半個取樣點(這個點我稱之為反轉相位，因為彼此相位差距PI)，
然後同樣展開成基數項與偶數項的和，比較兩者差異

 

從這裡，可以觀察到兩者的差異只是在奇數項和的正負值。

因此，可以得知當求得某個相位的傅立葉係數，那我可以立刻得知反轉相位的傅立葉係數

大幅減少計算總量，這就是快速傅立葉的原理。


那麼，在實作快速傅立葉轉換需要注意幾個點：

1.離散取樣數必須是2的冪次

只有這樣才能讓每個取樣點，他的反轉相位也是取樣點

2.為了加速找到每個取樣點的反轉相位，會建立一個索引表幫助搜尋

 這個索引表用圖來表示的話會看起來如下圖所示



因為形狀看起來就像蝴蝶，所以又被叫做 Butterfly Lookup Table

• 總結

本篇只是稍微介紹快速傅立葉轉換的原理，沒有打算去實作快速傅立葉轉換。

因此在最後實作快速傅立葉轉換的水渲染時會引用其他高手所公開的程式碼。

因為原理只是將奇數項偶數項分開計算，

所以不管提供的資料內容是什麼，都可以不用修改直接引用。

基於快速傅立葉轉換的水波生成(2)

• 前言

嚴格來說，這篇不是介紹如何用快速傅立葉轉換做水波生成，
而是用離散傅立葉轉換做水波生成。

為了方便理解如何用傅立葉轉換做水波生成，我將文章分為4篇，從簡而繁去做介紹。
想看主題的話可以直接看第4篇。(前提是有人看...)


本篇的文章是參考 Ocean simulation part one: using the discrete Fourier transform
(以下用作者名稱 Keith Lantz 簡稱)這個網頁的說明，製作出來的。

這個做法又稱為 spectrum-based Ocean，
Keith Lantz 根據<Simulating Ocean Water>[1] 的記載使用 Phillip's spectrum，
來模擬海浪波形。並使用二維反向傅立葉轉換將頻譜轉換到每個X-Z平面的浪高，藉此渲染出海洋形狀。

Keith Lantz 利用CPU運算實現了文獻上的做法，

本文則是挑戰修Keith Lantz的做法，改用GPU運算實現，並公開在GitHub上作為自己曾經走過遊戲業的痕跡。

不再廢話了，開始進入正題，

• 演算法分析

<Simulating Ocean Water>[1]定義了在時間t，海平面座標(x,z)的高度為


Keith Lantz 為了方便計算，將 n 和 m 正規化到 0 ~ N-1 和 0 ~ M - 1，並命名為 n '和 m'，
代入方程式得到



現在將焦點放在頻譜函式(註1)




最後，代入 n', m' 得到


就完成了整個演算法所需要的環節。

• Unity 程式碼實現

從演算法內容可以看出，初始頻率函式是個不隨時間變化的常數值，所以先計算每個(n',m')位置的初始頻率值，然後記錄在2維貼圖上，供GPU運算讀取

phillips spectrum 程式碼實現如下

float Phillips(int _n_prime, int _m_prime)
{
    Vector2 k = new Vector2(Mathf.PI * (2 * _n_prime - N) / LenN, Mathf.PI * (2 * _m_prime - M) / LenM);
    float k_length = k.magnitude;
    if (k_length < 0.000001f) return 0.0f;

    float k_length2 = k_length * k_length;
    float k_length4 = k_length2 * k_length2;

    k.Normalize();

    float k_dot_w = Vector2.Dot(k, Wind.normalized);
    float k_dot_w2 = k_dot_w * k_dot_w;

    float w_length = Wind.magnitude;
    float L = w_length * w_length / GRAVITY;
    float L2 = L * L;

    float damping = 0.001f;
    float l2 = L2 * damping * damping;

    //return Amp * Mathf.Exp(-1.0f / (k_length2 * L2)) / k_length4 * k_dot_w2;
    return Amp * Mathf.Exp(-1.0f / (k_length2 * L2)) / k_length4 * k_dot_w2 * Mathf.Exp(-k_length2 * l2);
}

本程式碼是參考 Keith Lantz 所改寫，和上個章節內容有著一些明顯的不同，

Keith Lantz 新增了 Amplitude 參數在他的演算法內，

並且回傳值多乘了  Mathf.Exp(-k_length2 * l2)

參考網路上的一些討論，似乎是為了修正在高頻區域會產生無法收束的狀況而做的一個修正

高斯隨機變數 程式碼實現如下

Vector2 GaussianRandomVariable()
{
    float x1, x2, w;
    do
    {
        x1 = 2.0f * UnityEngine.Random.value - 1.0f;
        x2 = 2.0f * UnityEngine.Random.value - 1.0f;
        w = x1 * x1 + x2 * x2;
    }
    while (w >= 1.0f);

    w = Mathf.Sqrt((-2.0f * Mathf.Log(w)) / w);
    return new Vector2(x1 * w, x2 * w);
}

得到了phillips spectrum 和 高斯隨機變數，初始頻率函式就跟著出來了

Vector2 H_Tilde0(int _n_prime, int _m_prime)
{
    Vector2 r = GaussianRandomVariable();
    return r * Mathf.Sqrt(Phillips(_n_prime, _m_prime) / 2.0f);
}

分母的根號2其實是在計算高斯隨機變數使用的，但為了節省計算，搬到外面和 phillips spectrum 一起處理

波浪的分散關係式

float Dispersion(int _n_prime, int _m_prime)
{
    // w(k) = (gk)^0.5
    float w_0 = 2.0f * Mathf.PI / 200.0f;
    float kx = Mathf.PI * (2 * _n_prime - N) / LenN;
    float kz = Mathf.PI * (2 * _m_prime - N) / LenM;

    return Mathf.Floor(Mathf.Sqrt(GRAVITY * Mathf.Sqrt(kx * kx + kz * kz)) / w_0) * w_0;
    //return Mathf.Sqrt(GRAVITY * Mathf.Sqrt(kx * kx + kz * kz));
}

在這裡也可以看出 Keith Lantz 做了部分改寫，但為了方便比較，
我以 Keith Lantz 為主

最後頻譜函數常數部分的實現方式

每個座標(n',m')都有5個常數值，將這5個常數值儲存在貼圖的5個channel(程式碼用到兩張貼圖)，剩下的就交由GPU在runtime處理。

GPU運算的實現使用shader來處理，程式碼如下

• Unity下的執行結果

專案位置

海水渲染shader是我自己寫的，可以依照專案需求更改

註1.頻譜函數是我自己取的名字，無論是在文獻[1]，或是 Keith Lantz  的網頁中，都沒有給予一個名稱，只說明是回傳一個複數振幅的函式

註2.這個也是我自己取的名字，但他的確是求出一個起始的頻譜函式

參考
[1] SIGGRAPH 2001, Tessendorf J. Simulating ocean water[J]. Simulating nature: realistic and interactive techniques.

基於快速傅立葉轉換的水波生成(1)

作為目前電影業界廣泛採用的海洋波浪生成方案，基於快速傅立葉轉換(Fast Fourier Transform, 之後簡稱 FFT)的水波生成方式其優點是高真實度和豐富的細節。

對於一個發現自己走錯行的遊戲工程師，決定把自己對FFT水水波生成的理解當作不務正業開始墮落的第一篇網誌。

在討論實際做法之前，必須先了解什麼是傅立葉轉換。
先從數學上來說，這個方程式就是傅立葉轉換

... 看得懂有鬼！

但是不要急，將這段方程式用中文來解釋就是：

....聽得懂有鬼！

看不懂也聽不懂怎麼辦？只剩下賣雞排或異世界轉生這兩條路了...

別急！讓我用一個簡單的範例來做說明。
這是一個看起來有點複雜的訊號波形 f(t)



但實際上，卻是可以用簡單的Sin波相加合成


簡單分析Sin波可以得到


因此，範例中的3個Sin波可以分析成



為了方便理解，將波形資料圖形化，取Sin波的振幅為縱軸，角周波數為橫軸，可以得到



到了這裡，回頭看看一開始說明的傅立葉轉換，先不管那複雜的時間積分之類的描述，
我們將範例波形由時間函數f(t)轉換成角周波函數F(w)來表示，而這就是傅立葉轉換在做的事情。

歸納一下上述總結，可以用下圖來解釋



但現實上，不是所有訊號都能像上述範例那麼簡單。
範例是用到3個週期波形所組合，所以可以簡單地被分解。
那什麼是週期波形？顧名思義，就是訊號的形狀呈現週期性(一段時間後，形狀重複)


非週期函數在計算上有其困難性。
因此，在水波生成計算上，會假設訊號都是由週期波形所組合



這叫做傅立葉級數展開。
然後整理一下傅立葉級數展開，可以得到


賽已、摳賽、甲賽  ...
看到那麼多個賽，都想跑廁所了...
還好 Eular 大老師幫我們建立了一個公式



然後，利用這個公式，連立求解得到



代入到傅立葉級數可以得到



將同類項合併



從結果可以觀察到



接下來，可以把整個方程式做個整理




比照之前的角周波函數F(w)



所以，只要求得所有的係數，就完成了傅立葉轉換！

... 回到了最根本的問題 ，係數怎麼求？


在談到如何計算係數之前，要先介紹三角函數的一個特性，正交
不同角周波數的三角函數彼此正交，證明如下



那正交有什麼好處，答案是求係數很方便

舉一個2維向量作範例

(2,5) = a*(1,0) + b*(0,1)， 求 a,b 值

看到題目可以直接知道 a =2. b = 5。
但是我們用數學的方式來求，就使用到內積

a = dot((2,5),(1,0)) = 2*1 + 5*0 = 2
b = dot((2,5),(0,1)) = 2*0 + 5*1 = 5

現在我們回去看看最後求得的時間函數 f(t)


由於複數內積要取共軛，所以可得




在水波生成的實作上，之前提到會假設波形是由週期函數所組合
因此，



這就是離散傅立葉轉換。



寫了那麼多，簡單的介紹傅立葉轉換，
下一篇會開始進入正題介紹如何用傅立葉轉換做水波生成。


註1. 本文對於angular frequency一詞翻譯成角周波，在國內常用角頻或角頻率來說明，但因為想強調其周期性和本身哈日(日文也是使用角周波一詞)，所以就直接使用角周波。
之後的文章也會繼續使用角周波來代替angular frequency