快速傅立叶变换(FFT) 4.1引言 快速傅立叶变换(FFT)并不是一种新的变换,而是离散傅立叶变换(DFT)的一种快速算法。 DFT的计算在数字信号处理中非常有用。例如在FIR滤波器设计中会遇到从h(n)求H(k)或由H(k)计算h(n),这就要计算DFT;信号的谱分析对通信、图像传输、雷达等都是很重要的,也要计算DFT。因直接计算DFT的计算量与变换区间长度N的平方成正比,当N较大时,计算量太大。 自从1965年图基(J. W. Tukey)和库利(T. W. Coody)在《计算数学》(Math. Computation , Vol. 19, 1965)杂志上发表了著名的《机器计算傅立叶级数的一种算法》论文后,桑德(G. Sand)-图基等快速算法相继出现,又经人们进行改进,很快形成一套高效运算方法,这就是快速傅立叶变换简称FFT(Fast Fourier Transform)。这种算法使DFT的运算效率提高1~2个数量级。 4.2 基2 FFT算法 一、直接计算DFT的问题及改进的途径 设x(n)为N点有限长序列,其DFT正变换为 其反变换(IDFT) x(n)= 二者的差别只在于 考虑x(n)为复数序列的一般情况,每计算一个X(k),需要N次复数乘法以及(N-1)次复数加法。因此,对所有N个k值,共需N2次复数乘法及 下面讨论减少运算工作量的途径。仔细观察DFT的运算就可看出,利用系数 (1) (2) (3) 由此可得: 二、时域抽取法基-2 FFT原理 先设序列点数为N=2M,M为整数。如果不满足这个条件,可以人为地加上若干零值点,使之达到这一要求。这种N为2的整数幂的FFT称基-2 FFT。 设输入序列长度为N=2M (M为正整数) ,将该序列按时间顺序的奇偶分解为越来越短的子序列,称为按时间抽取(DIT )的FFT算法。也称Cooley - Tukey算法。 将N=2M的序列x(n)先按n的奇偶分成以下两组: x(2r+1)=x2(r) 则x(n)的DFT为: = = = = 式中 由(4.2.4)式可看出,一个N点DFT已分解成两个N/2点的DFT,它们按(4.2.4)式又组合成一个N点DFT。 现讨论 式中用了 同理可得: 再考虑 所以前半部分: 后半部分: = 这样,只要求出0到(N/2-1)区间的所有 采用蝶形运算符号表示的图示法,可将上面讨论的分解过程表示于下图中。此图表示N=23=8的情况,其中输出值X(0)到X(3)是由(4.2.7)式给出的,而输出值X(4)到X(7)是由(4.2.8)给出。 既然如此,由于N=2M,因而N/2仍是偶数,可以进一步把每个N/2点子序列再按其奇偶部分分解为两个N/4点的子序列。 x1(2r+1)=x4(r) = = 且 同理, 且 其中, 进一步具体化:N=8=23, =x(0)+x(4) 注意上式中 同理, 三、DIT-FFT算法与直接计算DFT运算量的比较 由上图得,当N=2M时,其运算流图有M级蝶形,每一级都有N/2个蝶形运算构成。每一个蝶形运算需要1次复数乘和2次复数加。所以每一级运算都需要N/2次复数乘和N次复数加。 M级运算总共需要的复数乘法次数为: M级运算总共需要的复数加法次数为: 如直接计算DFT,复数乘法为 当N=1024时,复数乘之比: 显然,N越大,优越性就越明显。 四、按时间抽选的FFT算法的特点: 1、 原位运算 由图4.2.4可以看出,DIT-FFT的运算过程很有规律。N=2M点的FFT共进行M级运算,每级由N/2个蝶形运算组成。同一级中,每个蝶形的两个输入数据只对计算本蝶形有用,而且每个蝶形的输入、输出数据结点又同在一条水平线上,这就意味着计算完一个蝶形后,所得输出数据可立即存入原输入数据所占用的存储单元。这样,经过M级运算后,原来存放输入序列数据的N个存储单元中便依次存放 2、 倒序规律 由图4.2.4看出,按原位计算时,FFT的输出 造成倒序的原因是输入x(n)按标号n的偶奇的不断分组而造成。由于N=2M,所以倒序数可用M位二进制数(nM-1nM-2…n0)2(当N=8=23时,二进制为三位)表示。第一次分组,标为n0。 n为偶数在上半部分,用n0=0表示,n为奇数在下半部分,用n0=1表示;第二次分组,标为n1。偶数部分再分为偶(0)奇(1),奇数部分再分为偶(0)奇(1)…。依次类推,直到M次分组,最后所得二进制倒序数如图示。 下表列出了N=8时以二进制数表示的顺序数和倒序数,由表显而易见,只要将顺序数(n2n1n0)的二进制位倒置,则得对应的二进制倒序值(n0n1n2)。 3、 倒序的实现 设原输入序列x(n)先按自然顺序存入数组A中。例如N=8,A(0),A(1),A(2),A(3),A(4),A(5),A(6),A(7)中依次存放着x(0),x(1),x(2),x(3),x(4),x(5),x(6),x(7)。 顺序数用I表示,I=1~N-2。倒序数用J表示,与I对应分别为4,2,6,1,5,3。当I=J时不需要交换,I<J时调换存放内容。 I=1时,对应的倒序数是4;I=2时,对应的倒序数是2…。倒序数从4到2到…的关系可从表4.2.1得到:每次最高位加1。(注意J用十进制数表示)。如果最高位为0,J直接加N/2,如果最高位为1,则要将最高位归0,次高位加1。但次高位加1时也要判断是否为1或0。程序框图如下图虚线框里所示: 4、 蝶形运算两个输入数据的“距离” 以图4.2.4的8点FFT为例,其输入是倒位序的,输出是自然顺序的。N=2M,共有第一级蝶形运算,第二级蝶形运算,…,第M级蝶形运算。用L表示第某级运算,每个蝶形的两个输入数据的“距离”为B=2L-1。 5、 旋转因子的变化规律 仍观察图4.2.4,每级都有N/2个蝶形,每个蝶形都要乘以因子 L=1,第一级: L=2,第二级: L=3,第三级: 所以对于N=2M的一般情况,第L级的 由于2L=2M2L-M=N2L-M ,所以 所以,第L级的 例如: L=1,第一级:2L-1=1,J=0, L=2,第二级:2L-1=2,J=0,1, L=3,第三级:2L-1=4,J=0,1,2,3, 6、蝶形运算规律 设序列x(n)经时域抽选(倒序)后,存入数组X中。如果蝶形运算的两个输入数据相距B个点,应用原位计算,则蝶形运算可表示成如下形式: 式中 L=1,…,M。 DIT-FFT运算和程序框图如下: 同一旋转因子对应着间隔为2L点的2M-L个蝶形。 五、按频率抽选(DIF)的基2 FFT算法 设序列点数为N=2M ,M为整数。 先把输入按n 的顺序分成前后两半: = = = -1,k为奇数 将 当k取偶数即k=2r时(r=0,1,…,N/2-1), 当k取奇数即k=2r+1时(r=0,1,…,N/2-1), = (4.2.16)用下图表示,N=8。一次分解流图。 由于N=2M,N/2仍是偶数,继续将N/2点DFT分成偶数组和奇数组。图4.2.12表示N=8时二次分解运算流图。 最后完整的分解流图为下图: 这种算法是对 DIF-FFT算法与DIT-FFT算法类似,不同的是DIF-FFT算法输入序列为自然顺序,而输出为倒序排列。另外,蝶形运算略有不同,DIT-FFT蝶形先乘后加,而DIF-FFT蝶形先加后乘。 上述两种FFT的算法流图形式不是唯一的。只要保证各节点所连支路及其传输系数不变,改变输入与输出点以及中间结点的排列顺序,就可以得到其他变形的FFT运算流图。 六、IDFT的高效算法 离散傅立叶反变换 x(n)= 与离散傅立叶正变换( 如果希望直接调用FFT子程序计算IFFT,则可用下面的方法: 由于 x(n)= ∴ x*(n)= x(n)= = 这样可以先将 4.3 进一步减少运算量的措施 研究进一步减少运算量的途径,以程序的复杂度换取计算量的进一步提高 一、 多类蝶形单元运算 由图4.2.4已得出结论,N=2M点FFT共需要MN/2次复数乘法。 由 综上所述,先除去第一、第二两级后,所需复数乘法次数 进一步考虑各级中的无关紧要旋转因子。当L=3时,有两个无关紧要的旋转因子 因此 在基2 FFT程序中,若包含了所有旋转因子,则称该算法为一类蝶形单元运算;若去掉 二、 旋转因子的生成 在FFT运算中,旋转因子 三、 实序列的FFT算法 在实际工作中,数据x(n)一般都是实序列。如果直接按FFT运算流图计算,就是把x(n)看成一个虚部为零的复序列进行计算,这就增加了运算时间。处理这个问题有二种方法,一种是早期提出的用一个N点FFT计算N点实序列的FFT。第二种方法是用N/2点FFT计算一个N点实序列的DFT。 |
|