数字信号处理05.docx
- 文档编号:12100640
- 上传时间:2023-04-17
- 格式:DOCX
- 页数:23
- 大小:438.98KB
数字信号处理05.docx
《数字信号处理05.docx》由会员分享,可在线阅读,更多相关《数字信号处理05.docx(23页珍藏版)》请在冰豆网上搜索。
数字信号处理05
第五章快速傅立叶算法
DFT与卷积运算是信号处理中两个最基本、最常用的运算,而实际上两者可以相互转化,不但卷积可以化为DFT处理,相关、滤波、谱分析都可以化为DFT来实现。
对于N点序列x(n),其DFT变换对定义为:
其中,
,因此,求N点X(k)需要N2次复数乘法,N(N-1)次复数加法,而实现一次复数乘法需要4次实数乘法和2次实数加法,而实现一次复数加法,则需要两次实数加法。
当N较大时,计算量相当可观。
如N=1024,则需要1048576次复数乘法,也就是4194304次实数乘法,因此这种变换在时间和操作上是极不现实的,尤其在遇到2D图象处理时,计算量实在大得惊人。
;
但考察W矩阵,可以看出,DFT运算中其实包含了大量的重复计算,虽然其中有N2个元素,但其中只有N个是独立的。
即:
;且这N个值也有一定对称关系,这种周期性和对称性可以描述为:
如,对4点DFT,按公式直接计算需42=16次复数乘法,按上述周期性和对称性,可得:
;
将该矩阵第二列和第三列交换,得:
;
由此得:
;
这样,求四点DFT实际上只需要一次复数乘法即可。
问题的关键是如何利用W因子的周期性和对称性,导出一种高效的快速算法。
1965年J.W.Cooley和J.W.Tukey提出了快速傅立叶变换算法(FASTFOURIERTRANSFORM,FFT),使N点DFT乘法计算量由N2次下降为(Nlog2N)/2次,仍以N=1024为例,计算量降为5120次,仅为原来的4.88%,这一发现在数字信号处理领域具有里程碑式的意义,并且带动了其它新算法的不断涌现,使数字信号处理广泛应用于众多的技术领域。
近三十多年来,快速傅立叶变换的发展主要有两个方向:
一类是针对N等于2的整数次幂的算法,如基2算法、基4算法、实因子算法、分裂基算法等,可以证明前面提到的四点DFT可以不用乘法,而只用加法来实现,因此基4算法比基2算法更有效。
1984提出的分裂基(split-radix)算法同时使用基2和基4算法,是目前对N等于2的整数次幂的各类算法中最为理想的一种。
另一类是N不等于2得整数次幂的算法,它以Winograd(WFTA)算法为代表,是建立在下标映射和数论的基础上的一套全新的算法,WFTA算法基于中国剩余定理,素因子算法基于素数定理,其运算速度比分裂基算法还要快,但理论复杂、编程困难、数据长度也受到较大限制,且在程序运行中,数据所占内存和数据传递次数比Cooley-Tukey增加了很多,因此在技术实现上是否真的具有突出的优越性,值得怀疑,因此,在我们的课程里,对WFTA不做讨论。
§5.1基2FFT算法
1.时间抽取(DIT)基2FFT算法
(1)算法的推导
对于(2-5-1):
令N=2M,M为正整数,如果将x(n)按奇偶分成两组,即令n=2r,n=2r+1,r=0,1,…,N/2-1,于是):
其中:
;
令
,
则:
;
A(k),B(k)都是N/2点的DFT,X(k)是N点DFT,因此单用上式表示X(k)是不完全的,但:
;
这样用A(k),B(k)可完整地表示X(k)。
当N=8时,A(k),B(k)和X(k)关系如图:
A(k),B(k)仍是高复合数(N/2)的DFT,我们按照上述方法继续给以分解,令r=2l,r=2l+1,l=0,1,…,N/4-1,则A(k)和B(k)可分别表示为:
;
令:
;
那么:
;
同理,令:
;
那么:
;
若N=8,这时C(k),D(k),E(k),F(k)都是二点的DFT,无需再分,即:
C(0)=x(0)+x(4),E(0)=x
(1)+x(5),C
(1)=x(0)-x(4),E
(1)=x
(1)-x(5),
D(0)=x
(2)+x(6),F(0)=x(3)+x(7),D
(1)=x
(2)-x(6),F
(1)=x(3)-x(7);
若N=16,32,或2的更高次幂,可按这样的方法继续分下去,直到两点的DFT为止。
在以上算法中将时间下标n按奇偶分开,故称为时间抽取算法(DecimationinTime,DIT),计算过程可表述如下图:
基本运算单元如图所示:
2、算法的讨论
下面我们对上面的推导过程进行讨论,来找出FFT算法的一般性的规律。
(一)“级”的概念
上述过程中,将N点DFT先分成两个N/2点的DFT,再是4个N/4点DFT,进而八个N/8点DFT,直到N/2个两点DFT。
每分一次,称为一“级”运算。
因为M=log2N,所以,N点DFT可分成M级,如图,图中从左至右,依次为m=0级,m=1级,…,m=M-1级。
(二)蝶形单元
在算法的信号流图中存在大量几何形状像蝴蝶的运算结构,称为蝶形运算单元,在第m级有:
;
p,q是参与该蝶形单元运算的上、下节点的序号。
很明显,第m级序号为p,q的两个点只参与该蝶形单元的运算,其输出在第m+1级。
且这一蝶形单元也不再涉及别的点。
因此,在进行编程时,可以将蝶形单元的输出仍放在输入数组中,这一特点称为“同址运算”。
由于每一级都含有N/2个蝶形单元,每个蝶形单元又只需要一次复数乘法,两次复数加法,那么,完成M=log2N级共需要的复数乘法次数Mc和复数加法次数Ma分别为:
;
将算法的信号流图改写为下图:
在该图中,每一级的数据自上而下都按自然顺序排列,在第m级,上、下节点p,q之间的距离为:
q-p=2m
(三)“组”的概念
由上图可以看出,每一级的N/2个蝶形单元可以分为若干组,每一组有着相同的结构和Wr因子分布。
如m=0级分成了4组,m=1级分成了两组,m=M-1级分成了一组,因此,第m级的组数是N/2m+1,m=0,1,…,M-1。
(四)Wr因子分布
考察(2-6-1)到(2-6-10),可以看出,第一次将N点DFT分成两个N/2点DFT时,相当于上图中最后一级,这时出现的Wr因子是:
,r=0,1,…,N/2-1。
再往下分时,依次是
,因此,每一级Wr因子分布的规律是:
;
因此,不难得出Wr因子分布的一般规律:
第m级,
,r=0,1,…,2m-1。
(五)码位倒置
从上图还可以看出,变换后的序列X(k)按正序排列,但输入序列x(n)的次序不再是原来的自然顺序,这正是由于对x(n)作奇、偶分开所产生的。
对N=8,其自然序号是0,1,2,3,4,5,6,7;第一次奇、偶分开,得到两组N/2点DFT,x(n)的序号是:
0,2,4,6,||1,3,5,7;
对每一组再作奇、偶分开,这时应视每组仍按自然顺序排列,故抽取后得4组,每组的序号是:
0,4,||2,6,||1,5,||3,7;
这一顺序正好是上图中输入端序列x(n)的排列次序,掌握了这一规律,对N为2的更高次幂,也都可以得到正确的抽样次序。
如果我们将x(n)的序号n==0,1,…,N-1写成二进制,如N=8,x(0),…,x(7)对应的是:
x(000),x(001),x(010),x(011),x(100),x(101),x(110),x(111),
将二进制数码翻转,得:
x(000),x(100),x(010),x(110),x(001),x(101),x(011),x(111),它们对应的十进制序号分别是:
x(0),x(4),x
(2),x(6),x
(1),x(5),x(3),x(7),也正好是按奇偶抽取所得到的顺序,掌握这一规律,我们就可以正确、快捷地编程。
§5.2频率抽取(DIF)基2FFT算法
和DIT相对应,DIF算法是将频域X(k)的序号k按奇、偶分开,对DFT,先将x(n)按序号分成上、下两部分,得:
;其中
;分别令k=2r,k=2r+1,r=0,1,…,N/2-1,得:
;
令
;
则:
;
这样,将一个N点DFT分成两个N/2点的DFT,分的方法是将X(k)按序号k的奇偶分开,直到得出两点的DFT,如图所示16点DIF算法流图,可以看出,输入是正序的,输出是按奇、偶分开的倒序。
但不论是DIF还是DIT,都可以看作N×N的W矩阵作分解来实现。
§5.3进一步减少运算量的措施
1)多类蝶形单元运算
对N=2M,共需进行M级运算,每级N/2个蝶形单元,每个蝶形单元需要一次复数乘法,所以总共需要MN/2次复数乘法,由Wr因子分布的一般规律:
第m级,
,r=0,1,…,2m-1。
当m=0时,即对第零级,所有的W因子的指数全为零,所以Wr=1,这一级不需要乘法。
对m=1级,Wr=1,或Wr=-j。
我们知道,两个复数相乘时,若一个为纯虚数,则也不需要乘法。
在DFT中,Wr又称“旋转因子”(twiddlefactor),像
这样的旋转因子又称平凡的旋转因子(trivialtwiddlefactor)去掉前两级后,所需的复数乘法应是:
Mc=N(M-2)/2
进一步分析,在m=2级,每一组含有
这两个平凡的旋转因子,这一级共有(N+22+1)=N/8组,故这一级平凡的旋转因子数为N/4个。
依次类推,m=3时有N/8个,最后一级,即m=M-1时,有N/2M-1个。
这样,从m=2到最后一级共有:
个平凡的旋转因子,因此Mc应改为:
Mc=N(M-3)/2+2
同样,所需要的复数加法量也不是Ma=MN,而是:
Ac=3N(M-1)/2+2
前面已经谈到,实现一次复数乘法需要4次实数乘法,两次实数加法,但对于:
这样的特殊复数,因为:
;
所以,可两次实数乘法、两次实数加法来实现。
从m=2到m=M-1级,每一级都包含不同数目的
这样的因子,数量多少,也如上计算。
这样,为完成N=2M点DFT,所需要的实数乘法数:
MR=4[N(M-3)/2+2-(N/2-2)]+2(N/2-2)=N(2M-7)+12;
同理可导出所需的实数加法数是:
AR=3N(M-1)+4;
以上两式一般是比较各种算法的基础。
一个旋转因子对应一个蝶形单元,若在程序中包含了所有的旋转因子,则称该算法为含一类蝶形单元,若去掉Wr=±1,则说含二类蝶形单元,若再去掉Wr=±j,则称含三类蝶形单元,如果再特殊处理
这样的蝶形单元,我们则称该算法含有4类蝶形单元,下表给出了N在不同值时取不同类型蝶形单元所需要的实数乘法和实数加法,尽管该方法无根本的突破,但在N较大时,乘法的节约也是相当可观的。
如N=4096,用三类蝶形单元时的乘法数比仅用一类时节约25%,当然蝶形单元用得越多,编程就越复杂。
一类蝶形单元
二类蝶形单元
三类蝶形单元
四类蝶形单元
N
M1
A1
M2
A2
M3
A3
M4
A4
2
4
6
0
4
0
4
0
4
8
48
72
20
58
8
52
4
52
32
320
480
196
418
136
388
108
288
128
1792
2688
1284
2434
1032
2308
908
2308
512
9216
13824
7172
12802
6152
12292
5644
12292
2048
45056
67584
36868
63490
32776
61444
30732
61444
4096
98304
14756
81924
139266
73736
135172
69644
135172
2)W因子的生成
在FFT中,乘法主要来自旋转因子,因为
,
所以在对Wr相乘时,必须产生相应的正、余弦函数。
在编程时,正、余弦函数的产生一般有两个办法,一是在每一步直接产生,二是在程序开始前预先计算出Wr,将r=0,1,…,N-1这N个独立的值存于数组中,等效于建立了一个正、余弦函数“表”。
在程序执行时可直接“查表”得到。
这样提高了运算速度,但要占较多的内存。
3)实输入数据时的FFT算法
在实际工作中,输入数据x(n)一般都是实序列。
如果不采取特殊措施,往往是把x(n)视为一个虚部为零的复序列,这增加了运算的时间。
处理实数据输入问题有着不同的方法。
实际上,当每一位研究者提出一种新算法时,都相应地讨论在该算法下实数据输入问题。
但DFT的核函数毕竟是复序列,因此尚没有有效的方法来解决这一问题,最早提出的方法是用N点FFT同时计算两个N点实序列的DFT,一个作为实部,另一个作为虚部,计算完后再把输出按奇、偶、虚、实特性加以分离,另一个方法就是用N/2点FFT计算一个N点序列的DFT,将该序列的偶序号置为实部,奇序号置为虚数,同样在最后使之分离。
理论上讲,这样可以减少一半的计算量。
§5.4分裂基算法
分裂基算法也称基2/4算法或混合基算法。
1)频率抽取基4FFT算法
令N=4M,对N点DFT可按如下方法作频率抽取:
;
即
分别令:
k=4r,k=4r+2,k=4r+1以及k=4r+3,r=0,1,…,N/4-1,有:
若N=16,通过上述推导即可把一个16点的DFT分成4个4点的DFT,其信号流图如所示,该图分m=0,m=1两级,最右边的4个4点的DFT,每一个都是基4FFT的基本单元,如图中方框内所示。
可以看出,基4FFT的基本单元仅有一个纯虚数j需要做乘法。
由于基4算法使做FFT的级数减少一半,故所需乘法数量也相应减少。
仿照基2算法,可得到使用四类蝶形单元时所需的实数乘法和实数加法的数量:
MR=[3N/2]log2N-5N+8;AR=[11N/4]log2N-13N/6+8/3
2)分裂基算法
仔细观察基2频率抽取算法,可以发现在每一级中每一组的上半部的输出都不乘以旋转因子,它们对应偶序数的输出,而旋转因子都出现在奇序号上,这一点也可以从下式看出来:
;
从前面我们知道,基4算法比基二算法更有效,因此可以将二者结合起来,这就是在1984年第一次被提出来的“分裂基”算法。
DuhamelP,HoltmannH.Split-radixFFTalgorithm.ElectronicsLetters.1984,20
(1)14-16
该算法的基本思路是对偶序号输出使用基2算法,对奇序号输出使用基4算法。
由于分裂基算法在目前已知的所有针对N=2M的算法中具有最少的乘法次数和加法次数,并且具有和Cooley-Tukey算法同样好的结构,因此被认为是最好的快速傅立叶变换算法。
现在的研究表明,该算法也最接近理论上所需乘法次数的最小值。
a)算法的推导
对N=2M点DFT,重写上式的偶序号输出项,即:
对k的奇序号项用基4算法,即:
式中r=0,1,…,N/4-1。
上面三式构成了分裂基算法的L型算法结构。
如下图所示:
现在以N=16为例,推导其算法,并给出信号流图:
令:
a(n)=x(n)+x(n+8)n=0,1,2,…,7
b(n)=x(n)-x(n+8)n=0,1,2,3
c(n)=x(n+4)-x(n+12)n=0,1,2,3
;由奇、偶序号输出式,有:
;
其中后两式各是四点DFT,不需再分,对第一式,可继续做分裂基算法。
因为:
;所以,令r=2l,r=4l+1,r=4l+3,得:
;式中:
f(n)=a(n)+a(n+4)n=0,1,2,3
g(n)=a(n)-a(n+4)n=0,1h(n)=a(n+2)-a(n+6)n=0,1
;由此可得16点的分裂基算法的信号流图。
b)分裂基算法的计算量
由式(2-6-21)、(2-6-22)可以看出,一个N点DFT在第一级被分成了一个N/2点DFT和两个N/4点的DFT。
N/2点DFT对应偶序号输出,不包括W因子,两个N/4点DFT对应奇序号输出,共有N/2个旋转因子。
其中包含有两个W0因子,两个
因子,它们可以特殊处理,因此,这一级应需要(N/2-4)个一般复数乘法和两个乘以的特殊复数乘法。
若实现一次复数乘法需要四次实数乘法、两次实数加法,则实现这一级运算共需要[4(N/2-4)+2×2]=2N-12次实数乘法。
由此可得到递推公式:
Qn=Qn-1+2Qn-2+2×2n-12
式中n=3,4,…,M,M=log2N,Qn代表N=2n时所需的乘法量,初始条件是Q1=0,Q2=0,现在来求解这一差分方程,不妨将上式写成:
x(n)=x(n-1)+2x(n-2)+2n+1-12
假定N→∞,则M也为无穷,两边取Z变换(注意n值从3开始):
整理得:
,作部分分式分解,并求逆Z变换,注意到n≥3,可得:
,当n=M时,有:
,这就是分裂基算法在四类蝶形单元情况下所需的实数乘法次数的计算公式,同理可推出所需实数加法的计算公式为:
。
下表给出了N取不同值时使用四类蝶形单元时基2、基4和分裂基算法所需实数乘法和实数加法的数量,从表中可以看出,当N≥64时,分裂基算法所需的实数乘法约是基2算法的2/3。
N
基2
基4
分裂基
MR
AR
MR
AR
MR
AR
2
0
4
0
4
4
0
16
0
16
0
16
8
4
52
4
52
16
28
148
24
144
24
144
64
332
964
264
920
248
912
256
2316
5380
1800
5080
1656
5008
1024
13324
27652
10248
25944
9336
25488
4096
69644
135172
53256
126926
48248
123792
c)算法的特点
从前面的介绍可知,分裂基算法具有以下优点:
1在已知的N=2M的算法中,所需乘法数最少,并接近于理论上的最小值。
理论上已经证明,对于长度为N=2M的实序列x(n),其DFT所需的理论上的最小实数乘法次数为:
MR=2N-M2-M-2;但实际上要达到这一理论下界将需要过多的加法并且算法也将十分复杂,因而难以实现。
2具有与基2、基4算法一样的规则结构,可以进行同址运算,这在用IC芯片来实现这些算法时特别重要。
③
3分裂基算法合理地安排算法结构,使平凡的旋转因子最大程度地减少了。
§5.5输入\输出端仅取少数点的FFT算法(FFTPruning)
我们知道,DFT在时域和频域都是N点的周期序列,当数据较短时,有时希望在数据后面补零以提高频域的分点数,当然这要增加计算量。
当数据足够长时,我们有时并不需要得到频域的全部分点,而只要得到某个频带内的值即可,如对应低通或带通的值。
特别是当信号为一个或几个正弦信号时,因为是线谱,我们总希望在正弦信号频率处给以“细化”,而在其它处可以不取或分点较稀疏。
这样,在用FFT解决问题时,在输入端和输出端都可以简化。
近一步降低计算量。
1原始输入数据中含有较多零时的FFT算法
输入端作简化的频率抽取算法是在1971年首次被提出的。
如图所示:
该图是一个N=16点的DIFFFT信号流图,在输入端仅有两个点不为零,其余14点均为零。
按普通的FFT算法,每级应有8个蝶形单元,且共有4级,但仔细观察会发现,在第一级仅有两个蝶形单元,第二级有4个,第三、第四级都有8个。
但前三级都是形如下图的“半蝶形单元”。
我们知道计算一个完整的蝶形单元需要一次复数乘法,两次复数加法,而计算一次复数乘法都需要四次实数乘法和两次的实数加法,而计算一次复数加法需要两次实数加法。
所以,计算一个完整的蝶形单元需要四次实数乘法加六次实数加法,而计算一个半蝶形单元则需要四次实数乘法和两次实数加法。
一般说来,若N=2M,那么FFT运算共有M级,若x(n)中仅有NL=2L个非零值点。
那么零值点的个数等于:
2M-2L=2L(2M-L-1)个。
在FFT的流图中,可以简化运算的有M-L级,不能简化的有L级。
当L=1,M=4时,NL=2,那么有三级可以简化运算,一级不能简化。
在第一级,有2L个半蝶形单元,在第二级则为2L+1个,直到第M-L级,有2M-L个半蝶形单元。
令ta表示实现两次实数加法所需的时间,tm表示实现两次实数乘法所需的时间,
那么要实现一个半蝶形单元的时间为:
(ta+2tm),实现一个全蝶形单元所需的时间为
(3ta+2tm)。
对前M-L级的半蝶形单元,所需计算时间是:
;后L级为全蝶形单元所需的计算时间是:
L(3ta+2tm)2M-1
如不使用FFTPruning,那么M级全蝶形单元所需时间为:
M(3ta+2tm)2M-1。
令tf是使用FFTPruning和不使用FFTPruning所需时间之比,则:
,大部分计算机有:
tm>>ta,则上式可变为:
;其中Δ=M-L,如果我们把半蝶形单元也视为全蝶形单元来看待,那么节约时间的多少便取决于蝶形单元的减少量,于是:
;即:
;
若M=4,L=1,则tf=0.6875,这样使用FFTPruning可节约31%的计算时间。
观察半蝶形单元结构,若是DIT算法,则Wr因子应乘在输入端的下部,因为xm(q)为零,所以在时间抽取的DIT算法中,半蝶形单元既不需要乘法,也不需要加法,令N=2M,若有NL=2L个非零点,则将有M-L级可以将前一级的数据简单地复制到下一级相应的蝶形单元的输入端,这样又进一步节约了计算时间,在N=16,L=2时,算法的信号流图如下:
现在来看一下所需的时间。
因为非零点为2L个,故前M-L级可以由前一级复制过来,不需要乘法运算,这样实际所需的相对时间比为:
tt=L/M,节约的时间为:
tT=1-tt=(M-L)/M,对M=4,L=1,所需的相对时间为0.25,节约的时间为:
75%。
对DIFPruning,由前面tf的式子可知:
;
即
;因此DITPruning可以比DIFPruning节约更多的时间。
2输入输出端同时使用FFTPruning算法
我们知道,DIF和DIT的FFT信号流图在结构上是互补的,差别只是W因子所处的位置。
当输出端仅需少数点时,也可以将前面的算法用于输出端。
当输入、输出端都只有少数点时,自然可以在输入、输出端同时使用Pruning技术。
例如,在输入端使用DIFPruning,那么在输出端使用DITPruning。
设N=2M,输入端非零数据长度NL=2L,输出端所需的变换长度NF=2LF。
将输入端使用DITPruning,由此节约的时间tT=1-tt=(M-L)/M。
对于剩下的L级,若L+LF>M,那么DIFPruning可应用到最后的M-LF级,若L+LF<M,则可应用到最后的L级。
这两种情况由输出DIFPruning节约的时间分别是:
tF1=[M-LF-2(1-2LF-M)]/ML+LF≥M;tF2=2L+1-M(2LF-1)/ML+LF<M
因此,总的节约时间应为:
tT+tF1或tT+tF2,取决于L+LF是大于还是小于M。
例如,当M=4,L=2,LF=2时,总的相对节约时间为0.625,即节约了62.5%。
上面讨论的算法,在输入端x(n)的非零值序号n是从0到2L-1,从2L补零到N-1。
输出端X(k)的非零值序号也是从0开始,到2LF-1,这是对应于低通的情况,当输入端希望是窄带谱时,上述算法还需要进一步改进。
当N=16,L=3,输出仅需要两个点,即LF=1,一般要输出的X(k)的序号是从k1到k2,k1,k2由使用者指定,但两者要满足:
k2-k1=2LF-1。
由上面的讨论可知:
FFTPruni
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数字信号 处理 05