§7-4-几种特殊的FFT算法

•一、离散傅立叶反变换（IDFT）的快速变换（IFFT）•二、N为组合数的FFT算法•三、实数序列的FFT一、离散傅立叶反变换（IDFT）的快速变换（IFFT）IDFTFFTNWWDFTWnxnxDFTkXWkXNkXIDFTnxnkNnkNNknkNNknkN算法都可以拿来运算或频率抽取抽取）那么以上讨论的时间（将运算结果都除以改成运算中的每个系数只要把3)2()1()()]([)()(1)]([)(1010IDFTFFTNWWDFTWnxnxDFTkXWkXNkXIDFTnxnkNnkNNknkNNknkN算法都可以拿来运算或频率抽取抽取）那么以上讨论的时间（将运算结果都除以改成运算中的每个系数只要把3)2()1()()]([)()(1)]([)(101022020/3/432020/3/441**01()()NnkNkxnXkWN*1***011()()[()]NnkNkxnXkWDFTXkNN101()()NnkNnxnXkWN共轭FFT共轭乘1/N()Xk*()Xk()xn直接调用FFT子程序计算IFFT的方法：2020/3/4二、N为组合数的FFT算法5N为组合数的FFT（任意基数的FFT算法）以上讨论的都是以2为基数的FFT算法，即N=2M，这种情况实际上使用得最多，其优点是程序简单，效率很高，使用起来非常方便，实际应用时，有限长序列的长度N很大程度上由人为因素确定，因此多数场合可取N=2M，从而直接使用以2为基数的FFT算法。2020/3/46如N不能人为确定，N的数值也不是以2为基数的整数次方，处理方法有两种:•①补零，将x(n)补零，使N=2M。例如N=30，补上x(30)=x(31)=0两点，使N=32=25，这样可直接采用以2为基数M=5的FFT程序。•②采用任意数基数的DFT算法如要求准确的N点DFT值，可采用任意数为基数的DFT算法，计算效率低于以2为基数FFT算法。2020/3/47•设N可以表示为若干因子的乘积LppppN321我们定义Lpppq321则11qpN这样，我们可把输入序列分解成p1个子序列，每个子序列由q1个取样组成。2020/3/4810)()(NnnkNWnxkX10)1(1110)1(11011111111)1()1()(qrkprpNqrkrpNqrrkpNWprpxWrpxWrpx1,...,2,1,0)(10110111NkWlrpxWqrrkpNpllkN因而2020/3/49利用关系式rkqrkpNWW11令1,...,2,1,0)()(110111qkWlrpxkGqrrkpNl因此10110111)()(qrrkqpllkNWlrpxWkX1,...,2,1,0)(101NkkGWlpllkN上式说明，当N为组合数时，序列x(n)的离散傅立叶变换X（k）可借助于P1个序列长度为其q1的DFT来表示。2020/3/4例：讨论一个长N＝18点序列x（n）的DFT快速计算问题。解：因为N＝18＝3×3×2令p1＝3，q1＝6首先将原始序列分成p1＝3个子序列，每个子序列长度q1=6个取样序列1（L=0）;x(0),x(3),x(6),x(9),x(12),x(15)序列2（L=1）;x(1),x(4),x(7),x(10),x(13),x(16)序列3（L=p1-1=2）;x(2),x(5),x(8),x(11),x(14),x(17)102020/3/41110110111)()(qrrkpNpllkNWlrpxWkX503182018)3(rrkllkWlrxW2018)(lllkkGW)()()(22181180kGWkGWkGkk。、、，他们分别记为点的＝上式的内和式是三个)()()(DFT61q210kGkGkG2020/3/412然后将每一个6点子序列分成三个长度为2的序列。即有506)3()(rrklWlrxkG1036506)39(ttksskWlstxW则X（k）可表示为1022062018)39()(ttksskllkWlstxWWkX2020/3/4三、实数序列的FFT13•一个N点FFT算法可同时计算两个N点实数序列的DFT•一个N点FFT一次计算一个2N点实数序列的DFT2020/3/4一个N点FFT算法可同时计算两个N点实数序列的DFT14FFT计算为复数运算，所以输入序列x(n)在运算时可以为复数数据。如果是实序列一般是把虚部置另。如果利用虚部数据可以提高计算效率。2020/3/4设)(nx和)(ny为实序列，组合)()()(njynxna为复序列。计算：)()()]([)(kjAkAnaDFTkAir)()()]([)]([)()()()]([)(101010kjYkXnyjDFTnxDFTWnyjWnxWnanaDFTkAnkNNnnkNNnnkNNn152020/3/416)()()()())(()()()()()()(***10*101010)(10)(10)(10kjYkXWnyjWnxWnyjWnxWnyjWnxWnakNAnkNNnnkNNnnkNNnnkNNnkNnNNnkNnNNnkNnNNn)()()(*kjYkXkNA2020/3/417)]()([2)]()([21)]()([2)(*kAkNAjkNAkAkNAkAjkYrrii)()()()()()(*kjYkXkAkjYkXkNA)]()([21)]()([21)]()([21)(*kNAkAjkNAkAkNAkAkXiirr2020/3/4一个N点FFT一次计算一个2N点实数序列的DFT设)(nx为2N点实数序列，进行奇偶抽取。1,....,2,1,0)12()()2()(Nnnxnhnxng即：1,...2,1,0)()()()()()]([)(10NkkjHkGnjAkAWnanaDFTkAirnkNNn181,....,2,1,0)()()(NNnnjhngna点的复数序列将其组成一个2020/3/41,....,2,1,0)()()()()()(22NkkHWkGkNXkHWkGkXkNkN＋根据蝶形运算的特点1,....,2,1,0)]()([2)]()([21)()]()([2)]()([21)(NkkAkNAjkNAkAkHkNAkAjkNAkAkGrriiiirr由前面的讨论，可知192020/3/420例已知，是两个N点实序列，的值，今需要从，求，的值，为了提高运算效率，试用一个N点运算一次完成。XkYkxnynDFTXkYkxnynIFFT2020/3/421解：由题意XkDFTxnYkDFTyn，构造序列ZkXkjYk对作一次N点IFFT可得序列Zkzn又根据DFT的线性性质IDFTXkjIDFTYk而，都是实序列xnynReImxnznynzn()znIDFTZk()znIDFTZkIDFTXkjYkxnjyn2020/3/4MATLAB提供了一个称为fft的函数用于计算一个向量X的DFT。调用X=fft(x,N)就计算出N点的DFT。如果向量X的长度小于N，那么就将X补零。如果N略去，那么DFT的长度就是X的长度。如果X是一个矩阵，那么fft(x,N)计算X中每一列的N点的DFT。2020/3/422这个fft是用机器语言写成的，而不是用MATLAB命令（也就是不是作为一个.m文件来使用的），因此执行起来非常快。并且它是作为一种混合基算法写成的。如果N是2的某个幂，那么就能使用一个高速的基-2FFT算法。如果N不是2的某个幂，那么就将N分解为若干因子并用一个较慢的混合基FFT算法。最后，如果N就是某个质数，那么fft函数就退化为原始的DFT算法。应用ifft函数计算逆DFT，它与fft具有相同的特性。2020/3/423例2-1对连续的单一频率周期信号按采样频率采样，截取长度N分别选N=20和N=16，观察其DFT结果的幅度谱。解此时离散序列，即k=8。用MATLAB计算并作图，函数fft用于计算离散傅里叶变换DFT，程序如下页：2020/3/424()sin(2/)sin(2/8)asxnnffn8saffk=8;%所以k=8n1=[0:1:19];%向量n1为0,1,2,…,19，截取长度N为20xa1=sin(2*pi*n1/k);%生成正弦周期序列xa1subplot(2,2,1)%指定画图左上角plot(n1,xa1)%画图xa1xlabel(‘t/T’);ylabel(‘x(n)’);%标x,y轴说明xk1=fft(xa1);xk1=abs(xk1);%计算向量xa1的FFTsubplot(2,2,2)%指定画图右上角bar(n1,xk1)%画xk1的柱状图，要提取幅度用bar(n1,xk1/20*2)代替xlabel('k');ylabel('X(k)');%标x,y轴说明n2=[0:1:15];%向量n1为0,1,2,…,15，截取长度N为16xa2=sin(2*pi*n2/k);%生成正弦周期序列xa1subplot(2,2,3)%指定画左下角plot(n2,xa2)%画图xa2xlabel('t/T');ylabel('x(n)');%标x,y轴说明xk2=fft(xa2);xk2=abs(xk2);%计算向量xk2的FFTsubplot(2,2,4)%指定画右下角bar(n2,xk2)%画xk1的柱状图，要提取幅度用bar(n2,xk2/16*2)代替xlabel('k');ylabel('X(k)');%标x,y轴说明2020/3/4258saff计算结果示于上图，左上图和右上图分别是N=20时的截取信号和DFT结果，由于截取了两个半周期，频谱出现泄漏；左下图和右下图分别是N=16时的截取信号和DFT结果，由于截取了两个整周期，得到单一谱线的频谱。上述频谱的误差主要是由于时域中对信号的非整周期截断产生的频谱泄漏。2020/3/42605101520-1-0.500.51t/Tx(n)-100102000.20.40.60.8kX(k)051015-1-0.500.51t/Tx(n)-100102000.511.5kX(k)