关于快速傅里叶变换FFT应用中的几个问题

2.4关于FFT应用中的几个问题1、IDFT的高效算法上述FFT算法流图也可以用于离散傅里叶逆变换(InverseDiscreteFourierTransform，简称IDFT)。比较DFT和IDFT的运算公式：1010()[()]()1()[()]()NkNnNknNkXkDFTxnxnWxnIDFTxnXkWN图DIT―IFFT运算流图WN0WN1WN2WN3WN0WN0N1x(0)x(4)x(2)x(6)x(4)x(5)x(3)x(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)WN2WN2N1N1N1N1N1N1N1图DIT―IFFT运算流图(防止溢出)WN02121x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)212121WN121WN221WN3212121WN021WN2212121WN021WN22121WN02121WN0212121WN021WN021如果希望直接调用FFT子程序计算IFFT，则可用下面的方法：由于10101()()1()()NknNkNknNkxnXkWNxnXkWN对上式两边同时取共轭，得1011()[()][()]NknNkxnXkWDFTXkNN2、实数序列的FFT以上讨论的FFT算法都是复数运算,包括序列x(n)也认为是复数,但大多数场合,信号是实数序列,任何实数都可看成虚部为零的复数,例如,求某实信号x(n)的复谱,可认为是将实信号加上数值为零的虚部变成复信号(x(n)+j0),再用FFT求其离散傅里叶变换。这种作法很不经济,因为把实序列变成复序列,存储器要增加一倍,且计算机运行时,即使虚部为零,也要进行涉及虚部的运算,浪费了运算量。合理的解决方法是利用复数据FFT对实数据进行有效计算,下面介绍两种方法。（1）用一个N点FFT同时计算两个N点实序列的DFT设x(n)、y(n)是彼此独立的两个N点实序列,且X(k)=DFT[x(n)]，Y(k)=DFT[y(n)]则X(k)、Y(k)可通过一次FFT运算同时获得。首先将x(n)、y(n)分别当作一复序列的实部及虚部,令g(n)=x(n)+jy(n)通过FFT运算可获得g(n)的DFT值)利用离散傅里叶变换的共轭对称性通过g(n)的FFT运算结果G(k),由上式可得到X(k)的值。kjGkGkjYkjXkYkXkYkjYkjXkXkjYkXkGirriiririrkNGkGjkNGkGkNGkGkXngiirr212121ReDFT*通过g(n)的FFT运算结果G(k),由上式也可得到X(k)的值。kNGkGjkNGkGkNGkGkjYngjiirr212121ImDFT*kNGkGjkNGkGkYrrii2121作一次Ｎ点复序列的FFT，再通过加、减法运算就可以将X(k)与Y(k)分离出来。显然，这将使运算效率提高一倍。（2）用一个N点的FFT运算获得一个2N点实序列的DFT设x(n)是2N点的实序列,现人为地将x(n)分为偶数组x1(n)和奇数组x2(n)x1(n)=x(2n)n=0,1,…,N-1x2(n)=x(2n+1)n=0,1,…,N-1然后将x1(n)及x2(n)组成一个复序列：y(n)=x1(n)+jx2(n)通过N点FFT运算可得到：Y(k)=X1(k)+jX2(k)，N点根据前面的讨论，得到)]()([2)()]()([21)(*2*1kNYkYjkXkNYkYkX为求2N点x（n）所对应X（k），需求出X（k）与X1（k）、X2（k）的关系：而10)12(21022120)12()2()()(NnknNNnnkNNnnkWnxWnxWnxkX10210)12()2(NnnkNkNNnnkNWnxWWnx1011022101011)12()()()2()()(NnnkNNnnkNNnNnnkNnkNWnxWnxkXWnxWnxkX所以1）由x1(n)及x2(n)组成复序列，经FFT运算求得Y(k)，2）利用共轭对称性求出X1(k)、X2(k)，3）最后利用上式求出X（k）,达到用一个N点的FFT计算一个2N点的实序列的DFT的目的。X(k)=X1(k)+W2NkX2(k)3、线性卷积的FFT算法线性卷积是求离散系统响应的主要方法之一,许多重要应用都建立在这一理论基础上,如卷积滤波等。以前曾讨论了用圆周卷积计算线性卷积的方法归纳如下:将长为N2的序列x(n)延长到L,补L-N2个零，将长为N1的序列h(n)延长到L,补L-N1个零，如果L≥N1+N2-1,则圆周卷积与线性卷积相等,此时,可用FFT计算线性卷积，方法如下:a.计算X(k)=FFT[x(n)]b.求H(k)=FFT[h(n)]c.求Y(k)=H(k)X(k)k=0～L-1d.求y(n)=IFFT[Y(k)]n=0～L-1可见，只要进行二次FFT,一次IFFT就可完成线性卷积计算。计算表明,L32时,上述计算线性卷积的方法比直接计算线卷积有明显的优越性,因此,也称上述循环卷积方法为快速卷积法。上述结论适用于x(n)、h(n)两序列长度比较接近或相等的情况,如果x(n)、h(n)长度相差较多,例如,h(n)为某滤波器的单位脉冲响应,长度有限,用来处理一个很长的输入信号x(n),或者处理一个连续不断的信号，按上述方法,h(n)要补许多零再进行计算,计算量有很大的浪费，或者根本不能实现。为了保持快速卷积法的优越性,可将x(n)分为许多段,每段的长度与h(n)接近,处理方法有两种：(1)重叠相加法——由分段卷积的各段相加构成总的卷积输出h(n)x(n)假定xi(n)表示x（n）序列的第i段：则输入序列可表为：于是输出可分解为：其中01)1()()(22NiniNnxnxiiinxnx)()(iiiinynhnxnhnxny)()(*)()(*)()()(*)()(nhnxnyii1）先对h（n）及xi（n）补零，补到具有N点长度，N=N1+N2-1。一般选N=2M。2）用基2FFT计算yi（n）=xi（n）*h（n）。3）重叠部分相加构成最后的输出序列。由于yi（n）的长度为N，而xi（n）的长度为N2，因此相邻两段yi（n）序列必然有N-N2=N1-1点发生重叠。iinyny)()(计算步骤：a.事先准备好滤波器参数H（k）=DFT[h(n)]，N点b.用N点FFT计算Xi（k）=DFT[xi(n)]c.Yi(k)=Xi(k)H(k)d.用N点IFFT求yi（n）=IDFT[Yi(k)]e.将重叠部分相加(2)重叠保留这种方法和第一种方法稍有不同，即将上面分段序列中补零的部分不是补零，而是保留原来的输入序列值，这时，如利用DFT实现h(n)和xi(n)的循环卷积，则每段卷积结果中有N1-1个点不等于线性卷积值需舍去。重叠保留法与重叠相加法的计算量差不多，但省去了重叠相加法最后的相加运算。4、4、用FFT计算相关函数相关的概念很重要，互相关运算广泛应用于信号分析与统计分析，如通过相关函数峰值的检测测量两个信号的时延差等。两个长为N的实离散时间序列x（n）与y（n）的互相关函数定义为1010)()()()()(NnNnxynynxnynxr则可以证明，rxy(τ)的离散傅里叶变换为Rxy（k）=X*(k)Y(k)其中X(k)=DFT[x(n)],Y(k)=DFT[y(n)],Rxy(k)=DFT[rxy(τ)]，0≤k≤N-1互相关函数定义为x(n)及y(n)的卷积公式相比较，我们可以得到相关和卷积的时域关系：1010)()()()()(NnNnxymnynxnymnxmr)()()()()(10mymxnynmxmfNn)()()()]([)()()(1010mymxnynmxnymnxmrNnNnxy10102)()(1)()()(NnNkkNjxyekYkXNnynxr因得证毕。)(*)]())((DFT[kXnRnxNN当x（n）=y（n）时，得到x（n）的自相关函数为：维纳——辛钦定理：自相关函数与信号功率谱互为傅立叶变换对。10)()()(NnxxnxnxrkNjNkekXN2102|)(|1上面的推导表明，互相关和自相关函数的计算可利用FFT实现。由于离散傅里叶变换隐含着周期性，所以用FFT计算离散相关函数也是对周期序列而言的。直接做N点FFT相当于对两个N点序列x(n)、y(n)作周期延拓，作相关后再取主值（类似圆周卷积）。实际一般要求的是两个有限长序列的线性相关，为避免混淆，需采用与循环卷积求线性卷积相类似的方法，先将序列延长补0后再用上述方法。利用FFT求两个有限长序列的线性相关：（1）设x(n)和y(n)的长均为N，求线性相关；（2）为了使两个有限长序列的线性相关可用其圆周相关代替而不产生混淆，选择周期L≥2N-1，且L=2m，以使用FFT，将x(n)，y(n)补零至长为L。（3）用FFT计算X(k)，Y(k)（k=0，1…，L-1）（4）R（k）=X*（k）y（k）（5）对R（k）作IFFT，取后N-1项，得前N项，得11)(mNmrxy10)(Nmmrxyx=[13-112331];y=[21-1120-13];k=length(x);xk=fft(x,2*k);yk=fft(y,2*k);rm=real(ifft(conj(xk).*yk));rm=[rm(k+2:2*k)rm(1:k)];m=(-k+1):(k-1);stem(m,rm)xlabel('m');ylabel('幅度');-8-6-4-202468-6-4-2024681012m幅度两个序列的自相关函数例8例9-20-1001020-50050100150200250300m-20-1001020-50050100150200250300m幅度幅度(a)(b)延迟序列的互相关函数(a)和自相关函数(b)5、用FFT计算二维离散的傅里叶变换二维信号有图象信号、时空信号、时频信号等。二维离散傅里叶变换可用于处理二维离散信号。二维离散傅里叶变换的定义为：101022),(),(MmNnlmMjknNjeemnxlkX1010,MmNnlmMknNwwmnxMjMNjNewew22,二维离散傅里叶变换可通过两次一维离散傅里叶变换来实现:1）作一维N点DFT（对每个m做一次，共M次）k=0,1,…,N-1，m=0,1,…,M-12）作M点的DFT（对每个k做一次，共N次）k=0,1,…,N-1，l=0,1,…,M-1102),(),(NnknNjemnxmkA102),(),(MmlmMjemkAlkXNMNMMNNMNMAAAAxxxx11111111这两次离散傅里叶变换都可以用快速算法求得,若M和N都是2的幂,则可使用基二FFT算法,所需要乘法次数为而直接计算二维离散傅里叶变换所需的乘法次数为（M+N）MN，当M和N比较大时用用FFT运算，可节约很多运算量。)(log2log2log2222MNMNNNMMMN