Unity傅里叶变换的应用——真实海平面模拟

By losuffi
周二 16 十月 2018

傅里叶变换的定义描述，网上可以找到很多资源。而且讲得也很好懂。所以，在这篇文章中，不会大篇幅介绍傅里叶变换，假定了读者对傅里叶变换有了一定的了解。本篇的重点是，如何利用傅里叶变换的原理，在Unity游戏引擎中进行实际的应用。

[TOC]

傅里叶变换的算法

就笔者而言，傅里叶变换的算法了解两种。一种是基本的遵循离散傅里叶公式，使用穷举计算所有元素的（DFT）。另一种是在遵循该公式的前提下，利用复数的性质，使用分治法进行迭代的快速傅里叶算法（FFT）。

1. DFT方法

该方法没什么技巧可言。属于看到离散傅里叶公式，直接用代码描述的过程。

离散傅里叶公式如下：

$F(k)=\sum_{n=0}^{N}f(n)*e^{\frac{-j2\Pi nk}{N}}$

依照公式可以写下如下的python的DFT代码：

from cmath import exp,pi
ex=[1.0,1.0,1.0,0]
def dft(f):
    F=[]
    N=len(f)
    for k in range(N):
        s=complex(0)
        for n in range(N):
            s+=f[n]*exp(-2j*pi*n*k/N)
        F.append(s)
    return F
print(dft(ex))

显而易见这是一个时间复杂度为O(n^2)的算法，在理解该方法后，能对傅里叶有更好的认识。

2.FFT方法

该方法使用分治思想，能将复杂度降到O（nlogn）。FFT 算法分为两种

N为输入序列的长度

对N等于2的整数次幂的算法
基2算法
基4算法
......
N不等于2的整数次幂的算法
Winagrad算法
素因子算法
......

算法思路

1. 性质

首先定义：

$W_{N}^{nk}=e^{\frac{-j2\Pi nk}{N}}$

四个可用性质：

周期性
$W_{N}^{nk}=W_{N}^{(N+n)k}=W_{N}^{n(N+k)}$
对称性

$(W_{N}^{nk})^-(-代表共轭)=W_{N}^{-nk}=W_{N}^{(N-n)k}=W_{N}^{n(N-k)}$
可约性
$W_{N}^{nk}=W_{mN}^{mnk}=W_{N/m}^{nk/m}$
特殊点
$\begin{split} W_{N}^{0}&=1 \\W_{N}^{N/2}&=-1 \\W_{N}^{k+N/2}&=-W_{N}^{k} \end{split}$

四种性质的证明，不是太难，利用欧拉公式转换成三角函数就非常明了了。

首先欧拉公式如下：

$e^{jx}=\cos{x}+j\sin{x}$

以周期性为例。

$W_{N}^{nk}=e^{\frac{-j2\Pi nk}{N}}=cos{(\frac{-nk}{N}2\Pi)}+sin(\frac{-nk}{N}2\Pi)$

$\begin{align} W_{N}^{(N+n)k} & = cos{(\frac{-nk-Nk}{N}2\Pi)}+jsin{(\frac{-nk-Nk}{N}2\Pi)} \\\\ &=cos(\frac{-nk}{N}2\Pi-k2\Pi)+jsin{(\frac{-nk}{N}}2\Pi-k2\Pi) \\\\ &=cos(\frac{-nk}{N}2\Pi)+jsin{(\frac{-nk}{N}2\Pi)} \\\\ &=W_{N}^{nk} \end{align}$

其他性质证明大同小异，是来自三角函数的性质

2.思路

利用上述性质，进行迭代，将DFT分解成更小的DFT计算，合并重复计算。

算法有两条路线

时间抽取
频率抽取

按“基数”不同可以分成

基-2FFT算法
基-4FFT算法
混合基算法
分裂基算法

本篇描述是以时间抽取的基-2FFT算法，如果读者对其他算法感兴趣，可以在理解原理后，去查阅相关材料。

3.算法步骤：时间抽取的基-2FFT

保证输入序列的长度为2的整数幂次。如4、128、512等，序列长度不满足时，补空项。
将输入序列按索引分成奇偶两组
$\begin{align} \\\\ 偶序列 f_1(i) &=f(2i) \\\\ 奇序列 f_2(i) &=f(2i+1) \\\\ i&=0,1...,\frac{N}{2}-1 \end{align}$
利用性质，进行分治迭代，一直分奇偶两组，直到序列项数为1(N=1)，无法再分。此时单项的傅里叶变换F(0)=f(n)。

$\begin{align} F(k)&=\sum_{n=0}^{N-1}f(n)W_N^{nk}=\sum_{r=0}^{N/2-1}f(2r)W_N^{2rk}+\sum_{r=0}^{N/2-1}f(2r+1)W_N^{（2r+1）k}\\\\ &=\sum_{r=0}^{N/2-1}f(2r)W_{N/2}^{rk}+W_N^k\sum_{r=0}^{N/2-1}f(2r+1)W_{N/2}^{rk}\\\\ &=F_1(k)+W_N^kF_2(k) \end{align}$

到这一步，已经按照上述三个步骤。已经能写出代码了。示例 Python代码如下：

python from cmath import exp,pi from cmath import log ex=[1.0,1.0,1.0,0] def fft(f,l): N=len(f) if N<=1： return [f[0] for k in range(l)] even=fft(f[0::2],l) #偶数序列部分 odd=fft(f[1::2],l) #奇数序列部分 T=[exp(-2j*pi*k/N)*odd[k] for k in range(N)] return [even[k]+T[k] for k in range(l)] print(fft(ex,4))

但是显而易见，如果仅仅只是这样并没有改变时间复杂度。依旧与穷举法一致，但通过观察，该算法可知道，其实这里面很多数值是重复的。比如当N=1时，需要返回l个f[0]。仅仅是为了占位。观察后，如果我们算法将重复的值只计算一次。那么复杂度将降到O(nlogn)
我们还需要一条性质——周期性：当N=1时，k=2时。依照周期性则F(2)=F(1)

$\begin{align} 例子：N&=4;k=7 \\\\F(7)&=F(7-4)=F(3) \\\\=>F(N+l)&=F(l) \\\\F(k)&=F_1(k)+W_N^kF_2(k) \\\\k&=0,1,...,N-1(上述代码依此实现) \\\\利用周期性，&k可以只取到N/2-1，也可以得到同样的结果 \\\\F(k)&=F_1(k)+W_N^kF_2(k) \\\\W_N^{k+N/2}&=-W_N^k(化成三角函数，证明即可) \\\\F(k+N/2)&=F_1(k)-W_N^kF_2(k) \\\\k&=0,1,2...N/2-1 \end{align}$

如此，每一层的分治，计算量将能够分成前半部，和后半部。而只要求出前半部，后半部也可以表示出来。这样就可以得到整个序列。且，时间复杂度降到了O(nlogn)

该算法，因为其运算方式的流程图像，就像一只蝴蝶。所以又叫做蝶形运算。
该算法的Python代码例子如下：

python from cmath import exp,pi from cmath import log ex=[1.0,1.0,1.0,0] def fft(f): N=len(f) if N<=1： return f #每次减去一半长度，到单一序列，长度正好为1 even=fft(f[0::2]) #偶数序列部分 odd=fft(f[1::2]) #奇数序列部分 T=[exp(-2j*pi*k/N)*odd[k] for k in range(N//2)] #只需要一半 return [even[k]+T[k] for k in range(N//2)]+[even[k]-T[k] for k in range(N//2)] print(fft(ex))

4.逆变换IFFT

在计算出FFT后，推导IFFT就比较简单了。因为FFT和IFFT在计算公式上有很大的相似性

$\begin{align} & DFT: F(k)=\sum_{n=0}^{N-1}f(n)W_N^{nk} \\\\& IDFT:f(n)=\frac{1}{N}\sum_{k=0}^{N-1}F(k)W_N^{-nk} \end{align}$

显而易见，逆变换与正变换，在运算上只有两点不同，一是系数1/N,二是旋转因子为共轭变换。

$\begin{align} &依照对称性，可知 W_N^{-nk}=(W_N^{nk})^{-} \\\\&复数乘法已知 (A+Bi)●（C+Di)=(X+Yi) \\\\&那么如何得到（X-Yi)呢？即（X+Yi）的共轭 \\\\&显而易见，依照复数乘法性质。即（A-Bi）和（C-Di）相乘即可。 \\\\=>& f(n)=\frac{1}{N}\sum_{k=0}^{N-1}F(k)W_N^{-nk} \\\\=>& {f(n)}^-=\frac{1}{N}\sum_{k=0}^{N-1}{F(k)}^-W_N^{nk} \\\\=>& f(n)=\frac{1}{N}{\sum_{k=0}^{N-1}{fft(F(k)}^-)}^- \end{align}$

如上所示。只需要在输入傅里叶序列进行一次共轭。再进行fft，结果乘上系数以及再次共轭。即可求出序列的逆变换。

示例Python代码如下。承接上述fft代码:

def ifft(F):
    N=len(F)
    F=[k.conjugate() for k in F] #对输入的傅里叶变换序列 进行共轭运算
    F=fft(F) #利用已经写好的fft函数进行运算
    F=[(k.conjugate())/N for k in F] #得到结果
    return F
F=fft(ex)
print(F)
print(ifft(F))

Unity的傅里叶变换实现

到这一步，快速傅里叶的算法思路，就已经介绍完了。接下来就是，将这套算法移植到Unity中。使用GPU编程，对图像进行傅里叶变换。

ComputeShader

本篇的傅里叶实现，核心是使用ComputeShader，进行实时的傅里叶变换。CPU的傅里叶实现，基本上与Python的实现差不多，只需要将python的例子，移植到C#语言上就没有问题了。但是实时的话，还是利用GPU来运算效率会更高。

Unity的ComputeShader的基本使用，在本篇文章中假定了读者对其有了一定了解，不再做详细的介绍。不熟悉的读者可查阅相关关键词。

ComputeShader本身语法与Shader区别不大。不过ComputeShader比起Shader，设计上主要是针对并行计算。

CS的关键性质

CS的代码执行，首先分线程组。C#中对应API为ComputeShader.Dispatch(int i,int x,int y, int z)。x,y,z可由调用者，自行分配。
每一组中可并行执行多个线程，在CS里使用[numthread(x,y,z)]来分配线程组中的线程数。
在同一个线程组中的线程，可以抽象的认为是并行的，所以线程的调用顺序是不可知的，即6号线程可能先执行再执行1号线程。在Unity中可定义的线程数的极限是1024个
同一个线程组有办法做到共享内存，内存大小为32KB
函数中可以调用线程等待函数，即线程运行到等待函数后，被挂起，在同一个线程组中的线程都运行到该处后，才能继续执行下去

这五点，对于ComputeShader特别关键。本篇介绍得比较粗糙，如果读者不懂，必须要查阅资料弄明白，才能理解接下来的示例。

二维傅里叶变换

未完，待更！

tags: 算法

Losuffi's Blog