快速傅里叶变换（FFT）——按时间抽取DIT的基

【1】前言
- 1、DIF计算量
- 2、利用性质改善
【2】公式推导
- 1、N 到 2*N/2
- - a、分解原序列
  - b、分解后的DFT变换
  - c、一系列化简操作之后
  - d、蝶形信号流
  - e、计算量总结
- 2、N/2 到 2*N/4
- - a、分解X2(k)序列
  - b、蝶形信号流（2列）
- 3、N/4 到 2*N/8
- - a、蝶形信号流（3列）
【3】公式总结
【4】特点以及程序框架讲解
- 1、原址运算
- 2、倒位序规律
- 3、蝶形运算两节点的距离
- 4、WN^r的确定
- 5、程序框架
【5】代码实现

【1】前言

1、DIF计算量

更加清楚地了解计算步骤：

观察可知：
1、一次复数乘法需用四次实数乘法和二次实数加法；
2、一次复数加法需二次实数加法
3、整个 DFT 运算总共需要 4N^2 次实数乘法和 2N*(2N—1)次实数加法。
总结：
直接计算 DFT,乘法次数和加法次数都是和 N^2 成正比的。

2、利用性质改善

利用这些性质，将较大的N分解为若干个较小的N然后进行运算。

【2】公式推导

1、N 到 2*N/2

a、分解原序列

b、分解后的DFT变换

c、一系列化简操作之后

d、蝶形信号流

以N=8的序列为例：

e、计算量总结

因而通过第一步分解后，总共需要（N^2/2)+(N/2)=N(N+l )/2约等于 N^2/2 次复数乘法和 N( N/2-1 )+N = N^2/2 次复数加法。由此可见，通过这样分解后的运算工作量差不多节省了一半。

2、N/2 到 2*N/4

a、分解X2(k)序列

b、蝶形信号流（2列）

3、N/4 到 2*N/8

a、蝶形信号流（3列）

【3】公式总结

由于乘法的运算量较大，我们从乘法角度来探讨一下，DFT和FFT的运算量。
设N=2^M；有M列的蝶形信号运算。

从乘法角度：DFT需要N^2,FFT需要N*lbN;

当N=2048时，这一比值为372.4,即直接计算DFT的运算量是FFT运算量的372.4倍。
当点数N越大时，FFT的优点更为明显。

【4】特点以及程序框架讲解

1、原址运算

计算之后，将新的X(k)覆盖原本的X(k)。
注意：是将同一行的X进行覆盖（后面的列覆盖前面的列），不同行之间是没有覆盖关系的。
所以，最后只需要N个存储单元。（N个数据,N行）

2、倒位序规律

输出X(k),序列正常。
输入序列不正常。
原因：X(n)按照标号n的奇偶而不断分组。
例子：

步骤流程：
I+1，最低位+1，向左进位。
J在二进制最高位+1，逢2向右进位。
由此可以从当前的倒序值计算求得下一个倒序值。

观察变址处理，可以发现，只有当J>I时，才将X(I)和X(J)存储内容进行互换。

3、蝶形运算两节点的距离

输入为倒位序，输出为正位序，N=2^ M,在第m级运算，两个节点间的距离为2^(m-1);

4、WN^r的确定

r的变换规律：
1、运算两个节点中第一个节点标号为k，表示为M位的二进制数。
2、将此二进制数乘以2^(M-m)，相当于左移M-m位，把右边空出，此数位r的二进制数。

5、程序框架

【5】代码实现

没有验证代码的正确性，只是按照上面的流程图进行叙述。

#define PI 3.14159//数位倒读
int rev(int i, int m) {//i=0~2^m,m为二进制位数 int j = 0;while (m > 0) {j += (i & 0x01) * (0x01 << (m - 1));//j+=(i%2)*mypow(2,m-1);i >>= 1;//i/=2m -= 1;}return j;
}//快速傅里叶变换
//输入x(n)、N
//输出X(k)
void fft(const float real_in[], const float imag_in[], float real_out[], float imag_out[],int N)
{//【1】获取Mint M = log2(N);  //【2】倒序for (int i = 0;i < N;i++) {//数位倒读 int j;j = rev(i, M);real_out[j] = real_in[i];imag_out[j] = imag_in[i];}//【3】for (int m = 1;m <= M;m++){int B = 2 ^ m - 1;for (int J = 0;J <= B - 1;J++){int P = 2 ^ (M - m) * J;for (int k = J;k <= N - 1;k++){float tmpr1, tmpi1, tmpr2, tmpi2;//临时变量float theta = -2 * PI * P / N ;tmpr1 = real_out[k];tmpi1 = imag_out[k];tmpr2 = cos(theta) * real_out[k + B] - sin(theta) * imag_out[k + B];tmpi2 = cos(theta) * imag_out[k + B] + sin(theta) * real_out[k + B];real_out[k] = tmpr1 + tmpr2;imag_out[k] = tmpi1 + tmpi2;real_out[k + B] = tmpr1 - tmpr2;imag_out[k + B] = tmpi1 - tmpi2;}}}
}

参考资料：

《数字信号处理第三版.刘顺兰版》