RealFFT算法实现

FFT最后一步，解释RealFFT的算法原理。
RealFFT输入，对于实序列$x_n(0 \le n \le N - 1)$，以$(x_{2n }, x_{2n + 1})$的形式输入，输出为实序列FFT变换结果的前$N / 2$个点。

比如，以长为16的序列${1, 1, 1, 1, 0, \ldots 0 }$输入ComplexFFT和RealFFT，结果输出分别为

[  4.000000,   0.000000] 
[  3.013670,   2.013670] 
[  1.000000,   2.414214] 
[ -0.248303,   1.248303] 
[  0.000000,   0.000000] 
[  0.834089,  -0.165911] 
[  1.000000,   0.414214] 
[  0.400544,   0.599456] 
[  0.000000,   0.000000] 
[  0.400544,  -0.599456] 
[  1.000000,  -0.414214] 
[  0.834089,   0.165911] 
[  0.000000,   0.000000] 
[ -0.248303,  -1.248303] 
[  1.000000,  -2.414214] 
[  3.013670,  -2.013670] 
=========================
[  4.000000,   0.000000]
[  3.013670,   2.013670]
[  1.000000,   2.414214]
[ -0.248303,   1.248303]
[  0.000000,   0.000000]
[  0.834089,  -0.165911]
[  1.000000,   0.414214]
[  0.400544,   0.599456]

观察可以发现，利用序列的共轭对称性可以补全一些其余的点（有一个点无法补全，即对称中心点）
下面说明算法流程

$$\begin{aligned} X_k & = \sum_{n = 0}^{N - 1}x_nW_N^{kn} \\ & = \sum_{n = 0}^{N / 2 - 1}x_{2n}W_N^{2nk} + x_{2n + 1}W_N^{(2n + 1)k} \\ & = \sum_{n = 0}^{N / 2 - 1}x_{2n}W_{N/2}^{nk} + W_N^k\sum_{n = 0}^{N / 2 - 1}x_{2n + 1}W_{N/2}^{nk} \\ & = F_k + W_N^kG_k \end{aligned}$$

共轭对称证明：

$$\begin{aligned} X_{n-k}^* & = F_{n-k}^* + W_N^{*n-k}G_{n-k}^* \\ & = F_k + W_N^{*-k}G_k = F_k + W_N^kG_k = X_k \end{aligned}$$

之前的$F_k, G_k$已经计算完毕，只需要在上次程序的基础之上计算

$$F_k + W_N^kG_k$$

即可。

实现如下：

ComplexFFT(R, I, 0);
float FR, FI, GR, GI, YR, YI, CYR, CYI, XR, XI, cosr, sinr;
for (int r = 1; r <= MAX_LEN; r++)
{
	if(r == 1)
		FR = R[r], FI = 0.0, GR = I[r], GI = 0.0;
	else
	{
		YR  = R[r], YI = I[r];
		CYR = R[MAX_LEN + 2 - r], CYI = -I[MAX_LEN + 2 - r];
		FR  = (YR + CYR) / 2, FI = (YI + CYI) / 2;
		GR  = (YI - CYI) / 2, GI = (CYR - YR) / 2;
	}
	cosr = cos((r - 1) * PI / MAX_LEN);
	sinr = sin((r - 1) * PI / MAX_LEN);
	XR = FR + cosr * GR - sinr * GI;
	XI = FI + cosr * GI + sinr * GR;
    printf("[%12f, %12f]\n", XR, XI);
}

结果经过验证无误。