对于N个数,直接进行DFT计算时需要N2次复数乘法,以及N(N-1)次复数加法,这使得计算效率非常低下,这里介绍一种提高计算效率的算法,即基于时间抽取的快速傅里叶变换。
1.算法原理
设输入序列长度为N=2M(M为整数),将该序列按照时间顺序的奇偶分解为越来越短的子序列,称为基2按时间抽取的FFT算法。也成为Coolkey-Tukey算法。其中即2表示:N=2M,M为整数。若不满足这个条件,可以认为地加上若干零值(加零补长)使其达到N=2m。
2.算法步骤
①分组,变量置换
先将x(n)按n地奇偶分为两组,变量置换:
当n=偶数时,令n=2r;
当n=奇数时,令n=2r+1;
得到
②将上式代入DFT中,有以下结果:
由于旋转因子地性质
可以对X(k)进行变形,可以得到
X1(k),X2(k)只有N/2个点,以N/2为周期;而X(k)却有N个点,以N为周期。要用X1(k),X2(k)表达全部的X(k)的值,还必须要利用旋转因子的周期性:
又根据旋转因子的对称性:
因此可以得出下式:
该式子称之为蝶形运算,蝶形运算流图符号如下所示:
对蝶形运算进行说明:(1)左边两路为输出;(2)右边两路为输出;(3)中间以一个小圆表示加、减运算(右上为相加输出,右下为相减输出)。一次蝶形运算需要1次复乘,2次复加。
采用蝶形运算可以大大减少运算量,有下表做出比较:
例子:求N=23=8点FFT变换
按照基2时间抽取FFT算法,先将N=8点的DFT分解成2个4点DFT,再将4点的DFT分解成2个2点的DFT,再继续分解。
首先分解为两个4点的DFT,有如下分解图:
再分解为2点的DFT,有如下示意图:
最后分解为1点的DFT,有如下示意图:
以上就是整个的过程的推导,在实际代码实现中需要注意第一级蝶形运算时,x(k)打乱了顺序,可以注意到x(1)与x(4),x(3)与x(6)进行了替换,不难发现1的二进制与4的二进制数刚好倒序,因此在进行DFT前需要对原先的数组进行位逆序置换,这里贴上整个基2时间抽取的代码:
//位逆序置换
complex<double>*Rev(complex<double>*Data, int length)
{
int a = 1, num = 0;
while (a <= length)
{
a <<= 1;
num++;
}
vector<int>temp(length);
for (int i = 0; i < length; i++)
{
temp[i] = 0;
}
for (int i = 0; i<length; i++)
{
int m = i;
for (int j = 0; j < num - 1; j++)
{
temp[i] <<= 1;
temp[i] |= (m & 1);
m >>= 1;
}
if (i < temp[i])
{
complex<double> b = Data[i];
Data[i] = Data[temp[i]];
Data[temp[i]] = b;
}
}
return Data;
}
//基2时间抽取FFT,sign=-1时为正变换,sign=1时为逆变换
void dit2(complex<double>*Data, int Log2N, int sign)
{
int i, j, k, step, length;
complex<double> wn, temp, deltawn;
length = 1 << Log2N;
complex<double>*Data = Rev(Data, length);
for (i = 0; i<length; i += 2)
{
temp = Data[i];
Data[i] = Data[i] + Data[i + 1];
Data[i + 1] = temp - Data[i + 1];
}
for (i = 2; i <= Log2N; i++)
{
wn = 1; step = 1 << i; deltawn = complex<double>(cos(2.0*PI / step), sin(sign*2.0*PI / step));
for (j = 0; j<step / 2; j++)
{
for (i = 0; i<length / step; i++)
{
temp = Data[i*step + step / 2 + j] * wn;
Data[i*step + step / 2 + j] = Data[i*step + j] - temp;
Data[i*step + j] = Data[i*step + j] + temp;
}
wn = wn*deltawn;
}
}
if (sign == 1)
for (i = 0; i<length; i++)
Data[i] /= length;
}
参考资料:https://www.bilibili.com/video/BV1W7411c7Kc?p=2(讲的非常仔细地一个视频)