详解FFT快速傅里叶变换FFT.docx
- 文档编号:17647647
- 上传时间:2023-07-27
- 格式:DOCX
- 页数:27
- 大小:61.55KB
详解FFT快速傅里叶变换FFT.docx
《详解FFT快速傅里叶变换FFT.docx》由会员分享,可在线阅读,更多相关《详解FFT快速傅里叶变换FFT.docx(27页珍藏版)》请在冰点文库上搜索。
详解FFT快速傅里叶变换FFT
第四章快速傅里叶变换
有限长序列可以通过离散傅里叶变换(DFT)将其频域也离散化成有限长序列.但其计算量太大,很难实时地处理问题,因此引出了快速傅里叶变换(FFT).1965年,Cooley和Tukey提出了计算离散傅里叶变换(DFT)的快速算法,将DFT的运算量减少了几个数量级。
从此,对快速傅里叶变换(FFT)算法的研究便不断深入,数字信号处理这门新兴学科也随FFT的出现和发展而迅速发展。
根据对序列分解与选取方法的不同而产生了FFT的多种算法,基本算法是基2DIT和基2DIF。
FFT在离散傅里叶反变换、线性卷积和线性相关等方面也有重要应用。
快速傅里叶变换(FFT)是计算离散傅里叶变换(DFT)的快速算法。
DFT的定义式为
kn
N−1
X(k)=∑x(n)WN
RN(k)
n=0
N
在所有复指数值Wkn的值全部已算好的情况下,要计算一个X(k)需要N
次复数乘法和N-1次复数加法。
算出全部N点X(k)共需N2次复数乘法和N(N−1)次复数加法。
即计算量是与N2成正比的。
FFT的基本思想:
将大点数的DFT分解为若干个小点数DFT的组合,从而减少运算量。
WN因子具有以下两个特性,可使DFT运算量尽量分解为小点数的DFT
运算:
(1)周期性:
(k+N)n
N
=Wkn
=W(n+N)k
(2)对称性:
W(k+N/2)=−Wk
NN
利用这两个性质,可以使DFT运算中有些项合并,以减少乘法次数。
例子:
求当N=4时,X
(2)的值
3
∑44444
X
(2)=
n=0
x(n)W2n=x(0)W0+x
(1)W2+x
(2)W4+x(3)W6
=[x(0)+x
(2)]W0+[x
(1)+x(3)]W2
(周期性)
4
4
=[x(0)+x
(2)]-[x
(1)+x(3)]W0
4
(对称性)
通过合并,使乘法次数由4次减少到1次,运算量减少。
FFT的算法形式有很多种,但基本上可以分为两大类:
按时间抽取
(DIT)和按频率抽取(DIF)。
4.1按时间抽取(DIT)的FTT
为了将大点数的DFT分解为小点数的DFT运算,要求序列的长度N为复合数,最常用的是N=2M的情况(M为正整数)。
该情况下的变换称为基2FFT。
下面讨论基2情况的算法。
先将序列x(n)按奇偶项分解为两组
⎧x(2r)=x1(r)
⎨
⎩x(2r+1)=x2(r)
将DFT运算也相应分为两组
N
r=0,1,L,2−1
N−1
N
X(k)=DFT[x(n)]=∑x(n)Wkn
n=0
N−1
=∑x(n)Wkn+
N−1
∑x(n)Wkn
n=0
n为偶数
n=0
n为奇数
N/2−1
∑
(2)
2rk+
N/2−1
∑(2
+1)
(2r+1)k
=x
r=0
rWN
xrWN
r=0
N/2−1
∑
2rk+
N/2−1
k∑
2rk
=
r=0
x1(r)WN
WN
k
r=0
x2(r)WN
rk
N/2−1
=∑x1(r)W
r=0
N/2+WN
N/2−1
∑x2(r)W
r=0
rk
N/2
(因为W
2rk
N
rk
N/2
k
=X1(k)+WNX2(k)
其中X1(k)、X2(k)分别是x1(n)、x2(n)的N/2点的DFT
N/2−1
N/2−1
rkrk
X1(k)=∑x1(r)WN/2=∑x(2r)WN/2,0≤k≤−1
r=0
r=02
N/2−1
N/2−1
rk
rkN
X2(k)=∑x2(r)WN/2=∑x(2r+1)WN/2,0≤k≤−1
r=0
r=02
至此,一个N点DFT被分解为两个N/2点的DFT。
上面是否将全部N点的X(k)求解出来了?
分析:
X1(k)和
N
X2(k)只有N/2个点(k=0,1,L,
2
−1),则由
1N2
X(k)=X(k)+WkX(k)只能求出X(k)的前N/2个点的DFT,要求出
全部N点的X(k),需要找出X1(k)、X2(k)和X(k+N/2)的关系,其
N
中k=0,1,L,2−1。
由式子X(k)=X1
N
(k)+WkX
2(k)可得
1N
X(k+N/2)=X(k+N/2)+Wk+N/2
k
X2(k+N/2)化简得
N
X(k+N/2)==X1(k)−WNX2(k),k=0,1,L,2−1
这样N点DFT可全部由下式确定出来:
1N2
⎧⎪X(k)=X(k)+WkX(k)
1N2
⎩⎪⎨X(k+N/2)=X(k)−WkX(k)
k=0,1,L
N−1
2
(*)
上式可用一个专用的碟形符号来表示,这个符号对应一次复乘和两次复加运算。
aa+Wkb
ba−Wkb
-1
N
图蝶形运算符号
2
通过这样的分解以后,每一个N/2点的DFT只需要(N)2=N
次复数乘
24
2
法,两个N/2点的DFT需要2(N)2=N
次复乘,再加上将两个N/2点
22
DFT合并成为N点DFT时有N/2次与W因子相乘,一共需要
2
N+N
22
N2
≈次复乘。
可见,通过这样的分解,运算量节省了近一半。
2
因为N=2M,N/2仍然是偶数,因此可以对两个N/2点的DFT再
分别作进一步的分解,将两个N/2点的DFT分解成两个N/4点的DFT。
例如对x1(r),可以在按其偶数部分及奇数部分进行分解:
⎧x1(2l)=x3(l)
⎨
⎩x1(2l+1)=x4(l)
则的运算可相应分为两组:
N
l=0,1,L,4−1
N/4−1
X1(k)=∑x1(2l)W
l=0
2lk
N/2
N/4−1
+∑x1(2l+1)W
l=0
(2l+1)k
N/2
N/4−1
∑
3()/4
k
N/2
N/4−1
∑
4()
lk
N/4
=x
l=0
lWNW
k
xlW
l=0
N
=X3(k)+WN/2X4(k)
N/2
将系数统一为以N为周期,即Wk
k=0,1,L,4−1
N
=W2k,可得
13N
⎧⎪X(k)=X(k)+W2k
⎨
X4(k)
2k
k=0,1,L
N−1
X1(k+N/4)=X3(k)−WN
X4(k)4
同样,对X2(k)也可进行类似的分解。
一直分解下去,最后是2点的
DFT,2点DFT的运算也可用碟形符号来表示。
这样,对于一个N=23=8
的DFT运算,其按时间抽取的分解过程及完整流图如下图所示。
这种方法,由于每一步分解都是按输入序列在时域上的次序是属于偶数
还是奇数来抽取的,故称为“时间抽取法”。
分析上面的流图,N=2M,一共要进行M次分解,构成了从x(n)到
X(k)的M级运算过程。
每一级运算都是由N/2个蝶形运算构成,因此每一级运算都需要N/2次复乘和N次复加,则按时间抽取的M级运算后总共需要
复数乘法次数:
mF
=N⋅M
2
=NlogN
22
复数加法次数:
aF
=N⋅M=Nlog2N
根据上面的流图,分析FFT算法的两个特点,它们对FFT的软硬件
构成产生很大的影响。
(1)原位运算也称为同址运算,当数据输入到存储器中以后,每一级运算的结果仍然存储在原来的存储器中,直到最后输出,中间无需其它的存储器。
根据运算流图分析原位运算是如何进行的。
原位运算的结构可以节省存储单元,降低设备成本。
(2)变址分析运算流图中的输入输出序列的顺序,输出按顺序,输入是“码位倒置”
的顺序。
见图。
自然顺序
二进制表示
码位倒置
码位倒置顺序
0
000
000
0
1
001
100
4
2
010
010
2
3
011
110
6
4
100
001
1
5
101
101
5
6
110
011
3
7
111
111
7
码位倒置顺序
01234567
X(0)X(4)X
(2)X(6)X
(1)X(5)X(3)X(7)码位倒置的变址处理
在实际运算中,直接将输入数据x(n)按码位倒置的顺序排好输入很不
方便,一般总是先按自然顺序输入存储单元,然后通过变址运算将自然顺序
的存储换成码位倒置顺序的存储,这样就可以进行FFT的原位运算。
变质的功能如图所示。
用软件实现是通用采用雷德(Rader)算法,算出I的倒序J以后立即将输入数据X(I)和X(J)对换。
尽管变址运算所占运算量的比例很小,但对某些高要求的应用(尤其在实时信号处理中),也可设法用适当的电路结构直接实现变址。
例如单片数字信号处理器TMS320C25就有专用于FFT的二进制码变址模式。
4.2按频率抽取(DIF)的FTT
除时间抽取法外,另外一种普遍使用的FFT结构是频率抽取法。
频率抽取法将输入序列不是按奇、偶分组,而是将N点DFT写成前后两部分:
N−1
N
X(k)=DFT[x(n)]=∑x(n)Wkn
n=0
(N/2)−1
N−1
=∑x(n)Wkn+∑x(n)Wkn
n=0
NN
n=N/2
N/2−1
∑()
nk+
N/2−1
∑(+
nk
/2)
(n+N/2)k
=
n=0
xnWN
xnNWN
n=0
(N/2)k
N/2−1
=∑[x(n)+WN
x(n+N/2)]WN
n=0
因为WN/2=−1,W(N/2)k
=(−1)k,k为偶数时(−1)k
=1,k为奇数时
NN
(−1)k
=−1,由此可将X(k)分解为偶数组和奇数组:
N/2−1
∑N
X(k)=
n=0
[x(n)+(−1)kx(n+N/2)]Wnk
N/2−1
∑N
X(2r)=
n=0
[x(n)+x(n+N/2)]W2nr
N/2−1
∑
N/2
=
n=0
[x(n)+x(n+N/2)]Wnr
N/2−1
∑N
X(2r+1)=
N/2−1
n=0
[x(n)−x(n+N/2)]W(2r+1)n
∑NN/2
=
N
n=0
[x(n)−x(n+N/2)]WnWnr
⎧x1(n)=x(n)+x(n+N/2)
令⎨
n=0,1,L,N/2−1
⎩x2
(n)=[x(n)−x(n+N/2)]Wn
这两个序列都是N/2点的序列,对应的是两个N/2点的DFT运算:
N/2−1
∑1
N/2
X(2r)=
n=0
x(n)]Wnr
N/2−1
∑2
N/2
X(2r+1)=
x
n=0
(n)Wrn
这样,同样是将一个N点的DFT分解为两个N/2点的DFT了。
频率抽选法
对应的碟形运算关系图如下:
aa+b
(a−b)Wn
bN
n
-1WN
对于N=8时频率抽取法的FFT流图如下:
这种分组的办法由于每次都是按输出X(k)在频域的顺序上是属于偶数还是
奇数来分组的,称为频率抽取法。
与前面按时间抽取的方法相比,相同点问题:
如何利用快速算法计算IDFT?
分析IDFT的公式:
x(n)=IDFT[X(k)]=
1N−1
N
比较DFT的公式:
Nk=0
N−1
N
X(k)=DFT[x(n)]=∑x(n)Wnk,k=0,1,L,N−1
n=0
得知可用两种方法来实现IDFT的快速算法:
(1)只要把DFT运算中的每
一个系数Wnk该为W−nk,并且最后再乘以常数1,就可以用时间抽取法
NNN
或频率抽取的FFT算法来直接计算IDFT。
这种方法需要对FFT的程序和参
数稍加改动才能实现。
(2)因为
x(n)=
1N−1
N
1{DFT[X*(k)]},n=0,1,L,N−1,也就
[∑X*(k)Wnk]*=
Nk=0N
是说,可先将X(k)取共轭变换,即将X(k)的虚部乘以-1,就可直接调用
FFT的程序,最后再对运算结果取一次共轭变换并乘以常数1/N即可得到x(n)的值。
这种方法中,FFT运算和IFFT运算都可以共用一个子程序块,在使用通用计算机或用硬件实现时比较方便。
4.1.3混合基FFT算法
以上讨论的是基2的FFT算法,即N=2M的情况,这种情况实际上使用得最多,这种FFT运算,程序简单,效率很高,用起来很方便。
另
外,在实际应用时,有限长序列的长度N到底是多少在很大程度时是由人
为因素确定的,因此,大多数场合人们可以将N选定为N=2M,从而可以直接调用以2为基数的FFT运算程序。
如果长度N不能认为确定,而N的数值又不是以2为基数的整数次
方,一般可有以下两种处理方法:
(1)将x(n)用补零的方法延长,使N增长到最邻近的一个
N=2M数值。
例如,N=30,可以在序列x(n)中补进
x(30)=x(31)=0两个零值点,使N=32。
如果计算FFT的目的是为了了解整个频谱,而不是特定频率点,则此法可行。
因为有限长序列补零以后并不影响其频谱
X(ejw),只是频谱的采样点数增加而已。
(2)如果要求特定频率点的频谱,则N不能改变。
如果N
为复合数,则可以用以任意数为基数的FFT算法来计算。
快速傅里叶变换的基本思想就是要将DFT的运算尽量分小。
例如,N=6时,可以按照N=3×2分解,将6点DFT分解为3组2点DFT。
举例:
N=9时的快速算法。
凡是可以利用傅里叶变换来进行分析、综合、变换的地方,都可以利
用FFT算法及运用数字计算技术加以实现。
FFT在数字通信、语音信号处理、图像处理、匹配滤波以及功率谱估计、仿真、系统分析等各个领域都得到了广泛的应用。
但不管FFT在哪里应用,一般都以卷积积分或相关积分的具体处理为依据,或者以用FFT作为连续傅里叶变换的近似为基础。
4.2.1利用FFT求线性卷积—快速卷积
在实际中常常遇到要求两个序列的线性卷积。
如一个信号序列x(n)
通过FIR滤波器时,其输出y(n)应是x(n)与h(n)的卷积:
y(n)=x(n)∗h(n)=
∞
∑x(m)h(n−m)
m=−∞
有限长序列x(n)与h(n)的卷积的结果y(n)也是一个有限长序列。
假设x(n)与
h(n)的长度分别为N1和N2,则y(n)的长度为N=N1+N2-1。
若通过补零使x(n)与h(n)都加长到N点,就可以用圆周卷积来计算线性卷积。
这样得到用FFT运算来求y(n)值(快速卷积)的步骤如下:
(1)对序列x(n)与h(n)补零至长为N,使N≥N1+N2-1,并且
N=2M(M为整数),即
⎨
x(n)=⎧x(n),n=0,1,L,N1−1
⎩0,
n=N1,N1+1,L,N−1
⎨
h(n)=⎧h(n),n=0,1,L,N2−1
⎩0,
n=N2,N2+1,L,N−1
(2)用FFT计算x(n)与h(n)的离散傅里叶变换
x(n)⇔(FFT)⇔X(k)
(N点)
h(n)⇔(FFT)⇔H(k)
(N点)
(3)计算Y(k)=X(k)H(k)
(4)用IFFT计算Y(k)的离散傅里叶反变换得:
y(n)=IFFT[Y(k)](N点)
互相关及自相关的运算已广泛的应用于信号分析与统计分析,应用于连
续时间系统也用于离散事件系统。
用FFT计算相关函数称为快速相关,它与快速卷积完全类似,不同的是一个应用离散相关定理,另一个应用离散卷积定理。
同样都要注意到离散傅里叶变换固有的周期性,也同样用补零的方法来绕过这个障碍。
设两个离散时间信号x(n)与y(n)为已知,离散互相关函数记作Rxy(n),定义为
Rxy(n)=
∞
∑x(m)y(n+m)
m=−∞
如果x(n)与y(n)的序列长度分别为N1和N2,则用FFT求相关的计算步骤
如下:
(1)对序列x(n)与y(n)补零至长为N,使N≥N1+N2-1,并且
N=2M(M为整数),即
⎨
x(n)=⎧x(n),n=0,1,L,N1−1
⎩0,
n=N1,N1+1,L,N−1
⎨
y(n)=⎧y(n),n=0,1,L,N2−1
⎩0,
n=N2,N2+1,L,N−1
(2)用FFT计算x(n)与y(n)的离散傅里叶变换
x(n)⇔(FFT)⇔X(k)
(N点)
y(n)⇔(FFT)⇔Y(k)
(N点)
(3)将X(k)的虚部Im[X(k)]改变符号,求得其共轭X*(k)
(4)计算Rxy(k)=X*(k)Y(k)
(5)用IFFT求得相关序列Rxy(n)
Rxy(n)=IFFT[Rxy(k)](N点)
如果x(n)=y(n),则求得的是自相关序列Rxx(n)
4.2.3Chirp-Z变换
采用FFT算法可以很快的计算出全部DFT值,也即Z变换在单位圆上的全部等间隔采样值。
但是,很多场合下,并非整个单位圆上的频谱都是很有意义的,例如对于窄带信号过程,往往只需要对信号所在的一段频带进行分析,这是,希望采样能密集在这段频带内,而对频带以外的部分,则可以完全不管。
另外,有时也希望采样能不局限于单位圆上,例如,语音信号处理中,往往需要知道其Z变换的极点所在频率,如果极点位置离单位圆较远,,则其单位圆上的频谱就很平滑,如图所示,这是很难从中识别出极点所在的频率。
要是采样不是沿单位圆而是沿一条接近这些极点的弧线进行,则所得的结果将会在极点所在频率上出现明显的尖峰,从而可以较准确的测定极点频率。
螺线采样就是一种适用于这种需要的变换,并且可以采用FFT来快
速计算,这种变换也称为Chirp-Z变换,它是沿Z平面上的一段螺线作等分角的采样,这些采样点可表达为
k
z=AW−k,k=0,1,L,M−1
其中M为采样点的总数,A为起始点位置,这个位置可以进一步用它的半径
0
A0及相角θ0来表示
参数W可表示为
A=Aejθ0
0
W=WejΦ0
其中W0为螺线的伸展率,W0>1时螺线内缩(反时针方向);W0<1时螺
线外伸。
Φ0为螺线上采样点之间的等分角。
螺线采样点在Z平面上的分布可表示为下图。
下面分析这些点上采样值计算的特点。
假定x(n)是长度为N的有限
长信号序列,则其Z变换在采样点上的值为
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 详解 FFT 快速 傅里叶变换