《离散时间信号处理:第5章 离散傅立叶变换与快速算法2(FFT)》由会员分享,可在线阅读,更多相关《离散时间信号处理:第5章 离散傅立叶变换与快速算法2(FFT)(58页珍藏版)》请在金锄头文库上搜索。
1、15.2DFT的快速计算的快速计算(FFT)n一、一、DFT存在的问题及改进的途径存在的问题及改进的途径n二、戈泽尔算法二、戈泽尔算法n三、时域抽取基三、时域抽取基-2算法(库利算法(库利-图基算法)图基算法)n四、频域抽取基四、频域抽取基-2算法(桑德算法(桑德-图基算法)图基算法)n五、线性调频五、线性调频Z变换变换2nN点有限长序列,其点有限长序列,其DFT变换对为:变换对为:n n 可以看出,可以看出,变换与反变换变换与反变换的差别仅在于的差别仅在于 的指数符号和常数乘因子的指数符号和常数乘因子1/N。实际上:实际上:n 因而我们因而我们只讨论只讨论DFT正变换的运算量正变换的运算量,
2、反,反变换的运算量是完全相同的。变换的运算量是完全相同的。 5.2.1DFT存在的问题及改进途径存在的问题及改进途径3n一般来说一般来说 和和 都是复数,都是复数,n每一个每一个 值的计算,需要值的计算,需要N次复数乘次复数乘法法以及(以及(N-1)次复数加法,次复数加法,n完成整个完成整个DFT运算总需要运算总需要 次复数乘次复数乘法法和和N(N-1)次复数加法次复数加法。DFT直接计算的运算量4n复数运算实际上是由实数运算来完成的,如:复数运算实际上是由实数运算来完成的,如:n一次复数乘法需用:一次复数乘法需用:n四次实数乘法四次实数乘法n二次实数加法;二次实数加法;n一次复数加法则需:一
3、次复数加法则需:n二次实数加法。二次实数加法。n整个整个DFT运算总共需要:运算总共需要:n4N2次实数乘法次实数乘法n 次次实数加法实数加法。5n上述统计与实际的运算次数有上述统计与实际的运算次数有少许出入少许出入,某些,某些因子不能按照一般的复数来计算运算量因子不能按照一般的复数来计算运算量,如如:n当当N很大时,这种特例很小。很大时,这种特例很小。 n直接计算直接计算DFT运算量是很可观的运算量是很可观的:nN=8时,时,DFT需需64次复乘,次复乘,nN=1024时,时,DFT所需复乘为所需复乘为1 048 576次,次,n对对实时性很强的信号处理实时性很强的信号处理来说,对来说,对硬
4、件的计算硬件的计算速度要求是太高速度要求是太高了。为了实用的需要,了。为了实用的需要,需要改需要改进进DFT的计算方法,减少运算次数。的计算方法,减少运算次数。6如如何何才才能能减减少少运运算算量量?仔仔细细观观察察DFTDFT的的运运算算就就可可看看出,其出,其变换系数变换系数 的具有如下的具有如下特性特性:(1 1) 的的对称对称性性(2 2) 的的周期周期性性 = = =(3 3) 的的可约可约性性 由此可得由此可得: :变换基底的特性7n(1)利用)利用 的的对称性对称性,合并合并DFT运运算中某些项,并且可以计算反变换;算中某些项,并且可以计算反变换;n(2)由于)由于DFT的运算量
5、是与的运算量是与 成正成正比,利用比,利用 周期性周期性和和可约性可约性将长序将长序的的DFT分解为短序列分解为短序列的的DFT。nFFT基本可以分成基本可以分成两大类两大类n按时间抽取按时间抽取(Decimation-in-time,DIT)n按频率抽取按频率抽取(Decimation-in-frequency,DIF)。)。FFT的基本的基本出发点出发点85.2.1 戈泽尔算法戈泽尔算法戈泽尔算法是一个利用复变换基W的周期性减少运算量的典型例子可以用于单离散频点或少数任意离散频点场合(如DTMF辨识)的有效计算9变形10推论:n如果构造冲击响应h(n)为 的系统,则该系统对有限长输入在 时
6、刻的响应即为 。n 下图给出了该系统的计算流图:11n由上图可以看出计算一个X(k)需要N次复乘和N次复加,即4N次实乘和4N次实加。n比直接计算DFT的运算量稍大一些,但是却避免了计算或存储复系数 。n另外,立足于上面用卷积实现卷积实现DFT的方法,我们是有可能通过某种变形将上面的运算量减小一半。 12n注意到系统函数:13利用迭代实现了复因子的计算如果输入序列是复数,由于系数是实数并且-1不必作乘法运算,所以计算一个X(K)只需要2N次实数乘法和4N次实数加法。戈泽尔算法还有另外一个优点是只需要将前馈环节中的复因子取共轭就可以计算X(N-K),可以将运算量再次减半。戈泽尔算法虽然比直接运算
7、有效的多,但是无论如何变形,其运算量将仍然正比于N2。戈泽尔算法特点14n一、算法原理一、算法原理n 假设序列点数为假设序列点数为 ,L为整数。如果为整数。如果不满足,可以人为加上若干零值点,使不满足,可以人为加上若干零值点,使之达到这一要求。这种之达到这一要求。这种N为为2的整数幂的的整数幂的FFT也称基也称基-2 FFT。n 将将 的序列按的序列按n的奇偶的奇偶分成两组:分成两组: n n r=0,1,N/2-15.2.2 DIT的的2MFFT(库利库利-图基算法图基算法)1516n利用利用 的的可约性可约性:n式中式中 与与 均为均为N/2点点DFT:一分为二一分为二17n但但 以及以及
8、 , 都是都是N/2点的点的序列,即序列,即n而而 却有却有N点,因而点,因而计算全部计算全部的的 ,必须应用必须应用DFT隐含的周期性扩展:隐含的周期性扩展:定义的扩展18n同时考虑到同时考虑到 的以下的以下性质性质n n这样就可将这样就可将 表达为前后两部分:表达为前后两部分:n 1920 上面的运算可以抽象为下面的蝶形信号流图单元:上面的运算可以抽象为下面的蝶形信号流图单元: 可以看出,每个蝶形运算需要可以看出,每个蝶形运算需要一次复数乘法一次复数乘法 及及两次复数加两次复数加(减)法。(减)法。21nN/2点点DFT复乘复乘:n 复加复加:n蝶形运算蝶形运算 复乘复乘:n 复加复加:n
9、总运算量总运算量 复乘复乘:n 复加复加:n直接直接N点点DFT复乘复乘: N2 复加复加:N(N-1)运算量-减少一半22n既然如此,由于既然如此,由于 ,因此,因此N/2仍是仍是偶数,可以进一步偶数,可以进一步按时间按时间的的奇偶性分解奇偶性分解,并且并且一直进行下去一直进行下去。n由于这种方法每一步都是按输入序列由于这种方法每一步都是按输入序列时时间的奇、偶性间的奇、偶性分解为两个更短的子序列,分解为两个更短的子序列,所以称为所以称为“按时间抽选法按时间抽选法”(DIT)。)。n分解过程如下分解过程如下:231、原位运算,2、倒位序,3、叠型间距24n由按时间抽选法由按时间抽选法FFT的
10、流图可见,当的流图可见,当 时,共有:时,共有:nL级蝶形,每级都由级蝶形,每级都由N/2个蝶形个蝶形运算组成运算组成n每个蝶形有一次复乘、二次复加,因而每级每个蝶形有一次复乘、二次复加,因而每级运算都需运算都需N/2次复乘和次复乘和N次复加次复加nL级级FFT运算总共需要:运算总共需要:n 复乘数:复乘数: (DFT: )n 复加数:复加数: (DFT: )二、运算量二、运算量25n由于计算机上乘法运算所需时间比加法由于计算机上乘法运算所需时间比加法运算所需时间多得多,所以乘法是主要运算所需时间多得多,所以乘法是主要的运算量。的运算量。n直接计算直接计算DFT与与FFT算法计算量之比为算法计
11、算量之比为:2627三、三、DIT算法的特点算法的特点n1、原位运算(同址运算)原位运算(同址运算)n每级(每列)都由每级(每列)都由N/2个蝶形运算构成个蝶形运算构成n每一个蝶形结构完成下述基本迭代运算:每一个蝶形结构完成下述基本迭代运算: nm表示第表示第m级迭代,级迭代, k、j为数据所在行为数据所在行28n按原位计算时,按原位计算时,FFT的的输出输出 是按是按正常顺正常顺序排列序排列在存储单元中在存储单元中:n X(0),),X(1),),X(7),n但是但是输入输入 却不是按自然顺序存储的却不是按自然顺序存储的:n看起来好像是看起来好像是“混乱无序混乱无序”,实际是有规律,实际是有
12、规律的,称之为的,称之为倒位序倒位序。n造成倒位序的原因是输入造成倒位序的原因是输入 按标号按标号n的的奇、奇、偶性的不断分组偶性的不断分组造成。造成。n倒位序的实现倒位序的实现2、倒位序规律、倒位序规律293、蝶形运算两结点的、蝶形运算两结点的“距离距离”n由由8点点FFT的信号流图可以看出,其输入的信号流图可以看出,其输入是倒位序的,输出是自然顺序。是倒位序的,输出是自然顺序。n第一级每个蝶形的两节点间第一级每个蝶形的两节点间“距离距离”为为1,n第二级每个蝶形的两节点第二级每个蝶形的两节点“距离距离”为为2,n第三级每个蝶形的两节点第三级每个蝶形的两节点“距离距离”为为4,n由此类推得,
13、对由此类推得,对N2L ,当输入为倒位当输入为倒位序,输出为正常顺序时,其序,输出为正常顺序时,其第第m级运算,级运算,每个蝶形的两节点每个蝶形的两节点“距离距离”为为 。30四、四、DIT算法的其他形式流图算法的其他形式流图n对于对于任何流图,只要保持各节点所连任何流图,只要保持各节点所连的各支路及其传输系数不变的各支路及其传输系数不变,则不论,则不论节点位置怎么排列所得流图总是等效节点位置怎么排列所得流图总是等效n所得最后结果都是所得最后结果都是DFT的正确结果,只的正确结果,只是是数据的提取和存放的次序不同而已数据的提取和存放的次序不同而已。3132335.2.3 DIF的的2M FFT
14、(桑德桑德-图基算法图基算法)n频率抽选(频率抽选(DIF)的的FFT算法,是把频域算法,是把频域序列(也是序列(也是N点)按点)按K的奇偶性分解为越的奇偶性分解为越来越短的序列。来越短的序列。 n一、算法原理一、算法原理n假设序列点数为假设序列点数为 ,L为整数。为整数。n在将在将 按按k的奇偶分组之前,先把输入的奇偶分组之前,先把输入按按n的顺序分成前后两半的顺序分成前后两半3435n由于由于 ,故,故 ,可得:,可得:n n k=0,1,N-1n当当k为偶数时,为偶数时, ;n k为奇数时,为奇数时, 。n因此,按因此,按 k的奇偶可将的奇偶可将 分为两部分。分为两部分。 36n令令 n
15、 n则:则:37n令令 n则:则:38DIF运算关系的基本蝶形运算关系的基本蝶形39N=8点点DIF FFT结构结构40二、二、DIT与与DIF的异同的异同n形式上的区别是:形式上的区别是:DIF输入是自然顺序,输入是自然顺序,输出是倒位序的,与输出是倒位序的,与DIT正好相反。正好相反。n但这不是实质性的区别,因为但这不是实质性的区别,因为DIF与与DIT都可将输入或输出按照要求进行重排。都可将输入或输出按照要求进行重排。n实质性的区别是:实质性的区别是:nDIF的基本蝶形与的基本蝶形与DIT的基本蝶形不同。的基本蝶形不同。414243作业:n9.6n9.7 n9.14n9.17n9.26n
16、9.27五、chirp Z变换n对对非单位圆上的抽样感兴非单位圆上的抽样感兴趣趣n语音信号处理中往往需要知道极语音信号处理中往往需要知道极点所在处的复频率,如果极点位点所在处的复频率,如果极点位置离单位圆较远,只利用单位圆置离单位圆较远,只利用单位圆上的频谱,就很难知道极点所在上的频谱,就很难知道极点所在处的复频率。处的复频率。nz变换采用变换采用螺线抽样螺线抽样就适就适应于这些需要,称为线性应于这些需要,称为线性调频调频z变换变换(CZT)4445n已知已知 是有限长序列,其是有限长序列,其z变换为:变换为:n为适应为适应z可以沿可以沿Z平面更一般的路径取值,沿平面更一般的路径取值,沿z平面
17、上的一段平面上的一段螺线螺线作作等分角等分角的的抽样抽样,z的这的这些抽样点些抽样点zk为为:一、算法的基本原理一、算法的基本原理46n M为所要分析的复频谱的点数,不一为所要分析的复频谱的点数,不一定等于定等于N,A和和W都是任意复数都是任意复数:47n当当 各各zk就均匀等间隔地就均匀等间隔地分布在单位圆上,这时就是求序列的分布在单位圆上,这时就是求序列的DFT。48nCZT的快速算法的快速算法n直接计算这一公式,与直接计算直接计算这一公式,与直接计算DFT相相似,总共算出似,总共算出M个抽样点,需要个抽样点,需要NM次复次复数乘法数乘法与(与(N-1)M次复数加法,当次复数加法,当N,M
18、较大时,运算量将很大。较大时,运算量将很大。49n采用采用布鲁斯坦布鲁斯坦(Bluestein)提出的等式,)提出的等式,可以将以上运算可以将以上运算转换为卷积和转换为卷积和形式,进而形式,进而采用采用FFT算法,提高运算速度。算法,提高运算速度。n布鲁斯坦所提出的等式为:布鲁斯坦所提出的等式为:n 由此可得:由此可得:5051nCZT变换可以通过求变换可以通过求g(n)与与h(n)的线性的线性卷积,然后乘上卷积,然后乘上1/h(n)而得到,即:而得到,即:nh(n)为角为角频率频率随时间随时间线性增长线性增长的复指数序列,的复指数序列,称为线性调频信号(称为线性调频信号(Chirp sign
19、al),因此,),因此,称为线调频称为线调频z变换。变换。 52n线性系统线性系统h(n)是非因果是非因果的。的。n当当n取值为取值为0到到N-1,k取值为取值为0,1,M-1时,则时,则 h(n)的定义区间为的定义区间为 ,可以看作有限长序列可以看作有限长序列二、二、Chirp-Z变换的快速计算变换的快速计算53n由于输入信号由于输入信号g(n)也是长度为也是长度为N有限长序有限长序列,所以列,所以g(n)*h(n)的点数为的点数为2N+M-2n因而用因而用圆周卷积代替线性卷积圆周卷积代替线性卷积不产生混叠不产生混叠失真的条件是圆周卷积的点数(周期)应失真的条件是圆周卷积的点数(周期)应大于
20、或等于大于或等于2N+M-2n但是我们只需要前但是我们只需要前M个值个值 ,对其他值是否,对其他值是否有混叠失真并不感兴趣,这样有可能将有混叠失真并不感兴趣,这样有可能将圆圆周卷积的点数缩减到最小为周卷积的点数缩减到最小为N+M-1。54g(n)h(n)5556CZT快速运算的实现步骤:n(1)选择一个最小的整数L,使其满足LN+M-1,同时以便采用基-2 FFT;n(2)将 g(n)补零成为L点的序列; 利用FFT法求此序列的L点DFT: 57n(3)形成L点h(n)序列n n 利用FFT法求其L点的DFT,n (4)将H(r)和G(r)相乘,得Q(r)=H(r)G(r), n (5)用FFT法求Q(r)的L点IDFT,其中前M个值为CZT数值,nM值对我们没有意义。n (6)最后求X(zk):58n9.19n9.21n9.28n9.48