05-第五章-快速傅里叶变换(蝶形运算)

第五章第五章快速傅里叶变换快速傅里叶变换2本章目录本章目录直接计算DFT的问题及改进的途径按按时间抽取时间抽取的的基基2-FFT2-FFT算法算法按按频率抽取频率抽取的基的基2-FFT2-FFT算算法法快速快速傅里叶逆变换傅里叶逆变换(IFFT)(IFFT)算法算法MatlabMatlab实现实现35.15.1引言引言DFTDFT在实际应用中很重要在实际应用中很重要::可以计算信号的可以计算信号的频谱、功率谱和线性卷积等。频谱、功率谱和线性卷积等。直接按直接按DFTDFT变换进行计算，当序列长度变换进行计算，当序列长度NN很很大时，计算量非常大，所需时间会很长。大时，计算量非常大，所需时间会很长。FFTFFT并不是一种与并不是一种与DFTDFT不同的变换，而是不同的变换，而是DDFTFT的一种快速计算的算法。的一种快速计算的算法。45.25.2直接计算直接计算DFTDFT的问题及改进的途径的问题及改进的途径DFTDFT的运算量的运算量设复序列x(n)长度为N点，其DFT为10()()NnkNnXkxnWk=0，，…，N-1（1）计算一个X(k)值的运算量复数乘法次数：N复数加法次数：N－155.2.1DFT5.2.1DFT的运算量的运算量（2）计算全部N个X(k)值的运算量复数乘法次数：N2复数加法次数：N(N－1)（3）对应的实数运算量1100()()[Re()Im()][ReIm]NNnknknkNNNnnXkxnWxnjxnWjW10{[Re()ReIm()Im]NnknkNNnxnWxnW[Re()ImIm()Re]}nknkNNjxnWxnW6一次复数乘法：4次实数乘法2次实数加法＋一个X(k)：4N次实数乘法＋2N+2(N-1)=2(2N-1)次实数加法所以整个N点DFT运算共需要：N×2(2N-1)=2N(2N-1)实数乘法次数：4N2实数加法次数：7DFTDFT运算量的结论运算量的结论N点DFT的复数乘法次数举例NN2NN22464404941612816384864256655361625651226214432102810241048576结论：当N很大时，其运算量很大，对实时性很强的信号处理来说，要求计算速度快，因此需要改进DFT的计算方法，以大大减少运算次数。85.2.25.2.2减少运算工作量的途径减少运算工作量的途径nkNW主要原理是利用系数的以下特性对DFT进行分解：nkNW（1）对称性()nkNW()kNnNW（2）周期性()()nNknkNnkNNNWWW（3）可约性mnknkmNNWW//nknkmNNmWW另外，12/NNWkNNkNWW)2/(95.35.3按时间抽取的基按时间抽取的基2-FFT2-FFT算法算法算法原理算法原理按时间抽取基按时间抽取基-2FFT-2FFT算法与直接计算算法与直接计算DDFTFT运算量的比较运算量的比较按时间抽取的按时间抽取的FFTFFT算法的特点算法的特点按时间抽取按时间抽取FFTFFT算法的其它形式流程算法的其它形式流程图图105.3.15.3.1算法原理算法原理1(2)()xrxr设N＝2L，将x(n)按n的奇偶分为两组：2(21)()xrxrr=0，1，…，12N10()[()]()NnkNnXkDFTxnxnW则1010)()(NnnnkNNnnnkNWnxWnx为奇数为偶数11120)12(1202)12()2(NrkrNNrrkNWrxWrx1202212021)()(NrrkNkNNrrkNWrxWWrx)()(21kXWkXkN1010)()(NnnnkNNnnnkNWnxWnx为奇数为偶数式中，X1(k)和X2(k)分别是x1(n)和x2(n)的N/2的DFT。另外，式中k的取值范围是：0，1，…，N/2－1。12因此，只能计算出X(k)的前一半值。12()()()kNXkXkWXk后一半X(k)值，N/2，N/2＋1，…，N？rkNW2(2)2rNkNW利用可得到1()2NXk21(2)120()NrNkNrxrW21120()NrkNrxrW)(1kX同理可得22()()2NXkXk13考虑到kNkNNNkNNWWWW2)2(因此可得后半部分X(k))2()2()2(221NkXWNkXNkXNkN12()()()kNXkXkWXk及前半部分X(k))()(21kXWkXkNk=0，1，…，N/2－1k=0，1，…，N/2－114蝶形运算蝶形运算12()()()kNXkXkWXk12()()()kNXkXkWXk蝶形运算式蝶形运算信号流图符号因此，只要求出2个N/2点的DFT，即X1(k)和X2(k)，再经过蝶形运算就可求出全部X(k)的值，运算量大大减少。15以以88点为例第一次按奇偶分解点为例第一次按奇偶分解0NW1NW2NW3NW以N=8为例，分解为2个4点的DFT，然后做8/2=4次蝶形运算即可求出所有8点X(k)的值。16蝶形运算量比较蝶形运算量比较复数乘法次数：N2复数加法次数：N(N－1)复数乘法次数：2(N/2)2+N/2=N2/2+N/2复数加法次数：2(N/2)(N/2－1)+2N/2=N2/2N点DFTDFT的运算量的运算量分解一次后所需的运算量＝分解一次后所需的运算量＝22个个N/2N/2的的DFTDFT＋＋N/N/22蝶形：蝶形：因此通过一次分解后，运算工作量减少了差因此通过一次分解后，运算工作量减少了差不多一半。不多一半。17进一步按奇偶分解进一步按奇偶分解由于N＝2L，因而N/2仍是偶数，可以进一步把每个N/2点序列再按其奇偶部分分解为两个N/4点的子序列。以N/2点序列x1(r)为例1314(2)()0,1,,1(21)()4xlxlNlxlxl则有rkNNrWrxkX212011)()(klNNllkNNlWlxWlx)12(21401221401)12()2(lkNNlkNlkNNlWlxWWlx41404241403)()()()(42/3kXWkXkNk=0,1,…,14N18且13/24()()4kNNXkXkWXkk=0,1,…,14N由此可见，一个N/2点DFT可分解成两个N/4点DFT。同理，也可对x2(n)进行同样的分解，求出X2(k)。19以以88点为例第二次按奇偶分解点为例第二次按奇偶分解20算法原理算法原理13/40()lkNlxlW02(0)(4)xWx0(0)(4)NxWx对此例N=8，最后剩下的是4个N/4=2点的DFT，2点DFT也可以由蝶形运算来完成。以X3(k)为例。/4133/40()()NlkNlXkxlWk=0,1即03323(0)(0)(1)XxWx13323(1)(0)(1)XxWx12(0)(4)xWx0(0)(4)NxWx这说明，N=2M的DFT可全部由蝶形运算来完成。21以以88点为例第三次按奇偶分解点为例第三次按奇偶分解N=8按时间抽取法FFT信号流图225.3.25.3.2按时间抽取基按时间抽取基2-FFT2-FFT算法与直接计算算法与直接计算DFTDFT运算量的比较运算量的比较按时间抽取法FFT的信号流图可知，当N=2L时，共有级形运算；每级都由个蝶形运算组成，而每个蝶形有次复乘、次复加，因此每级运算都需次复乘次复加。LN/2N/212N23这样级运算总共需要：L复数乘法：NNLN2log22复数加法：NNLN2log直接DFT算法运算量复数乘法：复数加法：N2N(N－1)直接计算DFT与FFT算法的计算量之比为MNNNNNM222log2log224FFTFFT算法与直接算法与直接DFTDFT算法运算量的比较算法运算量的比较NN2计算量之比MNN2计算量之比M2414.01281638444836.641644.025665536102464.0864125.45122621442304113.816256328.0102410485765120204.83210288012.82048419430411264372.464404919221.4NN2log2NN2log2255.3.35.3.3按时间抽取的按时间抽取的FFTFFT算法的特点算法的特点序列的逆序排列同址运算（原位运算）蝶形运算两节点间的距离的确定rNW26序列的逆序排列序列的逆序排列)(01221)()(BINMMDECnnnnnn由于x(n)被反复地按奇、偶分组，所以流图输入端的排列不再是顺序的，但仍有规律可循：因为N=2M，对于任意n（0≤n≤N-1)，可以用M个二进制码表示为：10,,,,,01221nnnnnMMn反复按奇、偶分解时，即按二进制码的“0”“1”分解。序列的逆序排列27倒位序的树状图（倒位序的树状图（N=8N=8））28码位的倒位序码位的倒位序(N=8)(N=8)自然顺序n二进制数倒位序二进制数倒位序顺序数0000000010011004201001023011110641000011510110156110011371111117ˆn29倒位序的变址处理（倒位序的变址处理（N=8N=8））30同址运算同址运算（原位运算）（原位运算）某一列任何两个节点k和j的节点变量进行蝶形运算后，得到结果为下一列k、j两节点的节点变量，而和其他节点变量无关。这种原位运算结构可以节省存储单元，降低设备成本。)(kX)2(NkX)(1kX)(2kXkNW运算前运算后)2(NkA)(kA)2(NkA)(kA例同址运算（原位运算）31观察原位运算规律观察原位运算规律32蝶形运算两节点间的距离蝶形运算两节点间的距离以N=8为例：第一级蝶形，距离为：第二级蝶形，距离为：第三级蝶形，距离为：规律：对于共L级的蝶形而言，其m级蝶形运算的节点间的距离为12412m蝶形运算两节点间的距离33的确定的确定rNW以N=8为例：0,10224/jWWWWmjjNrNm时，1,0,2422/jWWWWmjjjNrNm时，3,2,1,0,382jWWWWmjjjNrNm时，级：第LNM,212,,2,1,0,12LjrNjWWLMLMLMLN2222LMLMMLMLjNjNjjNjjNrNWeeWW222222rNW的确定345.45.4按频率抽取的基按频率抽取的基2-FFT2-FFT算法算法算法原理再把输出X(k)按k的奇偶分组先把输入按n的顺序分成前后两半设序列长度为N=2L，L为整数前半子序列x(n)后半子序列)2(Nnx0≤n≤12N0≤n≤12N355.4.15.4.1算法原理算法原理10()()NnkNnXkxnW由DFT定义得/2110/2()()NNnknkNNnnNxnWxnW12/0)2(12/0)2()(NnkNnNNnnkNWNnxWnx12/02)2()(NnnkNkNNWWNnxnxk=0，1，…，N36由于1222jNNjNNeeW/2120()()()2NNknkNNnNXkxnxnWW所以kkNNW)1(2则12/0)2()1()()(NnnkNkWNnxnxkXk=0，1，…，N37然后按k的奇偶可将X(k)分为两部分221krkrr=0，1，…，12N则式12/0)2()1()()(NnnkNkWNnxnxkX可转化为nrNNnWNnxnxrX212/0)2()()2(12/02/)2()(NnnrNWNnxnx)12(12/0)]2()([)12(rnNNnWNnxnxrXnrNnNNnWWNnxnx212/0})]2()({[38/21/20(2)()()2NnrNnNXrxnxnW令nNWNnxnxnxNnxnxnx)2()()()2()()(21n=0，1，…，12N代入/2120(21){[()()]}2NnnrNNnNXrxnxnWWnrNNnnrNNnWnxrWnxr2120221201)()12()()2(r=0，1，…，12N可得为2个N/2点的DFT，合起来正好是N点X(k)的值。39蝶形运算蝶形运算nNWNnxnxnxNnxnxnx)2()()()2()()(21将称为蝶形运算与时间抽选基2FFT算法中的蝶形运算符号略有不同。40例按频率抽取例按频率抽取(N=8)(N=8)例按频率抽取，将N点DFT分解为两个N/2点DFT的组合(N=8)41与时间抽取法的推导过程一样，由于N=2L，N/2仍然是一个偶数，因而可以将每个N/2点DFT的输出再分解为偶数组与奇数组，这就将N/2点DFT进一步分解为两个N/4点DFT。N=8425.4.25.4.2频率抽取法与时间抽取法的异频率抽取法与时间抽取法的异同同频率抽取法输入是自然顺序，输出是倒位序的；时间抽取法正好相反。频率抽取法的基本蝶形与时间抽取法的基本蝶形有所不同。频率抽取法运算量与时间抽取法相同。频率抽取法与时间抽取法的基本蝶形是互为转置的。435.55.5快速傅里叶逆变换快速傅里叶逆变换(IFFT)(IFFT)算算法法MN2IDFT公式10)(1NknkNWkXNkXIDFTnxDFT公式nkNNnWnxnxDFTkX10)()()(比较可以看出，nkNWnkNWMN211IDFT多出M个1/2可分解到M级蝶形运算中。44例频率抽取例频率抽取IFFTIFFT流图流图(N=8)(N=8)45快速傅里叶逆变换另一种算法快速傅里叶逆变换另一种算法10)(1)]([NknkNWkXNkXIDFT10)(1NknkNWkXN10)(1NknkNWkXN1[()][()]IFFTXkFFTXkN()Xk求共轭()XkFFT求[()]FFTXk[()]FFTXkN除以()xn求共轭465.8Matlab5.8Matlab实现实现用FFT进行谱分析的Matlab实现用CZT进行谱分析的Matlab实现在Matlab中使用的线性调频z变换函数为czt，其调用格式为>>X=czt(x,M,W,A)其中，x是待变换的时域信号x(n)，其长度为N，M是变换的长度，W确定变换的步长，A确定变换的起点。若M=N，A=1，则CZT变成DFT。475.8.15.8.1用用FFTFFT进行谱分析的进行谱分析的MatlabMatlab实实现现例5.1设模拟信号，以t=0.01n(n=0:N-1)进行取样，试用fft函数对其做频谱分析。N分别为：(1)N=45；(2)N=50；(3)N=55；(2)N=60。()2sin(4)5cos(8)xttt程序清单如下%计算N=45的FFT并绘出其幅频曲线N=45;n=0:N-1;t=0.01n;q=n2pi/N;x=2sin(4pit)+5cos(8pit);y=fft(x,N);figure(1)subplot(2,2,1)plot(q,abs(y))title('FFTN=45')48例例5.15.1程序清单程序清单%计算N=50的FFT并绘出其幅频曲线N=50;n=0:N-1;t=0.01n;q=n2pi/N;x=2sin(4pit)+5cos(8pit);y=fft(x,N);figure(1)subplot(2,2,2)plot(q,abs(y))title('FFTN=50')49%计算N=55的FFT并绘出其幅频曲线N=55;n=0:N-1;t=0.01n;q=n2pi/N;x=2sin(4pit)+5cos(8pit);y=fft(x,N);figure(1)subplot(2,2,3)plot(q,abs(y))title('FFTN=55')50%计算N=60的FFT并绘出其幅频曲线N=60;n=0:N-1;t=0.01n;q=n2pi/N;x=2sin(4pit)+5cos(8pit);y=fft(x,N);figure(1)subplot(2,2,4)plot(q,abs(y))title('FFTN=60')51例例5.15.1程序运行结果程序运行结果02468050100150FFTN=4502468050100150FFTN=5002468050100150FFTN=5502468050100150FFTN=60从图中可以看出，这几种情况下均有较好的精度。52例例5.15.1程序运行结果分析程序运行结果分析分析：由t=0.01n进行取样可得，采样频率fs=100Hz。而连续信号的最高模拟角频率为Ω＝8π，由Ω＝2πf可得，最高频率为8π/2π=4Hz。因此，满足采样定理的要求。采样序列为()2cos(4)5cos(8)xnTnTn48()2cos5cos100100xnnn即为周期序列，周期N=50。将程序中plot改为stem函数，则可以更清楚地看出频谱。53例例5.15.1修改程序运行结果修改程序运行结果02468050100150FFTN=4502468050100150FFTN=5002468050100150FFTN=5502468050100150FFTN=60