快速傅里叶变换

Veröffentlicht am | in

感谢傅里叶同学做出的杰出贡献

什么是傅里叶变换？

请看这里

一些先修的概念

复数

就是$\sqrt{-1}$啦，我们写作$i$

单位根

这是整个变换的核心，也就是满足$z^n=1$的所有复数
我们定义$\omega_n=e^{\frac{2\pi i}{n}}$
那么所有的n次单位根就是
$$\omega_n=e^{\frac{2\pi i k}{n}},k=0,1,2\cdots n-1$$
这些单位根为啥长这样呢？
仔细想想发现这些单位根构成复平面上的单位圆的内接正$n$边形

欧拉公式

计算机没法计算复数，我们需要用到欧拉公式：
$$e^{i\theta}=cos(\theta)+isin(\theta)$$

离散傅里叶变换

DFT

定义一个多项式$A(x)$为$A(x)=\sum_{i=0}^{n-1}a_ix^i$
那么其系数向量$(a_0,a_1\cdots a_{n-1})$的离散傅里叶变换为：
将$n$次的单位根向量作为多项式的点值带入多项式计算出来的向量，换句话就是说，系数向量的DFT为，$(A(\omega_n^0),A(\omega_n^1)\cdots , A(\omega_n^{n-1}))$

IDFT

IDFT是逆变换，实际上就是解方程的过程
我们不难发现DFT的过程实际上是下面这个矩阵运算：
$$
\begin{bmatrix}
(\omega_{n}^0)^0 & (\omega_{n}^0)^1 & \cdots & (\omega_{n}^0)^{n-1}\\
(\omega_{n}^1)^0 & (\omega_{n}^1)^1 & \cdots & (\omega_{n}^1)^{n-1}\\
\vdots & \vdots & \ddots & \vdots \\
(\omega_{n}^{n-1})^0 & (\omega_{n}^{n-1})^1 & \cdots & (\omega_{n}^{n-1})^{n-1}
\end{bmatrix}
\begin{bmatrix}
a_0\\
a_1\\
\vdots\\
a_{n-1}
\end{bmatrix}
=
\begin{bmatrix}
A(\omega_n^0)\\
A(\omega_n^1)\\
\vdots\\
A(\omega_n^{n-1})
\end{bmatrix}
$$
我们令系数矩阵为$V$，为了解决这个线性方程，我们必须找到$V$的逆。
我们定义$D$，其中$D$中的每个元素都是$V$中对应位置的共轭。
然后我们考察一下$E=D·V$
仔细计算后发现$E$竟然是$nI_n$，其中$I_n$表示的是$n$阶单位矩阵
这样一来我们就发现$V$的逆矩阵$V^{-1}=\frac{1}{n}D$
所以逆变换和正变换在运算上只是一个$n$的倍数和共轭的关系，下面讨论的都是正变换，逆变换类似。

Cooley-Tukey算法

Cooley-Tukey算法是FFT的一种，也是比较好懂的一种。

算法思想

算法的思想和CDQ分治很类似，都是通过预先处理好的一半来求出另一半的答案
算法将$A$的每个分量的运算，按照指数奇偶分类：
$A(\omega_n^k)=\sum_{i=0}^{n-1}a_i\omega_n^{ki}=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_n^{2ki}+\omega_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_n^{2ki}$
酱紫有什么好处呢？
我们知道$\omega_n^{2k}=e^{\frac{2\pi ki}{n}·2}=e^{\frac{2\pi ki}{\frac{n}{2}}}=\omega_{\frac{n}{2}}^{k}$
再把这个放回和式：
$$A(\omega_n^k)=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{\frac{n}{2}}^{ki}+\omega_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{\frac{n}{2}}^{ki}$$
另外由于单位根刚好将单位圆分成$n$边形
所以在原点的另一侧
$\omega_n^{k+\frac{n}{2}}=-\omega_n^{k}$
刚好转了180度
说了这么多，到底这些有什么用呢
我们可以发现：
对于$k<\frac{n}{2}$
$$A(\omega_n^k)=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{\frac{n}{2}}^{ki}+\omega_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{\frac{n}{2}}^{ki}$$
$$A(\omega_n^{k+\frac{n}{2}})=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{\frac{n}{2}}^{ki}-\omega_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{\frac{n}{2}}^{ki}$$
这样一来如果我们处理好了$A(\omega_n^k)$，那么我们就能$O(1)$求出$A(\omega_n^{k+\frac{n}{2}})$
这样递归下去后，算法的复杂度就变为了$O(nlogn)$

递归实现方法

代码不是我写的，代码是Miskcoo写的，看这里（不是我懒，是我写出来的没他好。。。。
这里$n$必须扩展到2的幂

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
void fft(int n, complex<double>* buffer, int offset, int step, complex<double>* epsilon)
{
	if(n == 1) return;
	int m = n >> 1;
	fft(m, buffer, offset, step << 1, epsilon);
	fft(m, buffer, offset + step, step << 1, epsilon);
	for(int k = 0; k != m; ++k)
	{
		int pos = 2 * step * k;
		temp[k] = buffer[pos + offset] + epsilon[k * step] * buffer[pos + offset + step];
		temp[k + m] = buffer[pos + offset] - epsilon[k * step] * buffer[pos + offset + step];
	}
	for(int i = 0; i != n; ++i)
		buffer[i * step + offset] = temp[i];
}

epsilon是事先打好的表

1
2
3
4
5
6
7
8
9
void init_epsilon(int n)
{
	double pi = acos(-1);
	for(int i = 0; i != n; ++i)
	{
		epsilon[i] = complex<double>(cos(2.0 * pi * i / n), sin(2.0 * pi * i / n)); 
		arti_epsilon[i] = conj(epsilon[i]);
	}
}

迭代实现方法

考察一下整个计算过程，你就会发现，如果将计算的下标都做二进制翻转：
比如讲0100变成0010
整个计算就是挨着挨着算的，从0算到1，从0,1算到2,3，从0,1,2,3算到4,5,6,7。。。
所以我们只需要先整理好计算的顺序，然后就可以迭代计算出来了：
这回是我自己写的了。。
预处理单位根

1
2
3
4
5
6
void initUnitRoot(int n) {
	for (int i = 0; i < n; i++) {
		unitRoot[i] = std::complex<double>(cos(2.0 * pi * i / n), sin(2.0 * pi * i / n));
		artiUnitRoot[i] = std::conj(unitRoot[i]);
	}
}

二进制翻转

1
2
3
4
5
6
void bitReverse(int n, std::complex<double> *x) {
	for(int i = 0, j = 0; i != n; ++i) {
		if(i > j) std::swap(x[i], x[j]);
		for(int l = n >> 1; (j ^= l) < l; l >>= 1);
	}
}

变换

1
2
3
4
5
6
7
8
9
10
11
12
13
//FFT use transform(n, x, unitRoot)
//IFFT use transform(n, x, artiUnitRoot)
void transform(int n, std::complex<double> *x, std::complex<double> *w) {
	bitReverse(n, x);
	for(int i = 2; i <= n; i <<= 1) 
		for(int j = 0, m = i >> 1; j < n; j += i) 
			for(int k = 0; k != m; ++k) {
				std::complex<double> z = x[j + m + k] * w[n / i * k];
				x[j + m + k] = x[j + k] - z;
				x[j + k] += z;
			}
}