1、1离散傅里叶变换Discrete Fourier Transform2内容提要离散傅里叶变换(Discrete Fourier Transform,DFT)是时间函数是离散的,而且频谱函数也是离散的变换。离散傅里叶变换定义DFT物理意义DFT基本性质讨论频率取样理论。DFT的应用 3傅里叶变换的各种形式连续时间、离散频率的傅里叶变换对于周期为T的连续时间信号,可以采用傅里叶级数展开:连续时间、连续频率的傅里叶变换对于非周期的连续时间信号,可以进行傅里叶变换:它在时域和频域都是连续的。4离散时间、连续频率的傅里叶变换对于非周期的序列,其傅里叶变换在频域是以2为周期的连续函数。5设x(n)是一个长
2、度为M的有限长序列,则定义x(n)的N点离散傅里叶变换为10()()(),k=0,1,&,N-1 (3.1.1)NknNnX kDFT x nx n WX(k)的离散傅里叶逆变换为101()()(),k=0,1,&,N-1 (3.1.2)NknNnX kDFT x nX n WN式中,N称为DFT变换区间长度,NM,通常称(3.1.1)式和(3.1.2)式为离散傅里叶变换对。Note:有限长序列x(n)的DFT即X(k)仍是有限长序列。2jNNWe6例 3.1.1 x(n)=R4(n),求x(n)的8点和16点DFT.解:设变换区间N=8,则273880038()()sin()2,0,1,7s
3、in()8jknknnNj kX kx n Wekekk7对长度为M的序列x(n),其Z变换N点DFT进行对比,可以看出式中,表示z平面单位圆上辐角(k=0,1,N-1)的N个等间隔点。2,0,1,1jkNXkXekN8说明:序列x(n)的N点DFT是其Z变换在单位圆上的N点等角距取样,如图3.4(a)。序列x(n)的DFT是其FT在区间0,2上的N点等间隔取样。如图3.4(b)。9DFT变换对中(),kk mNNNWWk m N均为整数 所以式(3.1.1)中,X(k)满足1()010()()()()Nk mN nNnNknNnX kmNx n Wx n WX k同理可证明式(3.1.2)中
4、 x(n+mN)=x(n)10 x n任何周期为N的周期序列 都可以看做长度为N的有限长序列x(n)的周期延拓序列,而x(n)是 的一个周期 mNx nx nmNx nx n Rn x n(3.1.5)(3.1.6)11定义:为叙述方便,将式(3.1.5)该写成 表示x(n)以N为周期的周期延拓序列,符号(n)N表示n对模N的余数,即 这里k是商。的主值区间:周期序列 中从n=0到N-1的范围 的主值序列:主值区间上的序列 Nxn12由此对长度为N的序列x(n),且 ,则 的DFS为结论:与DFT定义比较,可见有限长序列x(n)的DFT即X(k)是x(n)的周期延拓序列 的离散傅里叶级数系数
5、的主值序列。例如,N=7,=x(n)7,则有 77770881xxxxxx Nx nxn 111000110011NNNknknknNNNNnnnNNknknNNnnX kx n WxnWx n Wx nX k WX k WNN NxnX K NX kX k Rk13解:因此得 X(0)=4.16114 X(1)=0.71063-j0.92558 X(2)=0.50746-j0.40597 X(3)=0.47017-j0.16987 X(4)=0.46235X(5)=0.47017+j0.16987X(6)=0.50746+j0.40597X(7)=0.71063+j0.92558Matlab
6、实现 fft1.m例3.1 求有限长序列的DFT,其中a=0.8,N=8。14关于离散傅里叶变换(DFT):序列序列x(n)在时域是有限长的在时域是有限长的(长度为长度为N),它的离散傅里叶变,它的离散傅里叶变换换X(k)也是离散、有限长的也是离散、有限长的(长度也为长度也为N)。n为时域变量,为时域变量,k为频域变量。为频域变量。DFT的物理意义:序列的物理意义:序列x(n)的的Z变换在单位圆上的等角距取变换在单位圆上的等角距取样。序列傅里叶变换在区间样。序列傅里叶变换在区间0,2上的等间隔取样。上的等间隔取样。离散傅里叶变换离散傅里叶变换(DFT)具有唯一性。具有唯一性。离散傅里叶变换与离
7、散傅里叶级数没有本质区别,离散傅里叶变换与离散傅里叶级数没有本质区别,DFT实实际上是离散傅里叶级数的主值,际上是离散傅里叶级数的主值,DFT也隐含有周期性。也隐含有周期性。15 DFT隐含着周期性,因此在讨论DFT的性质时,常与DFS的概念联系起来,并把有限长序列看作周期序列的一个周期来处理。设x1(n)和x2(n)的长度都为N,且它们对应的DFT分别为X1(k)和X2(k)。设x3(n)=ax1(n)+bx2(n),a和b都为常数,则 若它们长度不等,取长度最大者,将短的序列通过补零加长,注意此时DFT与未补零的DFT不相等。3.2 离散傅里叶变换的性质16 一个长度为N的序列x(n)的循
8、环移位定义为 循环移位分3步计算:(1)将x(n)延拓成周期为N的周期序列 ;(2)将 移位得 或 x(n+m)N;(3)对 x(n+m)N 取主值得 x(n+m)NRN(n)。这个过程如下图所示。a)序列的循环移位:17 从图中两虚线之间的主值序列的移位情况可以看出,当主值序列左移m个样本时,从右边会同时移进m个样本,而且好像是刚向左边移出的那些样本又从右边循环移了进来。因此取名“循环移位”。显然,循环移位不同于线性移位 181920对长度为N的有限长序列x(n),其循环移位后序列y(n)的DFT为 证明:b)时域循环移位定理时域循环移位定理 1100NNknknNNNNNnnY kDFT
9、y nxnmRn WxnmW令n+m=n,则有:11NmNmk nmkmknNNNNNnmnmY kDFT y nxnWWxnW knNNxnW 因为 以N为周期,上式中的求和区间任取一个周期即可,取主值区间为求和区间,得证。21若则c)频域循环移位定理频域循环移位定理 22 长度分别为N1和N2的有限长序列x1(n)和x2(n)的N点DFT分别为:(N=max N1,N2)。则由上式表示的卷积称为循环卷积X1(k)=DFTx1(n)X2(k)=DFTx2(n)如果X(k)=X1(k)X2(k)11201210IDFTNNNmNNNmx nX kx m xnmRnxm xnmRn23循环卷积的
10、过程:(1)周期延拓 x2(m)x2(m)N(2)折叠 x2(m)Nx2(-m)N(3)移位和取主值 x2(-m)Nx2(n-m)NRN(m)(4)相乘 x2(n-m)NRN(m)x1(m)x2(n-m)NRN(m)(5)相加 summ0,1,N-1循环反转序列Note:两个长度为N的序列的循环卷积长度仍为N,(与线性卷积不同),记为:112120121210NNNmNNNmx nx nxnx m xnmRnxnx nxm xnmRn242526循环卷积计算说明:x1(n)的N个值按顺时针方向均匀分布在内圆周上,x2(n)的N个值按反时针方向均匀分布在外圆周上,把内外圆周上对应的数值两两相乘,
11、然后把乘积相加就得到y(0)。若将外圆周顺时针方向转动一格,将内外圆周上对应的数值两两相乘并把乘积相加,便得到y(1)。依次类推,可以得出y(n)的其它值。因此循环卷积也叫做圆卷积。考虑到DFT关系的对偶性,可以证明,长为N的两序列之积的DFT等于它们的DFT的循环卷积除以N,即 频域循环卷积定理2728 是长度为N的序列x(n)的复共轭序列,xn X kDFT x n则,01DFT xnXNkkN 且 0X NX类似,01DFT xNnXkkNNote:对实序列有 X kXNk293.2.5 DFT的共轭对称性的共轭对称性 1.有限长共轭对称序列和共轭反对称序列分别用xep(n)和xop(n
12、)表示有限长共轭对称序列和共轭反对称序列,则二者满足如下定义式:xep(n)=x*ep(N-n),0nN-1 (3.2.9)xop(n)=-x*op(N-n),0nN-1 (3.2.10)当N为偶数时,将上式中的n换成N/2-n可得到()(),01222()(),01222epepopopNNNxnxnnNNNxnxnn 30图 3.2.3 共轭对称与共轭反对称序列示意图(图中*表示对应点为序列取共轭后的值)31 如同任何实函数都可以分解成偶对称分量和奇对称分量一样,任何有限长序列x(n)都可以表示成其共轭对称分量和共轭反对称分量之和,即 x(n)=xep(n)+xop(n),0nN-1 (3
13、.2.11)将上式中的n换成N-n,并取复共轭,再将式(3.2.9)和式(3.2.10)代入得到 x*(N-n)=x*ep(N-n)+x*op(N-n)=xep(n)-xop(n)(3.2.12)xep(n)=1/2x(n)+x*(N-n)(3.2.13)xop(n)=1/2x(n)-x*(N-n)(3.2.14)32 2.DFT的共轭对称性(a)若将序列x(n)分成实部xr(n)与虚部xi(n),即 x(n)=xr(n)+jxi(n)根据复共轭序列的DFT可得 11DFTDFT2211DFTDFT()22repiopxnx nxnX kXNkXkjx nx nxnX kXNkXk再由DFT的
14、线性性质可得 X=DFTDFTDFTDFT()ririepopkx nxnjx nxnjx nXkXk33(b)若将序列x(n)分成共轭对称部分xep(n)与共轭反对称部分xop(n),即 x(n)=xep(n)+xop(n)根据复共轭序列的DFT可得 11DFTDFTRe2211DFTDFTIm()22epopxnx nxNnX kXkX kxnx nxNnX kXkjX k 因此 X=DFTDFTepopRIkx nxnxnXkjXk34结论:如果序列x(n)的DFT为X(k),则 x(n)的实部和虚部(包括j)的DFT分别为X(k)的共轭对称分量和共轭反对称分量;x(n)的共轭对称部分和
15、共轭反对称部分的DFT分别为X(k)的实部和虚部乘以j.35有限长实序列DFT的共轭对称性:对长度为N的实序列,X(k)=DFTx(n),则 X(k)共轭对称,即 若x(n)=x(N-n),则X(k)实偶对称,即 若x(n)=-x(N-n),则X(k)纯虚奇对称,即 X=kX Nk,01X kXNkkN X=kX Nk36 这意味着,对于时间有限信号,可以像频带有限信号进行时域采样而不丢失任何信息一样,可以在频域上进行采 样而不丢失任何信息。这正是傅里叶变换中时域和频域对偶关系的反映,这有着十分重要的意 义。DFT实现了频域离散化,开辟了在频域采用数字技术处理的新领域。这使我们自然想到,对于任
16、意一个频率特性,是否均能用频域采样的办法来逼近,这是一个很吸引人的问题,因为用频率采样来逼近,可使问题大大简化。因此我们要讨论频率采样的可行性以及所带来的误差。3.3 频率域采样 37频率取样是指对序列的傅里叶变换或系统的频率特性进行取样。本节讨论在什么条件下能够用得到的频谱取样值无失真地恢复原信号或系统。设任意长序列x(n)绝对可和,其Z变换表示为如果在单位圆上对X(z)进行等角距取样,取样点数为N,则得 22()()()jkNjknNz enX kX zx n e根据DFT的定义,对X(k)求反变换 =IDFTpxnX k38根据上面两式可得:因为所以 上式表明,在z平面的单位圆上对序列的
17、Z变换进行等角距取样,将导致时间序列的周期延拓。这一结果与对连续时间信号取样导致频谱周期延拓类似。现在我们来考察xp(n)与原序列x(n)的关系,看它如何才能代表原序列x(n)。110011=NNk n mkmknpNNNkmmkxnx m WWx mWNN 101,1,0,Nk n mNkmnrNWrmnrNN为任意整数=prxnx nrN39 xp(n)是原非周期信号x(n)的周期延拓序列,因此xp(n)是一个周期序列,其主值为 在x(n)为有限长度M的情况下,如果取样点NM,那么x(n)周期延拓的结果不会产生混叠。这时,xp(n)的主值xN(n)与原序列x(n)一样,因此xN(n)完全能
18、代表原序列x(n),可由频域采样X(k)恢复x(n)。如果N=M,则有42令则 上式就是用X(z)在单位圆上的N个取样值X(k)表示X(z)的内插公式,内插函数为 。因为1kNNW所以 X z 1111NkkNzzNWz 10NkkX zX kz kz43如果 则可以得到傅里叶变换的内插公式,结论:长度为N的序列x(n),其N个频域取样值X(k)可以不失真地代表它,X(k)还能完整的表示序列的Z变换X(z)和傅里叶变换 。2/10111j Nkjk NNjkkeNeX eX k443.4 DFT应用举例 1.线性卷积实际应用中一般以线性卷积和相关运算处理为依据,如一个FIR数字滤波器的输出等于
19、输入与滤波器的单位取样响应的线性卷积。DFT计算循环卷积DFT的快速算法FFT的出现,使DFT在数字通信、语言信号处理、图像处理、功率谱估计、仿真、系统分析、雷达理论、光学、医学、地震以及数值分析等各个领域都得到广泛应用。2.谱分析45线性卷积 线性卷积不受主值区间限制循环卷积 在一定条件下与线性卷积相等。两个长度都为N的因果序列的循环卷积仍是一个长度为N的序列,而它们的线性卷积却是一个长度为2N-1的序列。3.4.1 利用循环卷积计算线性卷积 46 如果能将线性卷积转化成循环卷积,那么根据DFT的循环卷积性质,就能够用循环卷积来计算线性卷积,而循环卷积可以用FFT 进行快速计算。因此,首先需
20、要讨论在什么条件下,循环卷积与线性 卷积相等的问题。设h(n)和x(n)分别是长度为N和M的有限长序列,它们的线性卷积和循环卷积分别为 1010NlmLcLLmynh nx nh m x nmynh nx nh m xnmRn其中,L=maxN,M47所以对照线性卷积的公式,可以看出因 为 Lqxnx nqL 1010NcLmqNLqmynh mx nmqL Rnh m x nqLm Rn 10Nlmh m x nqLmynqL48 yc(n)是yl(n)以L为周期的周期延拓序列的主值序列。而yl(n)是个长度为N+M-1的序列,所以(1)如果L=N+M-1)计算L点DFTh(n),DFTxk
21、(n)计算 Yk(n)=Xk(n)H(K)用L点IFFT求yk(n)=IDFTYk(k)将yk(n)的重叠部分相加,最后输出为由于yk(n)的长度为N+M-1,而xk(n)长度为M,所以相邻两段yk(n)序列必然有N-1个点发生重叠,如图3.4.4所示。这些重叠部分应该相加起来才能构成最后的输出序列,这就是“重叠相加法”这一名称的由来。54图 3.4.4 重叠相加法卷积示意图 M0NMMx1(n)x0(n)x2(n)N M 1N M 1y0(n)y1(n)N M 1y2(n)2MM3M N 10N 1y(n)y0(n)y1(n)y2(n)nnnnnnh(n)553.4.2 用DFT对信号进行谱
22、分析 所谓信号的谱分析就是计算信号的频谱,包括振幅谱、相位谱和功率谱。处理对象:连续信号和离散信号1 用DFT对连续信号进行谱分析(近似分析)设连续信号xa(t)(经预滤波和截取处理的有限长带限信号)持续时间为Tp,最高频率为fc,如图3.4.5所示。xa(t)的傅里叶变换为()()()j taaaXjFT x tx t edt 56 对xa(t)以采样间隔T1/2fc(即fs=1/T2fc)采样得x(n)=xa(nT)。x(n)的傅里叶变换为X(ej).12()jamX eXjmTTT X(ej)与 的关系为(3.4.6)X(n)的长度为ppsTNT fT()aXj将 代入T 121()de
23、fj TaamX eXjmXjTTT57X(n)的N点DFT 2()DFT01jNkNX kx nX ekN 22()22aafaafXfXjXjfXfXXf21212(),01jkNaapX kX eXkXkkNTNTTT代入式(3.4.6),可得所以,令58 11(),0,1,1pdefaakkfNTTaNX kXfXkFkNTTXkFTX kT DFT x n可得F表示对模拟信号的频谱的采样间隔,称为频率分辨率,即频域取样中两相邻点间的频率间隔。说明:由上式可得,对连续信号采样并进行DFT再乘以T,就可近似得到模拟信号频谱的周期延拓函数在第一个周期0,Fs上的N点等间隔采样。如图3.4.
24、5.11spFFTNTN5960结论:对满足假设的持续时间有限的带限信号,在满足时域采样定理时,包含了模拟信号频谱的全部信号,可由X(k)恢复 。缺点:栅栏效应对持续时间有限的带限信号aXkF aaXjxt 或61截断效应用DFT分析理想低通滤波器单位冲激响应的频率响应特性截取理想低通滤波器的单位冲激响应的一段Tp:sin()()ath tt 假设Tp=8 s,采样间隔T=0.25 s(即采样速度fs=4 Hz),采样点数N=Tp/T=32。此时频域采样间隔F=1/NT=0.125 Hz。则Ha(kF)=TDFTh(n),0k16其中 h(n)=ha(nT)R32(n)注意:对实信号,其频谱函
25、数具有共轭对称性,只需分析正频率频谱62图 3.4.7 用DFT计算理想低通滤波器频响曲线 63从上图(c)可看出:低频部分近似理想低通频响特性;高频误差较大 截断效应整个频响都有波动减少截断效应的途径:加长信号分析时间Tp,增加采样点数 先用窗函数处理64连续信号谱分析的参数选择原则:关心的问题:谱分析范围0,Fs/2和频率分辨率F1)为避免频谱混叠现象,Fs 2fc2)谱分辨率F=Fs/N,若N不变,要提高频谱分辨率,必须降低Fs 若Fs不变,为提高频谱分辨率,可增加采样点数N,即增加观察时间Tp。21cpfNFTF选取原则:65例3.4.2 对实信号进行谱分析,要求谱分辨率F10 Hz,
26、信号最高频率fc=2.5 kHz,试确定最小记录时间Tpmin,最大的采样间隔Tmax,最少的采样点数Nmin。如果fc不变,要求谱分辨率增加一倍,最少的采样点数和最小的记录时间是多少?解:因此Tpmin=0.1 s,因为要求fs2fc,所以110.110PTsF3maxmin110.2 1022250022250050010ccTsffNF66为使频率分辨率提高一倍,F=5 Hz,要求minmin225001000510.25pNTs672 用DFT对序列进行谱分析单位圆上的Z变换序列的傅里叶变换序列的DFT:0,2上对傅里叶变换的N点等间隔采样。频谱分辨率:2/N可用DFT计算序列的FT周
27、期为N的周期序列 的频谱函数:x n21022()()()()()()()jkNjknNnX eFT x nX kkNNX kDFS x nx n e 其中 Note:周期序列的频谱结构可用其离散傅里叶级数系数表示68由DFT的隐含周期性,截取 的主值序列 并进行N点DFT得到()x n ()Nx nx n Rn()()()()()()NNX kDFT x nDFT x n RnX k Rk()x n()x n因此可用X(k)表示 的频谱结构。如果截取长度M等于 的整数个周期,即M=mN,m为正整数,则2102(1)0()()()()()()()0,1,1MMMknMMMnm NknmNnxn
28、x n RnXkDFT xnx n ex n ekmN69令 n=n+i N,i=0,1,m-1,n=0,1,N-1,则2()110221100210210()()()()()nrN kmNjmNMrnnmNjkjrkmNmrnmjrkmrmjrkmrXkx nrN ex n eekXemkXem 210,0,Mjkrmrme因为 k/m=整数k/m整数 70(),()0,MkmXXkmk/m=整数k/m整数 结论:k=im时,表示周期序列的第i次谐波谱线,幅度扩大m倍,故截取周期序列的整数个周期进行DFT可得到其频谱结构。周期序列的周期不知道时的处理:截取M点做DFT 截取长度扩大1倍做DF
29、T 分析主谱频率差别71 在实际应用中,有时只对信号的某一频段感兴趣,或只需计算单位圆上某一段的频谱值。例如,在对窄带信号进行分析时,常希望在窄频带内对频率的取样很密集,以便提高频率分辨率,而在窄频带外不予以考虑。对于这种情况,如果采用DFT方法,则需要在窄频带内外都增加频域取样点,而窄频带外的计算量是浪费的。此外,有时对非单位圆上的取样感兴趣,例如在语音信号处理中,常常需要知道其Z变换的极点所在处的复频率,这时就需要在这些极点附近的曲线上进行频域取样,这样,就要沿着螺旋线对Z变换取样。这种沿螺旋线上取样点计算的Z变换,称为线性调频Z变换(Chirp Z Transform,简称CZT)。72
30、图 3.4.7 单位圆与非单位圆采样 73要计算序列在半径为r的圆上的频谱,那么N个等间隔采样点为 ,k=0,1,2,N-1,zk点的频谱分量为 2jkNkzre210()()()kNjknNkz znX zX zx n r e()()nx nx n r210()()(),01NjknNknX zx n eDFT x nkN令 若要计算有限角度2/M内的N点等间隔频谱采样,可取L=MN,作N点DFT,只取分析角度内的N点等间隔采样。对曲线,分别计算很多弧线上的采样,运算量大,效率低。743 Chirp-Z变换一个长度为N的序列x(n),其Z变换为为了使z可以沿着z平面更一般的路径(不只是单位圆
31、)取值,可以沿一段螺旋线对z作等分角取样,这些取样点上的zk表示为其中M为所要分析的复频谱的点数,不一定等于N。W和A为任意复数,极坐标下可表示为得到A0,W0正实数75取样点zk所在的路径如图所示(1)A0表示起始取样点z0的矢量长度,通常A01,否则将处于单位圆|z|=1之外。(2)0表示起始取样点z0的矢量的相角,它可以是正值或负值。(3)0表示两相邻取样点矢量之间的角度差。0为正时,表示zk的路径沿逆时针方向旋转;0为负时,zk的路径沿顺时针方向旋转。(4)W0表示螺旋线的伸展率。W01时,随着k的增加螺旋线向内盘旋;W01时,则随k的增加螺旋线向外盘旋;W0=1对应于半径为A0的一段
32、弧线,在A0=1时这段弧线是单位圆的一部分。76将zk代入Z变换公式得到1010()()(),01NknknNnknnX zx nAWx n A WkM2221()2nknkkn直接计算上式,总共要算M个取样点,需要NM次复数乘法和(N-1)M次复数加法,这与DFT的直接计算类似。当N和M很大时,这个计算量是很大的,因而限制了运算速度。如果将式中的因子Wnk的幂nk作如下的变换,则可以将以上运算转换成卷积的形式,并可采用FFT来计算,从而大大提高运算速度。具体来说,将nk写成以下形式:77222222221()/201/2/2()/20/2/2()()()()()()Nknnkk nnNknn
33、k nnnnnX zx n A WWx n A WWy nx n A Wh nW令则21/20()()(),01NkknX zWy n h knkM所以,长度为N的序列x(n)的M点Chirp-z变换可通过计算y(n)与h(n的线性卷积,然后再乘上Wk2/2 来得到78图 3.4.10 Chirp z变换计算框图 h(n)x(n)y(n)X(zk)v(n)V(k)22kW22nnWA79用DFT计算Chirp-Z变换的原理和实现(1)选择满足条件L(N+M-1)和L=2M的整数L值。(2)对 y(n)添补零取样值,构成长度为L点的序列y(n),用FFT计算y(n)的L点DFT。(3)按下式构造
34、长为L的序列hL(n),计算其L点DFT2/2(),01()0,1nnx n A WnNy nNnL22()/201()1()()01nLn LLWnMh nWMnLH kDFT h nkL802/2()()(),01()(),01kkV nIDFT Y k H knLX zWV kkM(4)计算()()Y kH k(5)(6)Chirp-Z 变换特点:它的输入序列长度N和输出序列长度M可以不相等,且均可为任意数,包括素数;各zk点间的角度间隔0可以是任意的,因而频率分辨率可以调整;计算Z变换的取样点的轨迹可以不是圆而是螺旋线;起始点z0可任意选定,也就是可以从任意频率开始对输入数据进行分析。814 用DFT进行谱分析的误差问题处理对象:连续信号和数字信号 连续信号采样和截断 非时限数据序列截断误差现象:混叠现象:采样定理限制 栅栏效应:N点的DFT是在频率区间0,2上FT的N点等间隔采样,采样点之间的频谱不可见。提高精度的方法:对有限长序列补零,对无限长序列增加观察时间,对连续信号,Fs足够高82截断效应 泄露:频谱展宽,使频谱变模糊,谱分辨率降低。谱间干扰:旁瓣引起不同频率分量间的干扰。例如,x(n)=cos(0n),0=/4其频谱为