社区导航

 
快捷导航
  • 首页
  • 论坛
  • 查看新帖
  • 最新回复
  • 社区活动
  • 联系管理员
  • 消灭零回复
  • E金币兑换
  • 干货
搜索
查看: 3130|回复: 1

[资料分享] 【转】 数字信号处理 教案 ---- 4、快速傅立叶变换(FFT)

[复制链接]

1万

TA的帖子

33

TA的资源

版主

Rank: 6Rank: 6

发表于 2013-7-8 14:54:56 | 显示全部楼层 |阅读模式
第四章 快速傅立叶变换(FFT)
41引言
快速傅立叶变换(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个数量级。
42 2 FFT算法
一、直接计算DFT的问题及改进的途径
设x(n)为N点有限长序列,其DFT正变换为
= , k=0,1,…,N-1
其反变换(IDFT)
x(n)= ,n=0,1,…,N-1
二者的差别只在于的指数符号不同,以及差一个常数乘因子1/N,因而下面我们只讨论DFT正变换的运算量,反变换的运算量是完全相同的。
考虑x(n)为复数序列的一般情况,每计算一个X(k),需要N次复数乘法以及(N-1)次复数加法。因此,对所有N个k值,共需N2次复数乘法及
N(N-1)次复数加法运算。所以直接计算DFT,乘法次数和加法次数都是和N2成正比的,当N很大时,运算量是很可观的,因而需要改进对DFT的计算方法,以减少运算次数。
       下面讨论减少运算工作量的途径。仔细观察DFT的运算就可看出,利用系数以下固有特性,就可减小DFT的运算量:
(1)的对称性 ()*=
(2)的周期性 ==
(3)的可约性 ==
由此可得:===-1,=-。这样利用这些特性,可以将长序列的DFT分解为短序列的DFT。快速傅立叶变换算法正是基于这样的思路而发展起来的。它的算法基本上可以分成两大类:时域抽取法FFT(Decimation-In-Time FFT,简称DIT-FFT)和频域抽取法FFT(Decimation-In-Time FFT,简称DIF-FFT)。
二、时域抽取法基-2 FFT原理
先设序列点数为N=2M,M为整数。如果不满足这个条件,可以人为地加上若干零值点,使之达到这一要求。这种N为2的整数幂的FFT称基-2 FFT。
设输入序列长度为N=2M (M为正整数) ,将该序列按时间顺序的奇偶分解为越来越短的子序列,称为按时间抽取(DIT )的FFT算法。也称Cooley - Tukey算法。
将N=2M的序列x(n)先按n的奇偶分成以下两组:
x(2r)=x1(r)      r=0,1,…,N/2-1
x(2r+1)=x2(r)
则x(n)的DFT为:
=DFT[x(n)]= ,k=0,1,…,N-1
=+
=+
=+
=+
=+,k=0,1,…,N-1 (4.2.4)
式中分别是x1(r)及x2(r)的N/2点DFT:
      ==,k=0,1,…,N/2-1
==,k=0,1,…,N/2-1
由(4.2.4)式可看出,一个N点DFT已分解成两个N/2点的DFT,它们按(4.2.4)式又组合成一个N点DFT。
现讨论的周期性:
===
式中用了=
同理可得:=
再考虑的以下性质:
       ==-
所以前半部分:
   =+,k=0,1,…,N/2-1
后半部分:
      =+
=,(k=0,1,…,N/2-1)
这样,只要求出0到(N/2-1)区间的所有值,即可求出0到(N-1)区间内的所有值,这就大大节省了运算。
=+,k=0,1,…,    (4.2.7)
=,k=0,1,…, (4.2.8)
                     
采用蝶形运算符号表示的图示法,可将上面讨论的分解过程表示于下图中。此图表示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) = x3(r)      r=0,1,…,N/4-1
x1(2r+1)=x4(r)
=+
    =+
    =+,k=0,1,…,N/4-1
=-,k=0,1,…,N/4-1
同理,也可进行同样的分解:
      =+,k=0,1,…,N/4-1
=-,k=0,1,…,N/4-1
其中,==
           ==
进一步具体化:N=8=23,
     =
         =x(0)+x(4) ,k=0,1
    = x(0)+x(4)
        = x(0)+x(4) = x(0)-x(4)
注意上式中=
同理,= x(1)+x(5)
          = x(1)-x(5)
三、DIT-FFT算法与直接计算DFT运算量的比较
由上图得,当N=2M时,其运算流图有M级蝶形,每一级都有N/2个蝶形运算构成。每一个蝶形运算需要1次复数乘和2次复数加。所以每一级运算都需要N/2次复数乘和N次复数加。
M级运算总共需要的复数乘法次数为:
      
M级运算总共需要的复数加法次数为:
      
如直接计算DFT,复数乘法为次,复数加法次。
当N=1024时,复数乘之比:,这样使运算效率提高200多倍。下图为FFT算法和直接DFT算法所需运算量与计算点数N的关系曲线:
显然,N越大,优越性就越明显。
四、按时间抽选的FFT算法的特点:
1、        原位运算
由图4.2.4可以看出,DIT-FFT的运算过程很有规律。N=2M点的FFT共进行M级运算,每级由N/2个蝶形运算组成。同一级中,每个蝶形的两个输入数据只对计算本蝶形有用,而且每个蝶形的输入、输出数据结点又同在一条水平线上,这就意味着计算完一个蝶形后,所得输出数据可立即存入原输入数据所占用的存储单元。这样,经过M级运算后,原来存放输入序列数据的N个存储单元中便依次存放的N个值。这种利用同一存储单元存储蝶形计算输入、输出数据的方法称为原位计算。原位计算可节省大量内存,从而使设备成本降低。
2、        倒序规律
由图4.2.4看出,按原位计算时,FFT的输出是按正常顺序排列在存储单元中,即按,…,的顺序排列,但是这时输入x(n)却不是按自然顺序存储的,而是按x(0),x(4),…,x(7)的顺序存入存储单元,看起来好象是“混乱无序”的,实际上是有规律的,我们称之为倒序。
        造成倒序的原因是输入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个蝶形,每个蝶形都要乘以因子,称其为旋转因子,p称为旋转因子的指数。第L级蝶形运算(从左向右数),共有2L-1个
L=1,第一级:=,J=0
L=2,第二级:=,J=0,1
L=3,第三级:=,J=0,1,2,3
所以对于N=2M的一般情况,第L级的为:
,J=0,1,…,2L-1-1
由于2L=2M2L-M=N2L-M ,所以 ===
所以,第L级的为:,J=0,1,…,2L-1-1
例如:
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个点,应用原位计算,则蝶形运算可表示成如下形式:
   
   
式中;J=0,1,…,2L-1-1;
L=1,…,M。
DIT-FFT运算和程序框图如下:
         
同一旋转因子对应着间隔为2L点的2M-L个蝶形。
五、按频率抽选(DIF)的基2 FFT算法
设序列点数为N=2M ,M为整数。
,k=0,1,…,N-1
先把输入按n 的顺序分成前后两半:
,k=0,1,…,N-1
=+
    =+
    =,k=0,1,…,N-1
=,k=0,1,…,N-1
=(-1)k=   1,k为偶数
                  -1,k为奇数
分解成偶数组与奇数组,
当k取偶数即k=2r时(r=0,1,…,N/2-1),
==
当k取奇数即k=2r+1时(r=0,1,…,N/2-1),
=
=
         n=0,1,…,N/2-1
=,r=0,1,…,N/2-1
=        (4.2.16)
之间可用下图所示的蝶形运算符号表示:
(4.2.16)用下图表示,N=8。一次分解流图。
由于N=2M,N/2仍是偶数,继续将N/2点DFT分成偶数组和奇数组。图4.2.12表示N=8时二次分解运算流图。
   
最后完整的分解流图为下图:
   
    这种算法是对进行奇偶抽取分解的结果,所以称之为频域抽取法FFT。
    DIF-FFT算法与DIT-FFT算法类似,不同的是DIF-FFT算法输入序列为自然顺序,而输出为倒序排列。另外,蝶形运算略有不同,DIT-FFT蝶形先乘后加,而DIF-FFT蝶形先加后乘。
    上述两种FFT的算法流图形式不是唯一的。只要保证各节点所连支路及其传输系数不变,改变输入与输出点以及中间结点的排列顺序,就可以得到其他变形的FFT运算流图。
六、IDFT的高效算法
   离散傅立叶反变换
x(n)= ,n=0,1,…,N-1
与离散傅立叶正变换(= , k=0,1,…,N-1)相比,只要将DFT中的系数改变为,最后乘以1/N,就是IDFT的运算公式。流图输入为,输出为x(n)。因此原来的DIT-FFT改为IFFT后称为DIF-IFFT更合适;原来的DIF-FFT改为IFFT后称为DIT-IFFT更合适。
下图是由DIF-FFT运算流图改成的IFFT运算流图(DIT-IFFT):
    在实际中,有时为了防止运算过程中发生溢出,将1/N分配到每一级蝶形运算中。由于1/N=,所以每级的每个蝶形输出支路均有一相乘因子1/2。如下图示:
  
如果希望直接调用FFT子程序计算IFFT,则可用下面的方法:
由于    x(n)=
∴      x*(n)=
x(n)=
   ={DFT[]}*
这样可以先将取共轭,然后直接调用FFT子程序进行DFT运算,最后取共轭并乘以1/N得到序列x(n)。
43 进一步减少运算量的措施
研究进一步减少运算量的途径,以程序的复杂度换取计算量的进一步提高
一、   多类蝶形单元运算
由图4.2.4已得出结论,N=2M点FFT共需要MN/2次复数乘法。
=,J=0,1,…,2L-1-1得,当L=1时,只有一种旋转因子,所以第一级不需要乘法运算。当L=2时,共有两个旋转因子:=1和=-j,因此第二级也不需要乘法运算。在DFT中,称其值为±1和±j的旋转因子为无关紧要的旋转因子。
    综上所述,先除去第一、第二两级后,所需复数乘法次数
进一步考虑各级中的无关紧要旋转因子。当L=3时,有两个无关紧要的旋转因子。因为同一旋转因子对应着2M-L=N/2L个蝶形运算,所以,第三级共有2N/23=N/4个蝶形不需要复数乘法运算。依次类推,当L≥3时,第L级的2个无关紧要的旋转因子减少复数乘法的次数为2N/2L=N/2L-1。这样,L=3至L=M共减少复数乘法次数为
因此
        =
在基2 FFT程序中,若包含了所有旋转因子,则称该算法为一类蝶形单元运算;若去掉=±1的旋转因子,则称之为二类蝶形单元运算;若再去掉=±j的旋转因子,则称为三类蝶形单元运算;若再处理=,则称之为四类蝶形运算。我们将后三种运算称为多类蝶形单元运算。显然蝶形单元越多,编程就越复杂,但当N较大时,乘法运算的减少量是相当可观的。例如,N=4096时,三类蝶形单元运算的乘法次数为一类蝶形单元运算的75%。
二、   旋转因子的生成
在FFT运算中,旋转因子
=,求余弦和正弦函数值的计算量很大,所以编程时,一种方法是在每级运算中直接产生,另一种方法是在FFT程序开始前预先计算好,存放在数组中,作为旋转因子表,在程序执行过程中直接查表得到。
三、   实序列的FFT算法
在实际工作中,数据x(n)一般都是实序列。如果直接按FFT运算流图计算,就是把x(n)看成一个虚部为零的复序列进行计算,这就增加了运算时间。处理这个问题有二种方法,一种是早期提出的用一个N点FFT计算N点实序列的FFT。第二种方法是用N/2点FFT计算一个N点实序列的DFT。





此帖出自DSP论坛


回复

使用道具 举报

1262

TA的帖子

3

TA的资源

五彩晶圆(初级)

Rank: 7Rank: 7Rank: 7

荣誉会员勋章

发表于 2013-7-8 18:53:29 | 显示全部楼层
不错,学习一下,大热的天没法到处跑了


回复

使用道具 举报

您需要登录后才可以回帖 登录 | 注册

本版积分规则