FFT与游戏开发（一）

FFT与游戏开发（一）

傅里叶变换

傅里叶变换是信号与系统这门课中很重要的内容，它的功能是把信号从时域变成频域。

在计算机中，用到最多的离散傅里叶变换（DFT）。

这里推荐一些参考资料和在线内容：

1. 《Understanding Digital Signal Processing》
2. Discrete Time Signals and Systems

DFT

先来上公式：

[X(m) = sum_{n=0}^{N-1}x(n)e^{-j2pi nm/N} ]

说明：

1. x(n)代表时域输入，n从0到N-1，它是根据一个采样频率(f_s)从连续信号中采样得到的。
2. X(m)代表频域输出，m从0到N-1，相邻两个相邻频域输出之间的频率差为(f_s/N)。

DFT的时间复杂度是(O(N^2))，如果数据量特别大，这个复杂度是难以接受的，因此FFT应运而生。

FFT

这里介绍一种Radix-2 FFT，它的时间复杂度为(O(Nlog N))。它要求(N=2^k)，k为正整数。

推导过程

1. 把DFT中一个频率输出的计算分为奇数项和偶数项两部分
   [egin{aligned} X(m) & = sum_{n=0}^{(N/2)-1}x(2n)e^{-j2pi (2n)m/N} + sum_{n=0}^{(N/2)-1}x(2n+1)e^{-j2pi (2n+1)m/N} \ & = sum_{n=0}^{(N/2)-1}x(2n)e^{-j2pi (2n)m/N} + e^{-j2pi m/N} sum_{n=0}^{(N/2)-1}x(2n+1)e^{-j2pi (2n)m/N} end{aligned} ]

2. 设(W_N = e^{-j2pi /N})
   [X(m) = sum_{n=0}^{(N/2)-1}x(2n)W_{N/2}^{nm} + W_N^m sum_{n=0}^{(N/2)-1}x(2n+1)W_{N/2}^{nm} ]

3. 对于(m' geq N/2)的情况来说，可用(m' = m + N/2)进行带入
   [egin{aligned} X(m') &= X(m + N/2) \ &= sum_{n=0}^{(N/2)-1}x(2n)W_{N/2}^{n(m + N/2)} + W_N^{m+N/2} sum_{n=0}^{(N/2)-1}x(2n+1)W_{N/2}^{n(m+N/2)} \ end{aligned} ]
   1. 而
      [egin{aligned} W_{N/2}^{n(m+N/2)} &= W_{N/2}^{nm} W^{nN/2}_{N/2} \ &= W_{N/2}^{nm} W^n \ &= W_{N/2}^{nm} e^{-j2pi n} \ &= W_{N/2}^{nm} end{aligned} egin{aligned} W_N^{m+N/2} &= W_N^mW_N^{N/2} \ &= W_N^m e^{-jpi} \ &= -W_N^m end{aligned} ]

   2. 带入得
      [egin{aligned} X(m') &= X(m + N/2) \ &= sum_{n=0}^{(N/2)-1}x(2n)W_{N/2}^{nm} - W_N^m sum_{n=0}^{(N/2)-1}x(2n+1)W_{N/2}^{nm} end{aligned} ]
      因此，对于(m' geq N/2)，无需重新计算，只需要改变左侧DFT计算中的符号即可
   3. 简化公式得
      [0 leq m < N/2 egin{aligned} X(m) &= A(m) + W_N^mB(m) \ X(m+N/2) &= A(m) - W_N^mB(m) \ A(m) &= sum_{n=0}^{(N/2)-1}x(2n)W_{N/2}^{nm} \ B(m) &= sum_{n=0}^{(N/2)-1}x(2n+1)W_{N/2}^{nm} end{aligned} ]
      其中A(m)，B(m)为偶数列和奇数列分别进行FFT的结果的数组。

蝶形结构

实际上，一个DFT可以分解为分别把源数据的奇数项和偶数项抽出来分别进行DFT，而N/2尺寸的DFT又可以分解成N/4尺寸的DFT，直到输出数据数量为1。

递归版本FFT实现

可以得到递归版本的FFT实现

def RecursiveFFT(src:List[float])->List[complex]:
	n = len(src)
	if n == 1:
		return [complex(a) for a in src]
	wn = Euler(-2 * pi / n)
	w = complex(1)
	evenList = src[::2]
	oddList = src[1::2]
	evenRet = RecursiveFFT(evenList)
	oddRet = RecursiveFFT(oddList)
	ret = [None] * n
	for k in range(n//2):
		t = w * oddRet[k]
		ret[k] = evenRet[k] + t
		ret[k+n//2] = evenRet[k] - t
		w = w * wn
	return ret

bitwise-reversal

从蝶形结构可以看出，最左侧的输入顺序是乱序的，新的顺序可以通过bitwise-reversal获得，代码展示如下：

def ReverseBit(num:int, bitCount:int)->int:
	binary = bin(num)
	reverse = binary[-1:1:-1]
	reverse = reverse + '0' * (bitCount - len(reverse))
	return int(reverse, 2)

循环版本FFT实现

def IterativeFFT(src:List[float])->List[complex]:
	n = len(src)
	bitCount = int(log2(n))
	rev:List[float] = [src[ReverseBit(i, bitCount)] for i in range(n)]
	# 对每一层
	for s in range(1, bitCount+1):
		# 每一层进行FFT的输入数量（尺寸）
		m = 2 ** s
		# 单位W
		wm = Euler(-2 * pi / m)
		# 每层中每个需要进行FFT的初始索引
		for k in range(0, n, m):
			w = complex(1)
			# 对每个FFT进行处理
			for j in range(0, m//2):
				u = rev[k+j]
				t = w * rev[k+j+m//2]
				rev[k+j] = u + t
				rev[k+j+m//2] = u - t
				w = w * wm
	return rev