05-第五章快速傅里叶变换(蝶形运算)

1、第五章快速傅里叶变换本章目录直接计算DFT的问题及改进的途径 按时间抽取的基2-FFT算法 按频率挾取的基2-FFT算法 快速傅里叶逆变换(IFFT)算法 Matlab 实现5J引言 DFT在实际应用中很重要:可以计算信号的频谱、功率谱和线性卷积等。直接按DFT变换进行计算，当序列长度汕艮大时，计算量非常大，所需时间会很长。 FFT并不是一种与DFT不同的变换，而是DFT的一种快速计算的算法。竺亘斟算DFT的问题及改进的途径DFT的运算量5#设复序列MG)长度为N点，其DFT为N1X(Q = 2>(曲n=0k=0,，AM(1)计算一个X(k)值的运算量复数乘法次数：N复数加法次数：NT5

2、.2J DFT的运算量(2)计算全部N个X(k)值的运算量复数乘法次数：/V2复数加法次数：N(N1)(3)对应的实数运算量NlN7X 伙)=工兀()比篇=Re x(n) + jlmx(n)RcWk + jImWkn=0n=0N-l= Rex(n)ReWk Imx(n)ImWkw=0+yRe x(n) Im Wk + Im x(n) Re7一次复数乘法：4次实数乘法+ 2次实数加法一个X(k):4N次实数乘法+2A/+2(AM)= 2(2 AM)次实数加法所以整个N点DFT运算共需要：实数乘法次数：4/V2实数加法次数：NX2(2AM)=2N(2AM)9*2!算量的结论N点DFT的复数乘法次数

3、举例NN22441686416256321028NN26440491281638425665 536512262 14410241 048 576结论：当M艮大时，其运算量很大，对实时性很强的信号 处理来说，要求计算速度快，因此需要改进DFT的计算 方法，以大大减少运算次数。115.2.2减少运算工作量的途径主要原理是利用系数W护的以下特性对DFT进行分解:(1)对称性(wky = wnk =(2)周期性(3)可约性VVNwW+N)一 VVN另外,w;k 呼=w-nktmN / tn</2=-l 呼+ N/2)=肌5.3按时间抽取的基ZFFT算法算法原理 按时间抽取基-2FFT算法与直接

4、计算D FT运算量的比较按时间抽取的FFT算法的特点按时间抽取FFT算法的其它形式流程图5.3J算法原理设N=2J将x（门）按n的奇偶分为两组:x(2r)=召(厂) x(2r + 1) = x2(r)N-iX伙)=DFTxM=工兀昭*7?=0NN-=工兀（）*+工兀）”解/?=0n=0几为偶数料为奇数15n=0n为奇数2=2何呼+ 2（）呼A?=0n为偶数2=S 双2门叭"+ 工 X2r + l）W；2r+1U r=0r=07_1=2心）畔+观工勺（厂）",=X 伙）+呎X2伙）4伫1r=02r=022式中，E（Q和蜀仇）分别是re）和工2何的N/2的DFT。另外，式中氐的

5、取值范围是：0, 1,，N/2-1 o此，X(Q = X") +说/伙)只能计算出x(k)的前一半值。后一半X(k)1¥N/2 , N/2 +1, N ?利用 可得到J 2亠JyV/2-1N/21X(- + k)= £ 兀1(门叭；尸)=£ x“)W篇2 = X”)2r=0r=0同理可得X2(y+ £) = X2(k)考虑到W严k) = W閉2 ,Wk =_Wk及前半部分X X(k) = X(k) + WX2(k)k=0, 1,N/2-1因此可得后半部分X(QX(k + -) = Xl(k+) + W+N/2X2 (k + y)= Xl(k)-

6、WX2(k)k=Qf 1,N/2-119以8点为例第_次按奇偶分解蝶形运算X伙)=/伙)+叱X?伙)蝶形运算式X(k) = x(切-吩2 伙)蝶形运算信 号流图符号#以8点为例第_次按奇偶分解#以8点为例第_次按奇偶分解因此，只要求出2个N/2点的DFT,即E(k)和X2(k),再 经过蝶形运算就可求出全部X(k)的值，运算量大大减少。#以8点为例第_次按奇偶分解21以8点为例第_次按奇偶分解兀(0) =x(0) >xi(l)=x(2)>Xj(2) =x(4)>X|(3) =x(6) <x2(0) =x(1) x2(1)=x(3)<X2(2) =x(5).以3)

7、=x(7) 葺点罟点b(0)M)川2)DFT加3)总3)2(0)乂X(5)DFT06)X(7)X|(l)基X(0)i/2)以N=8为例,分解为2个4点 的DFT,然后 做8/2=4次蝶形运算即可求出 所有8点X(k)的 值。#以8点为例第_次按奇偶分解蝶形运算量比较触DFT的运算量复数乘法次数：N2复数加法次数：N(N1)分解一次后所需的运算量=2个N/2的DFT+23以8点为例第_次按奇偶分解N/2蝶形:复数乘法次数：2*(M2尸+N/2二N2/2+M2复数加法次数：1)+2*"/2=2/2因此通过一次分解后，运算工作量减少了差不多一半。#进_步按奇偶分解由于N=2S因而N/2仍是

8、偶数，可以进一步把每个N/2点 子序列再按其奇偶部分分解为两个N/4点的子序列。2 = 0,1,，理14以N/2点序列可(厂)为例再)=召x1(2Z + 1) = x4(Z)则有N/2-N/4-1N/4-1X(k)=血=此扇 + + 1)肌$讥r=0/=0/=0N/4-1N/41=陀4+W怎陀41=01=0= X3(k) + W/2X4(k) B0丄,J-125以8点为例第二次按奇偶分解(N A,N 5X - + k =X3(k)-W爲X4伙)f ,-1I 4 丿由此可见，一个N/2点DFT可分解成两个N/4点DFT。同理，也可对兀2(n)进行同样的分解，求出X2(k)027以8点为例第二次按

9、奇偶分解#以8点为例第二次按奇偶分解应(0)=兀1(1) =x(2) x4(l) =X|(3) =x(6) %3(0)=%|(0)=兀(0) 如1)x3(l)=xi(2)=x(4)e无3（°）益(0)W)孚点DFTML2WnX1Xi-1冷3)抡(0)E占4小、DFT疋X(0)x5(0) =x2(0) =x(l)x5(1)=x2(2) =x(5)/占4小、兀6(0) =X2(1)=X(3) 兀6(1)=七(3) =x(7) DFT屁花(0)总1)无(2)X(3)X(4)X(5)X(6)-1WN 赦1)WN XA2(3)x（0）X29以8点为例第二次按奇偶分解算法原理#以8点为例第二次按

10、奇偶分解对此例N=8,最后剩下的是4个N/4=2点的DFT, 2点 DFT也可以由蝶形运算来完成。以禺仇)为例。N/41N/411心伙)二工兀3W/4二工律/4上=°，1/=01=0即X3 (0) = x3 (0) + 必® (!)=班°)+ 必双4) = x(0) + W：x(4)X3=%3(°)+昭兀3=x(0) +岡兀(4) =%(0)-叱汝这说明，N=2M的DFT可全部由蝶形运算来完成。31N=8按时间抽取法FFT信号流图5.3.2按时间抽販基ZFFT算法与H接计算DFT运算豪的比较2次复加，因此每级运算都需M2次复乘和由按时间抽取法FFT的信号

11、流图可知，当N=2厶时，共有L级 蝶形运算；每级都由 M2个蝶形运算组成,而每个蝶形有 1次复乘、N次复加。33这样丄级运算总共需要：复数乘法：-L = -log22 2复数加法：N直接DFT算法运算量复数乘法：N1复数加法：N(N1)直接计算DFT与FFT算法的计算量之比为M“ N?2N10g2”=g2N+ FFT算法与直接DFTX法运算量的比较-NN2.N、,计算量 之比MNN2N、. 三呃N计算量 之比M2414.012816 38444836.641644.025665 5361 02464.0864125.4512262 1442 304113.816256328.010241 04

12、8 5765 120204.83210288012.820484 194 30411 264372.464404919221.4355.3.3按时间抽取的FFT算法的特点序列的逆序排列同址运算（原位运算）蝶形运算两节点间的距离 C的确定#I序列的逆序排列序列的逆序排列由于欢）被反复地按奇、偶分组，所以流图输入端的排列不再是顺序的，但仍有规律可循:37#因为N=2M ,对于任意n(0<n<A-l),可以用M个#二进制码表示为:n(DEC) = (“M-inM-2 A 2“l“0)(B/N) fo1反复按奇、偶分解时，即按二进制码的分解。#倒位序的树状图(N-8)#01x(切屮

13、6;)1"2-一x(OOO)-一 x(100)-一一 x(010)-x(110) 0-一一 x(001)班101)-一兀(011) x(lll)#码位的倒位序(Nm8)自然顺序n二鯉数倒位序二进制数 倒位序顺序数分000000001100010010011110100001101101110Oil11111139倒位序的变址处理(N8)41#存储单元自然顺序单元力兀(0)变址Tx(0)倒位序/(8)x(7)Tx(7)A(2) A(3) A(4) A(5) A(6) A(7) x(l) x(2) x(3) x(4) x(5) x(6)x(4)兀(2) x(6) x(l) x(5) x(

14、3)同址运算（原位运算）同址运算（原位运算）某一列任何两个节点氐和/的节点变量进行蝶形运算 后，得到结果为下一列肛/两节点的节点变量，而和其他 节点变量无关。这种原位运算结构可以节省存储单元， 降低设备成本。43#运算后X伙）一4伙）#观察原位运算规律双0)X4)-1双2)兀XI)双5)x(3)兀wl1-1总0)M)X(3)X(4)十 X(6)X(7)45#蝶形运算两节点间的距离蝶形运算两节点间的距离以N=8为例:第一级蝶形,距离为:第二级蝶形,距离为:第三级蝶形,距离为:规律：对于共L级的蝶形而言，其m级蝶形运算的节点间的距离为2心忙的确定 C的确定以N=8为例：加=10寸，Wr=Wj/4=

15、WJm =<,j = 0加=20寸，W=W/2=Wm =Wj = O,l加=3时，W=W =Wm =Wj = O丄2,3N = 2mL 级：肌二昭，丿=0丄2,A ,2i-l&2L =2M x 2lm =Nx 2lm5.4按频率抽取的基2-FFT算法算法原理先把输入按门的顺序分成前后两半再把输出X(k)按幺的奇偶分组设序列长度为N=2J L为整数47#f前半子序列只门)NI后半子序列血+亍N 1ox -10<n< 乡 T2495.4J算法原理由DFT定义得7V-1X伙） = 0（曲n=0N/2-1N-=£ x（）W/ + 工 x（）Wf/?=0n=N 12皆

16、 1, Ng N 5+3=£班砒叭+ 2兀（+亏）必2 n=0mn厶N/21/?=0n=0- N"TN kE心）+兀（十3）%2 叱yk=0 9 1,，N由于N/2-x(k)= £/?=0x(n) + x(n + )WjkN 2 =严=_51#所以N 12-X閃=Sn=0x(n) + (-1)A x(n + )Kk#k=0，1,，N#然后按氐的奇偶可将X(Q分为两部分则式k = 2rN .r=0, 19，一 1k = 2厂+ 12nkN/2-1Nx(k)= 2 x(n) + (-l)kx(n + -) W；53/?=0可转化为#N oX(2"二工 x(n

17、)x(n + )吒N/2-171=0nrN/21X"=0Nx(n) + xn + ) W.-nrN/2#n(2r+)NN/2-1NX(2r + 1)=工x(n)-x(n + )-h=o2N/2-1NX 兀-兀+ )%' 昭;2/?=o2#N(n) = x(n) + x(n + )N x2(n) = x(n) -xn + )代入N/2-1X(2r)=工打=0NgNX(2r + 1)=工Wn)-x(H + -W；.C；2x(h) + x(n + )必;2n=0可得心X(2r) = 2>(叫277=0 N/2-1X(2r + 1)=工也睥77=0r=0, 1为2个N/2点的D

18、FT,合起来正好是N点X(Q的值。55蝶形运算N£ (/7)= x(n) + x(n + )称为蝶形运算57#与时间抽选基2FFT算法中的蝶形运算符号略有不同。右例按频率抽取(N8)例 按频率抽取，将N点DFT分解为两个N/2点DFT的组合(N=8)59#N占-2小、DFTI/O)i/2)i/4)i/6)多点DFTi(l)f(3)i(5)i/7)#与时间抽取法的推导过程一样，由于N=2S N/2仍然是一个偶数, 与奇数组,因而可以将每个N/2点DFT的输出再分解为偶数组这就将N/2点DFT进一步分解为两个N/4点DFT。抽取法与时间抽取法的异同频率抽取法输入是自然顺序，输出是倒位序

19、的；时间抽取法正好相反。频率抽取法的基本蝶形与时间抽取法的基本蝶形有所不同。频率抽取法运算量与时间抽取法相同。频率抽取法与时间抽取法的基本蝶形是互为转置的。右55快速傅里叶逆变换（IFFT）算法IDFT公式1 N1兀（町=IDFT X（£）=帚工 X （k）WtjkN k=oDFT公式NlX伙）=DFT 兀（训二工心）呼71=0比较可以看出,IDFT多出 N 二F =>-=-M个1/2可分解到M级蝶形运算中。例频率抽取IFFT流图(28)快速傅里叶逆变换另一种算法NIDFTX（k） = XW1 "T/=-工对伙）（"广）N L-Ck=0厂NlN1nkNk=0

20、k=0IFFTX 伙）=X伙）_杏丿割一x *伙）FFTX*伙）求弊 FFFX*伙）除以何今 x（n）5.8 Matlab实现用FFT进行谱分析的Matlab实现用CZT进行谱分析的Matlab实现在Matlab中使用的线性调频z变换函数为czt, 其调用格式为 »X= czt(x3 M, W, A)其中，x是待变换的时域信号x(n),其长度为N, M是变换的长度，W确定变换的步长，A确定变换的起点。M=N, A=1,贝!|CZT变成DFT。655.8J用FFT进行谱分析的Matlab实现例51设模拟信号X0 = 2sin(4奴)+ 5cos(8奴)，以t= O.Oln («

21、;=0: -1)进行取样，试用fft函数对其做频谱分析。N分别 为:(1)N=45； (2)N=50； (3)N=55； (2)N=60。程序清单如下1=#%计算N=45的F FT并绘出其幅频曲线N=45;n=0:N-1 ;t=O.O1*n;q=n*2*pi/N; x=2*sin(4*pi*t)+5*cos(8*pi*t);y=fft(x,N); figure(1)subplot(2,2,1) plot(q5abs(y) title(TFT N=45J例51程序清单%计算N=50的FFT并绘出其幅频曲线N=50;n=0:N-1 ;t=0.01*n; q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t);y=fft(x5N); figure subplot(2,2