RealFFT算法铺垫

继续FFT，上回已经手写FFT了，但是也遗留了一个问题，实际应用中FFT处理的都是实信号，但是为了使用我上回写的ComplexFFT函数，需要给其虚部补0，这一操作不仅浪费空间，而且计算耗时啊。RealFFT算法就是处理实信号的快速算法。

同时，搞信号的重构也让我对实FFT的对称性加深了认识，确实是有用的……

实信号DFT的对称性:

$$\begin{aligned} W_N^k & =e^{-2kj\pi/N} \\ X_k & = \sum_{n = 0}^{N-1}x_nW_N^{nk} \end{aligned}$$

注意这里$x_n$是实数

$$\begin{aligned} X_{N - r} & = \sum_{n = 0}^{N-1}x_nW_N^{(N - r)k} \\ & = \sum_{n = 0}^{N-1}x_nW_N^{-rk}W_N^{Nk} \\ & = \sum_{n = 0}^{N-1}x_nW_N^{-rk} = \sum_{n = 0}^{N-1}[x_nW_N^{rk}]^* = X_r^* \end{aligned}$$

反之：

$$X_{N - r}^* = X_r \quad (1 \le r\le N-1)$$

在$0$这点：

$$X_0 = \sum_{n = 0}^{N-1}x_n$$

直观一点，看个例子，使用ComplexFFT函数，输入

[ 0.000078,    0.000000] => [  67.890442,    0.000000]
[ 1.315378,    0.000000] => [ -13.852791,    1.756794]
[ 7.556053,    0.000000] => [   0.532240,  -14.502794]
[ 4.586502,    0.000000] => [ -10.848876,    0.947389]
[ 5.327672,    0.000000] => [   8.034233,    8.668694]
[ 2.189592,    0.000000] => [  -8.089980,   12.082935]
[ 0.470446,    0.000000] => [ -14.220805,    6.269228]
[ 6.788647,    0.000000] => [   5.620103,    0.964417]
[ 6.792964,    0.000000] => [  -2.237431,    0.000000]
[ 9.346930,    0.000000] => [   5.620103,   -0.964417]
[ 3.835021,    0.000000] => [ -14.220805,   -6.269228]
[ 5.194164,    0.000000] => [  -8.089980,  -12.082935]
[ 8.309653,    0.000000] => [   8.034233,   -8.668694]
[ 0.345721,    0.000000] => [ -10.848876,   -0.947389]
[ 0.534616,    0.000000] => [   0.532240,   14.502794]
[ 5.297002,    0.000000] => [ -13.852791,   -1.756794]

注意，这种对称性不包括直流分量的。
使用对称性还可以做许多优化，包括RealFFT。
这里的铺垫是指：使用一次FFT计算的结果，得到实部，虚部序列的FFT结果，以下面这个为例子:

[ 0.000078,    1.315378] => [  76.081299,   78.423798]
[ 7.556053,    4.586502] => [  -7.566336,   10.051374]
[ 5.327672,    2.189592] => [  -2.170276,  -19.456476]
[ 0.470446,    6.788647] => [  -1.338717,    1.404751]
[ 6.792964,    9.346930] => [ -12.029690,   -0.553211]
[ 3.835021,    5.194164] => [  -1.947378,    5.299602]
[ 8.309653,    0.345721] => [ -19.588976,  -11.740066]
[ 0.534616,    5.297002] => [ -16.009092,    3.573137]
[ 6.711493,    0.076982] => [  11.795471,  -13.576748]
[ 3.834157,    0.668422] => [ -17.982683,   -8.924475]
[ 4.174860,    6.867727] => [ -11.516462,   -1.057013]
[ 5.889766,    9.304365] => [ -32.044506,  -12.575722]
[ 8.461669,    5.269288] => [  12.017742,   -0.259529]
[ 0.919649,    6.539190] => [  16.961304,    5.201901]
[ 4.159994,    7.011906] => [  -0.896533,  -20.641878]
[ 9.103209,    7.621980] => [   6.236089,    5.876601]

根据一个虚数序列$(X, Y)$得到的$(R, I)$，我们是可以得到$(X,0)$和$(Y,0)$的FFT结果的。
推导如下：

$$\begin{aligned} X_k & = \sum_{n = 0}^{N-1}x_nW_N^{nk} = \sum_{n = 0}^{N-1}(f_n+jg_n)W_N^{nk} \\ & = \sum_{n = 0}^{N-1}f_nW_N^{nk}+ j\sum_{n = 0}^{N-1}g_nW_N^{nk} = F_k + jG_k \end{aligned}$$

又$f_n,g_n$是实数序列，满足

$$\begin{aligned} F_{N - k}^* & = F_k \\ G_{N - k}^* & = G_k \end{aligned}$$

那么

$$X_{N - k}^* = F_{N - k}^* - jG_{N - k}^* = F_k - jG_k$$

要想从$X_k$中恢复出$F_k$和$G_k$，可以通过下式：

$$\begin{aligned} F_k & = \frac{1}{2}( X_{N - k}^* + X_k) \\ G_k & = \frac{j}{2}(X_{N - k}^* - X_k) \; 1 \le k \le N - 1 \end{aligned}$$

当$k = 0$时：

$$X_0 = F_0 + jG_0 = \sum_{n = 0}^{N-1}f_n + j\sum_{n = 0}^{N-1}g_n$$

所以：

$$\begin{aligned} F_0 & = R(X_0) \\ G_0 & = I(X_0) \end{aligned}$$

代码实现

ComplexFFT(R, I, 0);	
float FrR, FrI, GrR, GrI, YrR, YrI, CyrR, CyrI;
for (int r = 1; r <= MAX_LEN; r++)
{
	if(r == 1)
		FrR = R[r], FrI = 0.0, GrR = I[r], GrI = 0.0;
	else
	{
		YrR  = R[r], YrI = I[r];
		CyrR = R[MAX_LEN + 2 - r], CyrI = -I[MAX_LEN + 2 - r];
		FrR  = (YrR + CyrR) / 2, FrI = (YrI + CyrI) / 2;
		GrR  = (YrI - CyrI) / 2, GrI = (-YrR + CyrR) / 2;
	}
    printf("[%12f, %12f]\t[%12f, %12f]\n", FrR, FrI, GrR, GrI);
}

程序输出

[   76.081299,     0.000000]	[   78.423798,     0.000000]
[   -0.665123,     2.087387]	[    7.963987,     6.901212]
[   -1.533404,     0.592701]	[  -20.049177,     0.636872]
[    7.811293,    -1.898575]	[    3.303326,     9.150011]
[   -0.005974,    -0.146841]	[   -0.406370,    12.023716]
[  -16.995941,     8.937662]	[   -3.638060,   -15.048564]
[  -15.552719,    -5.341527]	[   -6.398539,     4.036257]
[  -16.995888,     6.248806]	[   -2.675669,    -0.986795]
[   11.795471,     0.000000]	[  -13.576748,     0.000000]
[  -16.995888,    -6.248806]	[   -2.675669,     0.986795]
[  -15.552719,     5.341527]	[   -6.398539,    -4.036257]
[  -16.995941,    -8.937662]	[   -3.638060,    15.048564]
[   -0.005974,     0.146841]	[   -0.406370,   -12.023716]
[    7.811293,     1.898575]	[    3.303326,    -9.150011]
[   -1.533404,    -0.592701]	[  -20.049177,    -0.636872]
[   -0.665123,    -2.087387]	[    7.963987,    -6.901212]

和以$(X,0)$和$(Y,0)$作为输入时的结果做对比

[    76.081299,    0.000000]    [   78.423798,    0.000000]
[    -0.665124,    2.087387]    [    7.963987,    6.901213]
[    -1.533404,    0.592701]    [  -20.049177,    0.636871]
[     7.811293,   -1.898576]    [    3.303326,    9.150011]
[    -0.005974,   -0.146841]    [   -0.406370,   12.023716]
[   -16.995941,    8.937661]    [   -3.638060,  -15.048563]
[   -15.552718,   -5.341526]    [   -6.398538,    4.036257]
[   -16.995888,    6.248806]    [   -2.675670,   -0.986795]
[    11.795471,    0.000000]    [  -13.576748,    0.000000]
[   -16.995888,   -6.248806]    [   -2.675670,    0.986795]
[   -15.552718,    5.341526]    [   -6.398538,   -4.036257]
[   -16.995941,   -8.937661]    [   -3.638060,   15.048563]
[    -0.005974,    0.146841]    [   -0.406370,  -12.023716]
[     7.811293,    1.898576]    [    3.303326,   -9.150011]
[    -1.533404,   -0.592701]    [  -20.049177,   -0.636871]
[    -0.665124,   -2.087387]    [    7.963987,   -6.901213]

基本无误。
熟悉了这部分，下面就可以正式啃RealFFT了。