音频采样中的FFT算法

参考博文

https://blog.csdn.net/ggn_2015/article/details/68922404
https://blog.csdn.net/jgj123321/article/details/96438301

参考此博文中对多项式乘法系数的求解DFT变换的推导过程，做进一步的分析说明

1.DFT多项式乘法

$B(x)$ = a[0] + a[1] $\times$ $x^1$ + a[2] $\times$ $x^2$+a[3] $\times$ $x^3$ + a[4]$\times$ $x^4$ + a[5] $\times$$x^5$+ a[6] $\times$ $x^6$ + a[7] $\times$ $x^7$
$C(x)$ = b[0] + b[1] $\times$ $x^1$ + b[2] $\times$ $x^2$+b[3] $\times$ $x^3$ + b[4]$\times$ $x^4$ + b[5] $\times$$x^5$+ b[6] $\times$ $x^6$ + b[7] $\times$ $x^7$
首先我们已知系数可以求得任意的点值；
那么如何用点值反求系数呢？为什么要用点值反求系数呢？
计算$B(x)$$\times$$C(x)$的表达式的复杂程度：O(n$\times$m)，就是每一个系数都要进行相乘来获取方程不同幂的系数。用点值法可以减少乘法的计算复杂度，请看如何用分治的思想求解

可将系数看作是一组方程的解，学过线性代数的都清楚$Ax=b$；
x就是所求的系数解，$b$为$f(x)$对应的值，要求所有系数有唯一解的条件就是
$r(A:b)$ = 系数$a$的个数$n$， 只需要构造一个矩阵$A$的秩为$n$就可以得到唯一的以组系数解，可以看出
$$\left[\begin{array}{ccccc}1& 1& 1& 1& 1\\ x_1^1& x_2^1& x_3^1& \cdots & x_n^1\\ x_1^2& x_2^2& x_3^2& \cdots & x_n^2\\ \cdots & \cdots & \cdots & \cdots & x_n^k\\ x_1^n& x_2^n& x_3^n& \cdots & x_n^n\end{array}\right]$$

其中$x_1$、$x_2$、$x_3$、$\cdots$、$x_n$为不相等的取值，根据范德蒙行列式的性质，该行列式的值为非零值，构成的矩阵可逆。$r(A:b)$= 系数$a$的个数$n$。

如何构造$n$个$x$的值，成了当前需要思考的问题，并且$x$取不同值的时候，尽可能的减少与系数的相乘次数，因此就需要一组有规律的$x$的值来尽可能减少求值次数。

$$e^{2\pi ki/n}=w_n^k$$其中$(k=1、2、3、\cdots 、n)$
将用$w_n^k$一组复数平面内的有规律的$x$的值求解。$w_n^k$有何规律？
$$w_n^k=w_{\frac{n}{2}}^{\frac{k}{2}}$$
欧拉公式$e^{ix}=cosx+i\times sinx$
$$w_n^{k+\frac{n}{2}}={\mathrm{e}}^{\frac{\pi k\mathrm{i}}{\mathrm{n}}}\times {\mathrm{e}}^{\pi \mathrm{i}}=-w_n^k$$

可以看出$n$个点的取值不同，获取的$w_n^k$不一样
并且按理我们可以通过$w_n^1$求得$w_n^{1+\frac{n}{2}}$，可以通过$w_{\frac{n}{2}}^1$求出$w_n^2$

$B(x)$ = a[0] + a[1] $\times$ $x^1$ + a[2] $\times$ $x^2$+a[3] $\times$ $x^3$ + a[4]$\times$ $x^4$ + a[5] $\times$$x^5$+ a[6] $\times$ $x^6$ + a[7] $\times$ $x^7$转换成
$A_0(x)$ = a[0] + a[2] $\times$ $x$ + a[4]$\times$ $x^2$ + a[6] $\times$ $x^3$
$A_1(x)$ = a[1] + a[3] $\times$ $x^1$ + a[5]$\times$ $x^2$ + a[7] $\times$$x^3$

请看如何求解
$B(x)$= $A_0(x^2)$+$x\times A_1(x^2)$
$B(w_8{}^1)$ = $A_0(w_4{}^1)$+$w_8{}^1\times A_1(w_4{}^1)$（$w_n^k=w_{\frac{n}{2}}^{\frac{k}{2}}$性质）
$B(w_8{}^5)$ = $A_0(w_4{}^1)$-$w_8{}^1\times A_1(w_4{}^1)$（$w_n^{k+\frac{n}{2}}={\mathrm{e}}^{\frac{\pi k\mathrm{i}}{\mathrm{n}}}\times {\mathrm{e}}^{\pi \mathrm{i}}=-w_n^k$性质）

对应只需要求解$B(w_8^1)$、$B(w_8^2)$、$B(w_8^3)$、$B(w_8^4)$即可
那么$A_0(w_4^1)$、$A_0(w_4^2)$、$A_0(w_4^3)$、$A_0(w_4^4)$又该如何求解呢？

因为$A_0(x)$也可以一样将系数分为奇偶项拆分

$A_0(x)=A_{00}(x^2)+x\times A_{01}(x^2)$ 得到第二层

$A_0(w_4^1)=A_{00}(w_2^1)+w_4^1\times A_{01}(w_2^1)$
同理只需要求解$A_0(w_4^1)、A_0(w_4^2)$
$A_{00}(x)=A_{000}(x^2)+x\times A_{001}(x^2)$ 得到第三层
$A_{00}(w_2^1)=A_{000}(w_1^1)+w_2^1\times A_{001}(w_1^1)$
通过欧拉函数将复数平面内不同二次项函数值的求解过程做了递归分析
这样就可以使得我们求解乘法减少为$nlogn$次，并计算出n个不同的点值

void fft_example( COMPLEX *a,int n,int dft)
{
     for(int i=0;i< n;i++)
	 {
         if(i < reverseVal[i])
			swap(A_ref[i],A_ref[reverseVal[i]]); //系数反转完毕
		 //系数翻转后的默认为a[0] a[4] a[2] a[6] a[1] [5] a[3] a[7]
		 //通过计算得到了A000,A001,A010,A011,A100,A101,A110,A111
		 //通过计算得到了A00(1/2),A00(1),A01(1/2),A01(1),A10(1/2),A10(1),A11(1/2),A(1)
		 //通过计算得到了A0(1/4),A0(2/4),A0(3/4),A0(1),A1(1/4),A1(2/4),A1(3/4),A1(1)
		 //通过计算得到了A(1/8).....
          for(int step=1;step <n;step<<=1)
		  {
             COMPLEX Wn= exp(COMPLEX(0,dft*PI/step)); //计算合并求y值的奇数项x取值;
             for(int j= 0;j<n;j+=step<<1)   //划分的n中的所有小组;
			 {
				 COMPLEX Cn(1,0);
				 for(int k=j;k<j+step;k++)
				 {
					 COMPLEX x = a[k];
					 COMPLEX y = Cn*a[k+step];   //这个就是wn的二分定理；
					 a[k] = x +y;
					 a[k+step] = x-y;     //蝴蝶算法的信息合并;
					 Cn *= Wn;   //每一个小组计算一次后需要放大
				 }
			 }
		  }
	 }
	 if(dft == -1)
	{
		for(int i=0;i<n;i++) a[i]/=n;
	}
}

void FFT(FFT_HANDLE *hfft, COMPLEX *TD2FD)
{
    int i,j,k,butterfly,p;

    int power = NumberOfBits(hfft->count); //对应采样点数对应最大的分治次数step
    /*蝶形运算*/
    for(k = 0; k < power; k++)   //k对应step层;
    {
        for(j = 0; j < 1<<k; j++)
        {
            butterfly = 1 << (power-k);   //step为0的时候就是3-0=3;  butterfly =8
            p = j * butterfly;            //p=0  
            int s = p + butterfly / 2; //s=4 ;
            for(i = 0; i < butterfly/2; i++)
            {
                COMPLEX t = TD2FD[i + p] + TD2FD[i + s];
                TD2FD[i + s] = (TD2FD[i + p] - TD2FD[i + s]) * hfft->wt[i*(1<<k)]; 
                TD2FD[i + p] = t;
                //a[0]与a[4]，a[1]与a[5]，a[2]与a[6] 当step为1的时候
            }
        }
    }
    /*重新排序*/
	for (k = 0; k < hfft->count; k++) 
    {
		int r = BitReverise(k, power);
		if (r > k) 
        {
            COMPLEX t = TD2FD[k];
            TD2FD[k] = TD2FD[r];
            TD2FD[r] = t;
		}
	}
}  //最初没有进行翻转，最终得出的结果为A(1/8)，A(5/8)，A(2/8)、、、、需要翻转。