引言 1.DFT是信号分析与处理中的一种重要变换。 2.1965年,Cooley, Tukey《机器计算傅里叶级数的一种算法》 4.本章主要讨论基2FFT算法。
第4章学习目标 1.理解按时间抽取的基-2FFT算法的算法原理、运算流图、所需计算量和算法特点 3.理解IFFT算法
第4章 快速傅里叶变换(FFT) 4.1 引言 4.2 基2FFT算法 4.3 进一步减少运算量的措施 4.4 分裂基FFT算法 4.1 引言 4.2 基2FFT算法 4.3 进一步减少运算量的措施 4.4 分裂基FFT算法 4.5 离散哈特莱变换(DHT)
4.2 基2FFT算法 一. 直接计算DFT的特点及减少运算量的基本途径 1.对N点有限长序列直接计算DFT的运算量
运算量 复数乘法 复数加法 一个X(k) N N – 1 N个X(k) (N点DFT) N 2 N (N – 1) 实数乘法 实数加法 一次复乘 4 2 一次复加 一个 X(k) 4N 2N+2 (N – 1)=2 (2N – 1) N个X(k) (N点DFT) 4N 2 2N (2N – 1)
2.减少DFT运算量的基本途径 (1)把长序列分解成短序列计算,然后再组合出结果。
FFT算法分类: 时间抽取法 DIT: Decimation-In-Time 频率抽取法 DIF: Decimation-In-Frequency
二.时域抽取法(库利-图基算法) 1、算法原理 设序列点数 N = 2M,M 为整数。 若不满足,则补零 N为2的整数幂的FFT算法称基-2FFT算法。 (1)第一次分解 将序列x(n)按n的奇偶分成两组:
其中X1(k)和X2(k)分别为x1(r)和x2(r)的N/2点DFT 则x(n)的DFT: 其中X1(k)和X2(k)分别为x1(r)和x2(r)的N/2点DFT
再利用周期性求 的后半部分 先乘因子再加减 图4.2.1 蝶形运算符号
图4.2.2 N点DFT的一次时域抽取分解图(N=8)
分解后的运算量: 复数乘法 复数加法 一个N / 2点DFT (N / 2)2 N / 2 (N / 2 –1) 两个N / 2点DFT 一个蝶形 1 2 N / 2个蝶形 N / 2 N 总计 运算量减少了近一半
(2)第二次分解 N / 2仍为偶数,进一步分解:N / 2 N / 4
由X3(k)和X4(k)的周期性和 的对称性
同理: 其中: 这样逐级分解,直到2点DFT
2.基2时间抽取FFT算法流图 N=2 x[n]={x[0], x[1]}
4点基2时间抽取FFT算法流图 X1[0] x[0] X [0] 2点DFT X1[1] x[2] X [1] -1 X2[0] x[1]
4点基2时间抽取FFT算法流图
8点基2时间抽取FFT算法流图 x[0] X [0] x[2] X [1] 4点DFT x[4] X [2] X [3] x[6] x[1] -1 X2[1] X [5] -1 X2[2] X [6] -1 X2[3] X [7] -1
8点基2时间抽取FFT算法流图 x[0] X [0] x[2] X [1] 4点DFT x[4] X [2] X [3] x[6] x[1] -1 X2[1] X [5] -1 X2[2] X [6] -1 X2[3] X [7] -1
基2时间抽取FFT算法 第三级 第一级 第二级
三 DIT―FFT算法与直接计算DFT运算量的比较 每一级运算都需要N/2次复数乘和N次复数加(每个蝶形需要两次复数加法)。所以,M级运算总共需要的复数乘次数为 复数加次数为
2. FFT算法与直接计算DFT运算量的比较曲线 复乘次数 N N 2 图4.2.5 FFT算法与直接计算DFT所需乘法次数的比较曲线
四 DIT―FFT的运算规律及编程思想 1.原位计算 2.旋转因子的变化规律 中间数据的存储,可采用原位存储法。即每次蝶形运算的结果可以存储在原数据的同一个存储单元。这样在高速硬件实现时,可节省存储器。 以蝶形运算为基础进行组合计算,旋转因子WNp的指数p与运算所在的级数和组内位置有关。如下所示: 继续
继续
第L级共有2L-1个不同的旋转因子。N=23=8时的各级旋转因子表示如下: L=1时, , J=0 L=2时, , J=0,1 对N=2M的一般情况,第L级的旋转因子为 由于 所以 继续
3. 蝶形运算规律 第m级运算每个蝶形的两节点距离为 2m–1 第m级运算: 其中,m表示第m列,且m =1,… ,L 例如N=8=23 ,第一级(列)距离为21-1=1, 第二级(列)距离为22-1=2, 第三级(列)距离为23-1=4。
4. 序列的倒序 输入序列的混序。因为DFT输入序列是顺序采样的,所以在计算FFT 之前需要进行序列按混序要求排序。排序算法很多,较常用的计算混序号的方法有二进制序号反转算法。
k k k 1 x [00 0] 1 1] x [ k 0] x [10 0] 1 x [01 0] x [11 0] x [ k ] 1 1 2 1 x [00 0] 1 1] 1 2 x [ k 0] x [10 0] 1 x [01 0] x [11 0] x [ k 2 1 ] 1 x [001 ] 1 x [101 ] 1 x [01 1] x [111 ]
n0 n1 n2 1 将原有二进制数倒序,即原来的高位变为低位,低位变高位 倒位序 自然序 000 100 4 1 001 010 2 100 4 1 001 010 2 110 6 3 011 101 5 111 7 n0 n1 n2 1
x(0) x(4) x(2) x(6) x(1) x(5) x(3) x(7) 倒序处理方法 存储单元 A(1) A(2) A(3) A(4) A(5) A(6) A(7) A(8) x(0) x(1) x(2) x(3) x(4) x(5) x(6) x(7) x(0) x(4) x(2) x(6) x(1) x(5) x(3) x(7) 自然顺序 变址 倒位序
五、按频率抽选的基-2FFT算法 1、算法原理 设序列点数N=2M,M为整数。 (1)将X(k)按k的奇偶分组前,先将输入x(n)按n的顺序分成前后两半: 常系数线性差分方程
(2)按k的奇偶将X(k)分成两部分: 当k取偶数(k=2r,r=0,1,…,N/2-1)时
则X(2r)和X(2r+1)分别是x1(n)和x2(n)的 N / 2点DFT,记为X1(k)和X2(k) 令 则X(2r)和X(2r+1)分别是x1(n)和x2(n)的 N / 2点DFT,记为X1(k)和X2(k)
该蝶形运算的特点:先加减后乘因子
x1(0) x1(1) -1 x1(2) x1(3) x2(0) x2(1) x2(2) x2(3) N/2点 DFT x(0) x(7) x(1) x(2) x(3) x(4) x(5) x(6) X1(0)=X(0) X2(0)=X(1) X1(1)=X(2) X1(2)=X(4) X1(3)=X(6) X2(1)=X(3) X2(2)=X(5) X2(3)=X(7)
(3)N /2仍为偶数,进一步分解:N /2 N /4
x3(0) x3(1) -1 x4(0) x4(1) N/4点 DFT x1(0) x1(1) x1(2) x1(3) X3(0)=X1(0)=X(0) X4(0)=X1(1)=X(2) X3(1)=X1(2)=X(4) X4(1)=X1(3)=X(6)
同理: 其中:
(4)逐级分解,直到2点DFT 例如:当N=8时,即分解到x3(n),x4(n),x5(n),x6(n),n=0,1
2.基2频域抽取FFT算法流图 N=2 x[n]={x[0], x[1]}
4点基2频域抽取FFT算法流图 x1[0] x[0] X [0] 2点DFT x1[1] x[1] X [2] -1 x2[0] x[2]
4点基2时间抽取FFT算法流图
8点基2频域抽取FFT算法流图 x1(0) X1(0)=X(0) x(0) 4点 DFT X1(1)=X(2) x1(1) x(1) -1 x2(1) X2(1)=X(3) -1 x2(2) X2(2)=X(5) -1 x2(3) X2(3)=X(7) -1
8点基2频域抽取FFT算法流图 x1(0) X(0) x(0) x1(1) X(4) x(1) x1(2) X(2) x(2) x1(3) -1 x1(2) x1(3) x2(0) x2(1) x2(2) x2(3) x(0) x(1) x(2) x(3) x(7) x(4) x(5) x(6) x3[0] x3[1] x4[0] x4[1] 2点DFT -1 X(0) X(4) X(2) X(6) x5[0] x5[1] x6[0] x6[1] 2点DFT -1 4点 DFT X2(0)=X(1) X2(1)=X(3) X2(2)=X(5) X2(3)=X(7) X(1) X(5) X(3) X(7)
第一级 第二级 第三级
3. 计算量和算法特点 (1)中间数据的存储,可采用原位存储法。即每次碟形运算的结果可以存储在原数据的同一个存储单元。这样高速硬件实现时,可节省存储器。 (2)输出序列为混序。因为DFT输出序列要求是顺序的,所以在计算FFT 之后需要进行序列按混序要求排序。输出序列混序与时间抽取方法的输入序列混序相同。所以排序算法相同。 (3)以碟形运算为基础进行组合计算,计算因子Wr的指数r与运算所在的级数和组内位置有关。
蝶形运算 对N=2M点FFT,输入自然序,输出倒位序, 两节点距离:2M-m=N / 2m 第m级运算: 继续
即将地址k乘以2m-1。 蝶形运算两节点的第一个节点为k值,表示成L位二进制数,左移m-1位,把右边空出的位置补零,结果为r的二进制数。 继续
4、DIT与DIF的异同 基本蝶形不同 输入顺序不同 运算量相同 都可原位运算 DIT和DIF的基本蝶形互为转置 DIT: 先复乘后加减;
六 、IFFT算法 比较: IDFT: DFT: 常系数线性差分方程
图4.2.17 DIT―IFFT运算流图(防止溢出)
直接调用FFT子程序计算IFFT的方法: 共轭 FFT 乘1/ N
4.3 进一步减少运算量的措施 一、多类蝶形单元运算 1.无关紧要旋转因子 在DFT中,称其值为 的旋转因子为无关紧要旋转 4.3 进一步减少运算量的措施 一、多类蝶形单元运算 1.无关紧要旋转因子 在DFT中,称其值为 的旋转因子为无关紧要旋转 因子,例如:WN0,WNN/2,WNN/4等,这些因子计算过程 中,不需要复数乘法。 因此,DIT-FFT的复乘次数降至
2.特殊因子WNN/8 结论:对应的每个蝶形运算节 省两次实数乘 从实数运算考虑,计算N=2M点 DIT―FFT所需实数乘法次数为
3.蝶形运算的分类 一类蝶形运算;二类蝶形运算;三类蝶形运算; 四类蝶形运算
二、旋转因子的生成 在FFT运算中,旋转因子WmN=cos(2πm/N)-jsin(2πm/N),求正弦和余弦函数值的计算量是很大的。 编程时,产生旋转因子的方法直接影响运算速度。 共有两种产生方法: 1.在每级运算中直接产生 2.在FFT程序开始前预先计算,存放在数组中,作为旋转 因子表,在程序执行过程中直接查表得到所需旋转因子
三、实序列的FFT算法 设x(n)为N点实序列,取x(n)的偶数点和奇数点分别作为新构造序列y(n)的实部和虚部,即 对y(n)进行N/2点FFT,输出Y(k),则
根据DIT―FFT的思想及蝶形运算规则,可得到 由于x(n)为实序列,所以X(k)具有共轭对称性,X(k)的另外N/2点的值为
4.4 分裂基FFT算法 4.4.1 分裂基FFT算法原理 1.频率抽取基4FFT算法 令N= 4M,对N点DFT可按如下方法作频率抽取: 即
分别令k=4r, k=4r+2,k=4r+1及k=4r+3,而r=0,1,…..,N/4-1,于是
2.分裂基算法
图4.4.1 分裂基第一次分解L形流图
图4.4.2 分裂基FFT算法L形排列示意图与结构示意图 (a)分裂基FFT算法L形排列示意图; (b)分裂基FFT算法运算流图结构示意图
图4.4.6 16点分裂基FFT运算流图
3.各种算法运算量的比较
4.5 离散哈特莱变换(DHT) 4.5.1 离散哈特莱变换定义 设x(n),n=0,1,…,N-1,为一实序列,其DHT定义为 4.5.1 离散哈特莱变换定义 设x(n),n=0,1,…,N-1,为一实序列,其DHT定义为 式中,cas(α)=cosα+sinα。其逆变换(IDHT)为
逆变换证明如下: 三角函数性质
可得 =x(n)
4.5.2 DHT与DFT之间的关系 容易看出,DHT的核函数 是DFT的核函数 的实部与虚部之和。
将XH(k)分解为奇对称分量XHo(k)与偶对称分量XHe(k)之和 由DHT定义有
所以,x(n)的DFT可表示为 同理,x(n)的DHT可表示为 因此,已知x(n)的DHT,则DFT可用下式求得:
4.5.3 DHT的主要优点 (1)DHT是实值变换,在对实信号或实数据进行谱分析时避免了复数运算,从而提高了运算效率,相应的硬件也更简单、更经济; (2)DHT的正、反变换(除因子1/N外)具有相同的形式,因而,实现DHT的硬件或软件既能进行DHT,也能进行IDHT; (3)DHT与DFT间的关系简单,容易实现两种谱之间的相互转换。
4.5.4 DHT的性质 1. 线性性质 2. x(N-n)的DHT 其中,当k=0时,XH(N-k)=XH(N)=XH(0)。
证明: 由DHT定义 而
3. 循环移位性质 (4.5.16) (4.5.17) 证明: 由DHT定义有
4. 奇偶性 奇对称序列和偶对称序列的DHT仍然是奇对称序列或偶对称序列,即DHT不改变序列的奇偶性。 5.循环卷积定理
证明: 下面利用DFT的循环卷积定理和DHT与DFT之间的关系来证明 其中,X1(k)=DFT[x1(n)],X2(k)=DFT[x2(n)],根据DHT与DFT之间的关系,则有
将上面两式代入式(4.5.20)并整理后,得 所以式(4.5.18)成立。同理可证明式(4.5.19)亦成立。 当x1(n)或x2(n)是偶对称序列时,则由DHT的奇偶性有
仿照快速DFT的分解方法,可通过时域抽取或频域抽取的方式实现快速DHT。 4.5.5 DHT的快速算法(FHT) 1.基2DIT―FHT算法及运算流图 仿照快速DFT的分解方法,可通过时域抽取或频域抽取的方式实现快速DHT。 x(n)的N=2M点DHT如下式: (4.5.22) 对x(n)进行奇偶抽取 (4.5.23)
与DFT的时域抽取分解比较, 不是 一个指数函数,所以处理要比W(2r+1)kN麻烦一些。 根据三角公式有 (4.5.24) (4.5.25)
令X0H(k)=DHT[x0(n)],X1H(k)=DHT[x1(n)],并考虑DHT的周期性,(4.5.25)式可写成 为了使以下推导中公式简明,记 C(k)=cos (2π/N)k ,S(k)=sin (2π/N)/k 。将式(4.5.26)中的k分别取k,N/2+k,N/2-k和N-k四个值,并考虑X0H(k)和X1H(k)以N/2为周期,得到
(4.5.27)
当k=0时,式(4.5.27)中有重复,可单独写成 (4.5.28) 同理,在k=N/4时有 (4.5.29)
图4.5.1 基2DIT-FHT原理和哈特莱碟形
图4.5.2 4点DFHT蝶形图
图4.5.3 16点基2DIT―FHT流图
2. 基2DIT―FHT的运算量 观察图4.5.3可知,运算流图中都是实数运算。N=2M点FHT流图共分为M级。用L表示流图自右向左的蝶形级数。第L级乘法次数为2L(N/2L-2),但最后二级无乘法运算,所以,总的乘法次数MH为 (4.5.30) 同理可求得总的加法次数AH为 (4.5.31)
4.5.6 实信号的快速循环卷积 用DIT―FFT计算两个实信号x1(n)和x2(n)的循环卷积时,最直接的方法是将信号的虚部都置为零,再按复序列用FFT计算。显然这样做是很浪费的。现在我们可用基2DIT―FHT直接进行实正交变换来处理实信号的循环卷积问题。 计算式 (4.5.18)