第4章快速傅里叶变换(FFT)1

1、第第4章章 快速傅里叶变换快速傅里叶变换(FFT) n4.1 引言引言n4.2 基基2FFT算法算法n4.3 进一步减少运算量的措施进一步减少运算量的措施n4.4 其他快速算法其他快速算法4.1 引言引言n DFT是信号分析与处理中的一种是信号分析与处理中的一种重要变换。因直接计算重要变换。因直接计算DFT的计算量与的计算量与变换区间长度变换区间长度N的平方成正比，当的平方成正比，当N较大较大时，计算量太大，所以在快速傅里叶变时，计算量太大，所以在快速傅里叶变换换(简称简称FFT)出现以前，直接用出现以前，直接用DFT算法算法进行谱分析和信号的实时处理是不切实进行谱分析和信号的实时处理是不切实

2、际的。直到际的。直到1965年发现了年发现了DFT的一种快的一种快速算法以后，情况才发生了根本的变化。速算法以后，情况才发生了根本的变化。 4.2 基基2FFT算法算法n 4.2.1 直接计算直接计算DFT的特点及减少运算量的的特点及减少运算量的基本途径基本途径 n 长度为长度为N的有限长序列的有限长序列x(n)的的DFT为为n 考虑考虑x(n)为复数序列的一般情况，对为复数序列的一般情况，对某一个某一个k值，直接计算值，直接计算X(k)值需要值需要N次复数乘法、次复数乘法、(N-1)次复数加法。次复数加法。 10( )( ),0,1,1NknNnX kx n WkNn 如前所述，如前所述，N

3、点点DFT的复乘次数的复乘次数等于等于N2。显然，把。显然，把N点点DFT分解为几个较分解为几个较短的短的DFT，可使乘法次数大大减少。另，可使乘法次数大大减少。另外，旋转因子外，旋转因子WmN具有明显的周期性和具有明显的周期性和对称性。其周期性表现为对称性。其周期性表现为22()jm lNjmm lNmNNNNWeeW其对称性表现为2mN mN mmNNNNNmmNNWWWWWW 4.2.2 时域抽取法基时域抽取法基2FFT基本原理基本原理 FFT算法基本上分为两大类：时域算法基本上分为两大类：时域抽取法抽取法FFT(Decimation In Time FFT,简简称称 D I T - F

4、 F T ) 和 频 域 抽 取 法和 频 域 抽 取 法FFT(Decimation In Frequency FFT,简简称称DIFFFT)。下面先介绍。下面先介绍DIFFFT算算法。法。设序列设序列x(n)的长度为的长度为N，且满足，且满足2 ,MNM为自然数 按按n的奇偶把的奇偶把x(n)分解为两个分解为两个N/2点的子序列点的子序列12( )(2 ),0,1,12( )(21),0,1,12Nx rxrrNx rxrrn 则x(n)的DFT为/2 1/2 12(21)00/2 1/2 121200( )( )( )(2 )(21)( )( )knknNNnnNNkrkrNNrrNNk

5、krNNrrX kx n Wx n Wxr WxrWx rWx r W由于222222/2jkrNjkrkrkrNNNWeeW所以 /2 1/2 11/22/21200( )( )( )( )( )NNkrkkrkNNNNrrX kx r WWx r WX kW Xkn 其中X1(k)和X2(k)分别为x1(r)和x2(r)的N/2点DFT，即 /2 111/210/2 122/220( )( )( )( )( )( )NkrNrNkrNrX kx r WDFT x rXkx r WDFT x r 由于X1(k)和X2(k)均以N/2为周期，且 ，所以X(k)又可表示为2NkkNNWW 121

6、2( )( )( )0,1,12()( )( )0,1,122kNkNNX kX kW XkkNNX kX kW Xkk图4.2.1 蝶形运算符号 CABA BCA BC图4.2.2 N点DFT的一次时域抽取分解图(N=8) N/2点DFTWN0N/2点DFTWN1WN2WN3x(0)X1(0)x(2)x(4)x(6)x(1)x(3)x(5)x(7)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)n 与第一次分解相同，将x1(r)按奇偶分解成两个N/4长的子序列x3(l)和x4(l)，即3241( )(2 )

7、,0,1,1( )(21)4x lxlNlx lxl那么，X1(k)又可表示为 /4 1/4 12(21)11/21/200/4 1/4 13/4/24/4003/24( )(2 )(21)( )( )( )( ),0,1,/21NNklklNNiiNNklkklNNNiikNX kxl WxlWx l WWx l Wx kWXk kN式中 /4 133/430/4 144/440( )( )( )( )( )( )NklNiNklNix kx l WDFT x lxkxl WDFT xl 同理，由X3(k)和X4(k)的周期性和Wm N/2的对称性 Wk+N/4 N/2=-Wk N/2 最后

8、得到： 13/2413/24( )( )( ),0,1,/41(/4)( )( )kNkNX kXkWXkkNX kNXkWXkn 用同样的方法可计算出25/2625/26( )( )( ),0,1,/41(/4)( )kNkNXkXkWXkkNXkNX kWXk其中 /4 155/450/4 166/4605262( )( )( )( )( )( )( )(2 ),0,1,/41( )(21)NklNiNklNiXkx l WDFT x lXkx l WDFT x lx lxllNx lxl图4.2.3 N点DFT的第二次时域抽取分解图(N=8) N/4点DFTWN12WN12WN0WN1W

9、N2WN3X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)x(0)X3(0)X3(1)X4(0)X4(1)x(4)x(2)x(6)x(1)x(5)x(3)x(7)N/4点DFTN/4点DFTN/4点DFTWN02WN02 图4.2.4 N点DITFFT运算流图(N=8) WN0WN1WN2WN3WN0WN2WN0WN2WN0WN0WN0WN0 x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(1)A(2)A(3

10、)A(4)A(5)A(6)A(7)A(0)A(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)n 4.2.3 DITFFT算法与直接计算DFT运算量的比较n 每一级运算都需要N/2次复数乘和N次复数加(每个蝶形需要两次复数加法)。所以，M级运算总共需要的复数乘次数为22(2)log22(2)logMANNCMNCN MNN复数加次数为 例如，N=210=1024时221048576204.8(/2)log5120NNN图4.2.5 FFT算法与直接计算DFT所需乘法次数的比较曲线 由于N=2L，因此N/2仍为偶

11、数，可以依照上面方法进一步把每个N/2点子序列，再按输入n的奇偶分解为两个N/4点的子序列，按这种方法不断划分下去，直到最后剩下的是2点DFT，两点DFT实际上只是加减运算加减运算。蝶形结即蝶式计算结构也即为蝶式信号流图上面频域频域中前/后半部分表示式可以用蝶形信号流图表示。X1(k)X2(k)kNW)()(21kXWkXkN)()(21kXWkXkN作图要素：作图要素：(1)左边两路为输入左边两路为输入(2)右边两路为输出右边两路为输出(3)中间以一个小圆表示加、中间以一个小圆表示加、减运算（右上路为相加输减运算（右上路为相加输出、右下路为相减输出出、右下路为相减输出)(4)如果在某一支路上

12、信号需要进行如果在某一支路上信号需要进行相乘运算，则在该支路上标以箭头，相乘运算，则在该支路上标以箭头，将相乘的系数标在箭头旁。将相乘的系数标在箭头旁。(5)当支路上没有箭头及系当支路上没有箭头及系数时，则该支路的传输比数时，则该支路的传输比为为1。例子：求 N=23=8点FFT变换 （1）先按）先按N=8-N/2=4，做做4点的点的DFT：)()() 2/()()()(2121kXWkXNkXkXWkXkXkNkN12/, 0Nk先将N=8 DFT分解成2个4点DFT：可知：时域上：x(0),x(2),x(4),x(6)为偶子序列 x(1),x(3),x(5),x(7)为奇子序列 频域上：X

13、(0)X(3)，由X(k)给出 X(4)X(7),由X(k+N/2)给出(a)比较N=8点直接DFT与分解2个4点DFT的FFT运算量N=8点的直接DFT的计算量为：N2次(64次）复数相乘,N(N-1)次(8(8-1)=56次)复数相加.共计120次。（b)求 一个蝶形结需要的运算量要运算一个蝶形结，需要一次乘法 ， 两次加法。)(2kXWkN)()() 2/()()()(2121kXWkXNkXkXWkXkXkNkN（c）分解为两个N/2=4点的DFT的运算量分解2个N/2点（4点）的DFT：) 12()2()rxDFTrxDFTkX（偶数其复数相乘为复数相加为奇数其复数相乘为复数相加为2

14、2）（N）（122NN22）（N）（122NN次。共计为次，（次，复数相加为复数相乘为5624) 123222NNN(d)用2个4点来求N=8点的FFT所需的运算量再将N/2点（4点）合成N点(8点)DFT时，需要进行N/2个蝶形运算)()() 2/()()()(2121kXWkXNkXkXWkXkXkNkN12/, 0Nk还需N/2次（4次）乘法 及 次(8次）加法运算。22N)(2kXWkN计算量。分解就可以节省约一半次。看出仅做一次共计为次复数加法次复数乘法）共需求点所以68,322) 12(,36221(22)(8222NNNNNNNNNkXFFT(e)将N=8点分解成2个4点的DFT

15、的信号流图 4点DFTx(0)x(2)x(4)x(6)4点DFTx(1)x(3)x(5)x(7)08W18W28W38WX(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)偶数序列奇数序列3821282118210821) 3 () 3 (3) 2() 2(2) 1 () 1 (1) 0() 0(0WXXXWXXXWXXXWXXX）（）（）（）（如：X(4)X(7)同学们自已写x1(r)x2(r)(2)N/2(4点)-N/4(2点)FFT(a)先将先将4点分解成点分解成2点的点的DFT：n因为4点DFT

16、还是比较麻烦，所以再继续分解。n若将N/2(4点)子序列按奇/偶分解成两个N/4点（2点）子序列。即对将x1(r)和x2(r)分解成奇、偶两个N/4点（2点）点的子序列。奇序列、偶序列、) 6() 2() 4() 0(: )(1xxxxrx奇序列、偶序列、同理：)7() 3 () 5 () 1 (: )(2xxxxrx1014012()()2(4131，），在此（）奇序列（）偶序列若设：LNLLXLxLxLx1014012()()2(5252，），在此（）奇序列（）偶序列同理：LNLLXLxLxLx(b)求2点的DFT)()4/()()()()()()()4/10)()() 12()2()()

17、(62/5262/52242/3142/314/0) 12(2/114/022/1111LXWLXNkxLXWLXkxrxLXWLXNkXLLXWLXWLXWLXkXkXrxkNkNkNkNNLkLNNLLkNDFT）也可分解为：同样，（同理：，其中（可分解为：（c)一个2点的DFT蝶形流图2点DFT2点DFTx(0)x(4)x(2)x(6)X3(0)X3(1)X4(0)X4(1)X1(0)X1(1)X1(2)X1(3)04W14W) 1 () 1 () 3 () 0() 0() 2() 1 () 1 () 1 () 0() 0() 0(41431404314143140431XWXXXWXX

18、XWXXXWXX其中（d)另一个2点的DFT蝶形流图2点DFT2点DFTx(1)x(5)x(3)x(7)X5(0)X5(1)X6(0)X6(1)X2(0)X2(1)X2(2)X2(3)04W14W) 1 () 1 () 3 () 0() 0() 2() 1 () 1 () 1 () 0() 0() 0(61452604526145260452XWXXXWXXXWXXXWXX其中同理：(3)将N/4（2点）DFT再分解成2个1点的DFT(a)求2个一点的DFT021022120202120230202020231021 , 0; 1 , 0)4()0()4()0() 1 ()4()0()4()0

19、()0()()(2WWknWWWWWxWxWxWxXWxWxWxWxXWnxkXnknknkNnkNnnkN，其中，则这里用到对称性这是一蝶形结代入上面流图可知： 最后剩下两点DFT,它可分解成两个一点DFT，但一点DFT就等于输入信号本身，所以两点DFT可用一个蝶形结表示。取x(0)、x(4)为例。(b)2个1点的DFT蝶形流图 1点DFTx(0)1点DFTx(4)X3(0)X3(1)02W进一步简化为蝶形流图：02WX3(0)X3(1)x(0)x(4)4()0()4()0() 1 ()4()0()4()0()0(023023xxxWxXxxxWxX其中：(4)一个完整N=8的按时间抽取FF

20、T的运算流图 x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)38W28W18W08W08W08W08W08W08W08W28W28WX(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)12/0NNNWW其中旋转因子，共有m=0m=1m=2n 4.2.4 DITFFT的运算规律及编程思想n 1.原位计算n 由图4.2.4可以看出，DITFFT的运算过程很有规律。N=2M点的FFT共进行M级运算，每级由N/2个蝶形运算组成。同一级中，每个蝶形的两个输入数据只对本蝶形有用，且每个蝶形的输入、输出数据结点又在同一水平线上，计算完一个蝶形后，所得输出数据可立即存入原数据占用

21、的存储单元，这就是位计算。原位计算可大节省内存单元，降低成本。n 原位运算（in-place)n原位运算的结构，可以节省存储单元，降低设备成本。n定义：当数据输入到存储器以后，每一组运算的结果，仍然存放在这同一组存储器中直到最后输出。02Wx(0)x(4)X3(0)X3(1)中放在一个暂存器中放在单元中，将放在单元例：将RWAxAx08) 1 ()4()0()0(单元送回单元送回将) 1 ()4()0()0()4()0(0808AxWxAxWxn例：N=8 FFT运算,输入:x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)

22、A(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)=x(0)A(1)=x(1)A(2)=x(2)A(3)=x(3)A(4)=x(4)A(5)=x(5)A(6)=x(6)A(7)=x(7)382818084,3,2,1WRWRWRWR暂存器R1R1R1R1R1R3R1R1R3R2R3R4看出：用原位运算结构运算后，A(0)A(7)正好顺序存放X(0)X(7),可以直接顺序输出。2.旋转因子的变化规律 如上所述，N点DITFFT运算流图中，每级都有N/2个蝶形。每个蝶形都要乘以因子WpN，称其为旋转因子，

23、p称为旋转因子的指数。 观察图4.2.4不难发现，第L级共有2L-1个不同的旋转因子。N=23=8时的各级旋转因子表示如下： L=1时，WpN=WJ N/4=WJ2L， J=0 L=2时， WpN =WJ N/2=WJ2L， J=0，1 L=3时， WpN =WJN=WJ2L， J=0，1，2，3 382818082814080402080,;2, 1 ,0,2,; 1 ,0,2;0, 1,8:WWWWJLWWWWJLWWWJLNN如LMLJNJNpNMLMLMLLppNJPJWWWNJWWLMMLL212, 2 , 1 , 0,222212, 2 , 1 , 0,12212对N=2M的一般情

24、况，第L级的旋转因子为 3. 蝶形运算规律 设序列x(n)经时域抽选(倒序)后，存入数组X中。如果蝶形运算的两个输入数据相距B个点，应用原位计算，则蝶形运算可表示成如下形式： XL (J) XL-1(J)+X L-1(J+B)WpN XL(J+B) XL-1(J)-X L-1(J+B)WpN 式中 p=J2 M-L；J=0,1,，2 L-1-1；L=1,2,，M 下标L表示第L级运算，XL(J)则表示第L级运算后数组元素X(J)的值。如果要用实数运算完成上述蝶形运算，可按下面的算法进行。 设: T=X L-1(J+B)WpN=TR+jTI X L-1(J)=XR(J)+jXI(J) 式中下标R

25、表示取实部，I表示取虚部，pNBJXpNBJXTpNBJXpNBJXTpNjpNBJjXBJXWBJXTRIIIRRIRpNL2sin)(2cos)(2sin)(2cos)(2sin2)cos()()(1IIIRRRIIRRIRIRpNLLIRLIIIRRRIIRRIRIRpNLLIRLTJXBJXTJXBJXTJXjTJXjTTJjXJXWBJXJXBJjXBJXBJXTJXJXTJXJXTJXjTJXjTTJjXJXWBJXJXJjXJXJX)()()()()()()()()()()()()()()()()()()()()()()()()()(1111则WN0WN1WN2WN3WN0WN2

26、WN0WN2WN0WN0WN0WN0 x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)n 4. 编程思想及程序框图图4.2.6 DITFFT运算和程序框图 开 始送入x(n)， MN2 M倒 序L1 , M0 , B 1P2 M LJk J , N1 , 2LpNpNWBkXkXBkXWBkXkXkX)()()

27、()()()(输 出结 束B 2 L1第第L级中，每个蝶形的两个输级中，每个蝶形的两个输入数据相距入数据相距B=2L-1个点；同一个点；同一旋转因子对应着间隔为旋转因子对应着间隔为2L点点的的2M-L个蝶形。个蝶形。L=1,B=1L=2,B=2L=3,B=4n 5. 序列的倒序n DITFFT算法的输入序列的排序看起来似乎很乱，但仔细分析就会发现这种倒序是很有规律的。由于N=2M，所以顺序数可用M位二进制数(nM-1nM-2n1n0)表示。 图4.2.7 形成倒序的树状图(N=23) 01010101010101(n2n1n0)200004261537100010110001101011111

28、表4.2.1 顺序和倒序二进制数对照表 图4.2.8 倒序规律 x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)倒序倒序子程序具体执行时，只须将1与4对调送入，3与6对调送入，而0，2，5，7不变，仅需要一个中间存储单元。I01234567J04261537 在实际运算时，先按自然顺序将输入序列存入存储单元，再通过变址运算将自然顺序变换成按时间抽取的FFT算法要求的顺序。变址的过程可以用程

29、序安排加以实现-称为“整序”或“重排”（采用码位倒读）且注意：(1)当I=J时，数据不必调换；(2)当IJ时，必须将原来存放数据x(n)送入暂存器R，再将x(J)送入x(I),R中内容送x(J).进行数据对调。(3)为了避免再次考虑前面已调换过的数据，保证调换只进行一次，否则又变回原状。IN/2=4，做做4点的点的DIF：nNWNnxnxnxNnxnxnx)2()()()2()()(21先将N=8DFT分解成2个4点DFT：可知：时域上：x(0),x(1),x(2),x(3)为偶子序列 x(4),x(5),x(6),x(7)为奇子序列 频域上：X(0),X(2),X(4),X(6)由x1(n)

30、给出 X(1),X(3),X(5),X(7),由x2(n)给出将N=8点分解成2个4点的DIF的信号流图 4点DFTx(0)x(1)x(2)x(3)4点DFTx(4)x(5)x(6)x(7)X(0)X(2)X(4)X(6)X(1)X(3)X(5)X(7)X1(k)前半部分序列后半部分序列nNnNnNnNWxxxWxxxWxxxWxxxxxxxxxxxxxxx)7()3(3,)6()2(2,)5() 1 (1,)4()0(0)7()3(3),6()2(2),5() 1 (1),4()0(022221111）（）（）（）（）（）（）（）（如：08W18W28W38Wx1(n)x2(n)X2(k)(

31、2)N/2(4点)-N/4(2点)FFT(a)先将先将4点分解成点分解成2点的点的DIF：n因为4点DIF还是比较麻烦，所以再继续分解。后半部分序列、前半部分序列、) 3 () 2() 1 () 0(: )(11111xxxxnx10140)2()()()2()(411311，），在此（）（若设：LNLLxWNnxnxLxNnxnxnN后半部分序列、前半部分序列、同理：) 3 () 2() 1 () 0(: )(22222xxxxnx10140)2()()()2()(622522，），在此（）（同理：LNLLxWNnxnxLxNnxnxnN（b)一个2点的DIF蝶形流图2点DFT2点DFTx1

32、(0)x1(1)x1(2)x1(3)x3(0)x3(1)x4(0)x4(1)X(0)X(4)X(2)X(6)08W28W) 3 () 1 () 1 (,) 2() 0() 0(113113xxxxxx其中(c)另一个2点的DIF蝶形流图2点DFT2点DFTx2(0)x2(1)x2(2)x2(3)x5(0)x5(1)x6(0)x6(1)X(1)X(5)X(3)X(7)08W28W(3)将N/4（2点）DFT再分解成2个1点的DFT(a)求2个一点的DFT021022120202120230202020231021 , 0; 1 , 0)4()0()4()0() 1 ()4()0()4()0()0

33、()()(2WWknWWWWWxWxWxWxXWxWxWxWxXWnxkXnknknkNnkNnnkN，其中，则这里用到对称性这是一蝶形结代入上面流图可知： 最后剩下两点DFT,它可分解成两个一点DFT，但一点DFT就等于输入信号本身，所以两点DFT可用一个蝶形结表示。取x3(0)、x3(1)为例。(b)2个1点的DFT蝶形流图 1点DFTx3(0)1点DFTx3(1)X(0)X(4)02W进一步简化为蝶形流图：02WX(0)X(4)x3(0)x3(1)(4)一个完整N=8的按频率抽取FFT的运算流图 x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)38W28W18W08W08

34、W08W08W08W08W08W28W28WX(0)X(4)X(2)X(6)X(1)X(5)X(3)X(7)12/0NNNWW其中旋转因子，共有m=0m=1m=2图4.2.14 DITFFT的一种变形运算流图WN0WN0WN2WN0X(0)X(4)X(2)X(6)X(1)X(5)X(3)X(7)x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)WN0WN2WN1WN3WN2WN0WN0WN0按时间抽取的FFT算法的若干变体 对于任何流图，只要保持各节点所连续的支路及其传输系数不变，则不论节点位置怎么排列所得流图总是等效的，最后所得结果都是x(n)的离散付里叶变换的正确结果。只是数

35、据的提取和存放的次序不同而已。图4.2.15 DITFFT的一种变形运算流图WN0WN0WN2WN0X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)WN0WN2WN1WN3WN2WN0WN0WN0DIF的特点(a)输入自然顺序，输出乱序且满足码位倒置规输入自然顺序，输出乱序且满足码位倒置规则。则。(b)根据时域、频域互换，可知：根据时域、频域互换，可知：输入乱序，输出自然顺序。输入乱序，输出自然顺序。DIF与DIT比较n相同之处：n（1）DIF与DIT两种算法均为原位运算。n（2）DIF与DIT运算量相同。n它们都

36、需要算法是两种等价的与次复加次复乘FFTDITDIFNaNmNFNF22lglg2n不同之处：n（1）DIF与DIT两种算法结构倒过来。nDIF为输入顺序，输出乱序。运算完毕再运行“二进制倒读”程序。nDIT为输入乱序，输出顺序。先运行“二进制倒读”程序，再进行求DFT。n（2）DIF与DIT根本区别：在于蝶形结不同。nDIT的复数相乘出现在减法之前。nDIF的复数相乘出现在减法之后。n 4.2.6 IDFT的高效算法 n 上述FFT算法流图也可以用于离散傅里叶逆变换(Inverse Discrete Fourier Transform，简称IDFT)。比较DFT和IDFT的运算公式： IDF

37、TFFTNWWDFTWnxnxDFTkXWkXNkXIDFTnxnkNnkNNknkNNknkN算法都可以拿来运算或频率抽取抽取）那么以上讨论的时间（将运算结果都除以改成运算中的每个系数只要把3)2() 1 ()()()()(1)()(1010图4.2.16 DITIFFT运算流图 WN0WN1WN2WN3WN0WN0N1x(0)x(4)x(2)x(6)x(4)x(5)x(3)x(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)WN2WN2N1N1N1N1N1N1N1图4.2.17 DITIFFT运算流图(防止溢出) WN02121x(0)x(4)x(2)x(6)x(1)x

38、(5)x(3)x(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)212121WN121WN221WN3212121WN021WN2212121WN021WN22121WN02121WN0212121WN021WN021n 如果希望直接调用FFT子程序计算IFFT，则可用下面的方法：n 由于 10101( )( )1( )( )NknNkNknNkx nX k WNx nXk WN对上式两边同时取共轭，得1011( )( )( )NknNkx nXk WDFT XkNN可知：只须将频域成份一个求共轭变换，即（1）将X(k)的虚部乘以-1，即先取X(k)的共轭，得X*(k)。

39、(2)将X*(k)直接送入FFT程序即可得出Nx*(n)。(3)最后再对运算结果取一次共轭变换，并乘以常数1/N，即可以求出IFFT变换的x(n)的值。直接利用FFT流图方法的注意点n（1）FFT与IFFT连接应用时，注意输入输出序列的排列顺序，即应注意是自然顺序还是倒序。n（2）FFT和IFFT共用同一个程序时，也应注意利用FFT算法输入输出的排列顺序，即应注意自然顺序还是例位序n（3）当把频率抽取FFT流图用于IDFT时，应改称时间抽取IFFT流图。n（4）当把时间抽取FFT流图用于IDFT时，应改称频率抽取IFFT流图。4.3 进一步减少运算量的措施n 4.3.1 多类蝶形单元运算n 由

40、DITFFT运算流图已得出结论，N=2M点FFT共需要MN/2次复数乘法。n 由(4.2.12)式，当L=1时，只有一种旋转因子W0N=1，所以，第一级不需要乘法运算。 n 综上所述，先除去第一、二两级后，所需复数乘法次数应是n 从L=3至L=M共减少复数乘法次数为(2)(2)2MNCM13312( )2222MMLLLLNNN因此，DITFFT的复乘次数降至 (2)(2)(2)(3)2222MNNNCMM22(1)()()222()()22()222()()22defjxjyxjyjxyxyj xyRjIRxyIxyyx n 从实数运算考虑，计算N=2M点DITFFT所需实数乘法次数为(2)

41、4(3)2(2)2213(2)102MNNRMNMn 4.3.2 旋转因子的生成n 在FFT运算中，旋转因子WmN=cos(2m/N)-jsin(2m/N)，求正弦和余弦函数值的计算量是很大的。 n 4.3.3 实序列的FFT算法 n 设x(n)为N点实序列，取x(n)的偶数点和奇数点分别作为新构造序列y(n)的实部和虚部，即1212( )(2 ),( )(21),0,1,12( )( )( ),0,1,12Nx nxnx nxnnNy nx njx n n对y(n)进行N/2点FFT，输出Y(k)，则1122( )( )( ),0,1,1( )( )( )2epopX kDFT x nYkN

42、kXkDFT x njYk 根据DITFFT的思想及式(4.2.7)和(4.2.8)，可得到 12( )( )( ),0,1,2kNNX kX kW Xk kn实信号的快速循环卷积实信号的快速循环卷积n 用DITFFT计算两个实信号x1(n)和x2(n)的循环卷积时，最直接的方法是将信号的虚部都置为零，再按复序列用FFT计算。显然这样做是很浪费的。现在我们可用基2DITFHT直接进行实正交变换来处理实信号的循环卷积问题。1122( )( )( )( )FHTHFHTHx nXkx nXk计算式 )()()()(21nxnxnykYIFHTH n%conv2nfunction y=circonv

43、2(x1,x2,N)nif length(x1)Nn error(N must not be less than length of x1)nendnif length(x2)Nn error(N must not be less than length of x2)nendX1k=fft(x1,N);X2k=fft(x2,N);Yk=X1k.*X2k;y=ifft(Yk);if(all(imag(x1)=0)&(all(imag(x2)=0) y=real(y);endxn=1,2,1,0;n=size(xn);N=n(2); %取得X(k)的长度subplot(2,1,1);ste

44、m(xn)Xk=fft(xn) %计算FFTX(k)magXk=abs(Xk)subplot(2,1,2)stem(real(Xk); x1=1,1,1,1,0,0,0,0;n=0:16;x2=cos(pi*n/4)+cos(pi*n/8);i=0:7;FFT的应用n凡是利用付里叶变换来进行分析、综合、变换的地方，都可以利用FFT算法来减少其计算量。nFFT主要应用在：n（1）快速卷积n（2）快速相关n（3）频谱分析一、快速卷积n设x(n)的长度为N点，h(n)的长度为M点，则：ny(n)的长度为N+M-1点。所以直接计算y(n)线性卷积运算量为NM。1.利用圆周卷积代替线卷积n用圆周卷积（F

45、FT）代替线卷积的必要条件：x(n)、h(n)补零加长至L=N+M-1.n然后计算圆卷积。求出y(n)代表线卷积。1010)()(1010)()(LnNNnnhnhLnNNnnxnx2、用FFT计算y(n)的步骤n(1)求H(k)=DFTh(n) (L点） n(2)求X(k)=DFTx(n) (L点）n(3)计算Y(k)=X(k)H(k) (L点）n(4)求y(n)=IDFTY(k) (L点）2、用FFT计算y(n)与直接计算y(n)的运算量比较称分段过滤的方法。这时须采用分段卷积或的优点。则体现不出其圆周卷积长度差很多时，与）当当（的优点。可以体现出其圆周卷积长度差不多时，与）当（计算运算量

46、用直接计算运算量求线卷积运算量用)()(2)()(1log231)1(log23)1(22nhnxnhnxMNNMFFTKLLFFTMNmL3、分段卷积的方法n（1）重叠相加法n（2）重叠保留法二、快速相关1.方法1010)()(1010)()(21)()()(:)()(10*LnMMnnynyLnNNnnxnxmLMNLFFTmnymxnrMnyNnxmLmyx为整数），令且，来代替线性相关，选择圆周相关法来求线性相关，是用利用相关点，要求线性的长度为点，的长度为设2.步骤自相关序列。则求得若（而得到。后，取共轭再乘的可以利用即：，来求程序计算同样可以用已有的点求求乘积点求并求点求),()(

47、),()()1)()(1)(1)()()()()()()()()()()(),()()(*)(*10*)()(nrnrnynxeNFFTkRrWkRNWkRNrrIFFTFFTkRIDFTnrIFFTNdkYkXkRcnyDFTkYFFTNbkXnxDFTkXFFTNaxxxyyxnyxNknkNyxnkNxynyxnyxyxyxyx三、频谱分析n这里我们仅介绍FFT的仪器设备：快速付里叶分析仪。n其功能为：n（1）能分析确定性信号的频谱。n（2）估计随机信号的功率谱n（3）并对信号进行快速卷积滤波n（4）计算相关函数1.频谱分析仪的框图数据选择器A/D保护滤波器输入电路输入数据存储器运算器数

48、据选择器控制器变址单元函数发生器输出缓冲器D/A输出2.部件说明n（1）保护滤波器：对输入信号进行模拟滤波，滤掉噪声，提取感兴趣的有用信号，以便分析频谱。n（2）运算器：完成时间抽取FFT或Chirp-Z变换所要求运算。n（3）RAM:存储输入数据。n（4）ROM(函数发生器）。以表格形式存放蝶形运算因子W.n(5)变址单元：根据输入数据，进行码位倒置排序。.2)()()2()()1.(2)()(2)(.41IDFTNnxIFFTNkXkXFFTNDFTNnxkXNnx点的实现求点设计用一次，若已知高效算法。完成计算点的设计用一次点的为的有限长序列，是长度为设：习题例)()(21)()()()

49、()(21)()()(1, 2 , 1 , 0),()()()()().()()()(1, 2 , 1 , 0),12()(1, 2 , 1 , 0),2()()()()(122112121212121kNYkYkYnjxDFTkjXkNYkYkYnxDFTkXNknyDFTkYnjxnxnykXkXFFTNDFTnxnxNnnxnxNnnxnxnxnxNnxFFTDITopep令：和求得点可用一次的共轭对称性，均为实序列，根据和和点实序列得到两个点和奇数点）在时域分别抽取偶数（解：解题思路是),()()(1, 1 , 0),()()(221221kXWkXNkXNkkXWkXkXkNkN1,

50、 2 , 1 , 0),()(1, 2 , 1 , 0),()(1, 2 , 1 , 0),12()(1, 2 , 1 , 0),2()(2221121NknxDFTkXNknxDFTkXNnnxnxNnnxnx）（kNkNkNWNkXkXkXNkXkXkXkXWkXNkXNkkXWkXkX221221221)()(21)()()(21)(),()()(1, 1 , 0),()()(1, 1 , 0)(Im)(Re)()()()()()()()()()()()()(21)()()(21)()()(2121221NnnyjnykYIFFTnykYkYkjXkXkYkYNkXkXbWNkXkXkX

51、NkXkXkXkXaopepkN点频域序列构成和由计算由运算过程如下：120)21()2()()()()()()()()()(21)(Im)()()()(21)(Re212121NnnnxnnxnxnxnxnxcnjxkYDFTnynynyjnxkYDFTnynynyopep，奇数，偶数，合成和由组元素中即可。的数组的偶数和奇数数依次放入的两个数组元素分别和编程时，只要将存放)()()(21nxnxnx习题习题4.10%Xk:被变换数据被变换数据X(k)%XN:IFFTX(k)结果结果x(n)%N:x(n)/X(k)长度长度Xk=X(0) X(1) X(N-1);n=size(Xk);N=n(

52、2); %取得取得X(k)的长度的长度Xk=conj(Xk); %取复共轭取复共轭XN=fft(Xk); %计算计算FFTX(k)XN=conj(XN)/N;stem(real(XN); %绘制绘制x(n)序列波形图序列波形图数据向量数据向量Xk也可通过键盘、数据文件或函数计算等多种方法也可通过键盘、数据文件或函数计算等多种方法建立，本程序中直接赋值建立，本程序中直接赋值Xk向量。向量。n 如果希望直接调用FFT子程序计算IFFT，则可用下面的方法：n 由于 10101( )( )1( )( )NknNkNknNkx nX k WNx nXk WN对上式两边同时取共轭，得1011( )( )( )NknNkx nXk WDFT XkNN可知：只须将频域成份一个求共轭变换，即（1）将X(k)