快速傅里叶变换课件.ppt

1、快速傅里叶变换快速傅里叶变换2n直接计算直接计算DFT的问题及改进的途径的问题及改进的途径 3 4 设复序列设复序列x(n)长度为长度为N点，其点，其DFT为为10()()NnkNnX kx n Wk=0，N-1（1）计算一个）计算一个X(k)值的运算量值的运算量复数乘法次数：复数乘法次数：N复数加法次数：复数加法次数：N1（2）计算全部）计算全部N个个X(k)值的运算量值的运算量复数乘法次数：复数乘法次数：N2复数加法次数：复数加法次数：N(N1)5N点点DFT的复数乘法次数举例的复数乘法次数举例NN2NN22464404941612816384864256 65 536 16256512

2、262 144 3210241024 1 048 576 结论结论：当：当N很大时，其运算量很大，对实时性很强的信号很大时，其运算量很大，对实时性很强的信号处理来说，要求计算速度快，因此需要改进处理来说，要求计算速度快，因此需要改进DFT的计算的计算方法，以大大减少运算次数。方法，以大大减少运算次数。6 nkNW主要原理是利用系数主要原理是利用系数 的以下特性对的以下特性对DFT进行分解：进行分解：nkNW（2）对称性）对称性()nkNW()k N nNW（1）周期性）周期性()()n N kn kNnkNNNWWW（3）可约性）可约性 mnknkmNNWW/nknk mNN mWW另外，另外

3、，12/NNWkNNkNWW)2/(7DIT-FFT(Decimation-In-Time)1(2)()xrx r设设N2M，将，将x(n)按按 n 的奇偶分为两组：的奇偶分为两组：2(21)()xrx rr=0，1，12N10()()()NnkNnX kDFT x nx n W则则1010)()(NnnnkNNnnnkNWnxWnx为奇数为偶数8120)12(1202)12()2(NrkrNNrrkNWrxWrx1202212021)()(NrrkNkNNrrkNWrxWWrx)()(21kXWkXkN1100()()()NNnknkNNnnnnX kx n Wx n W为偶数为奇数式中，式

4、中，X1(k)和和X2(k)分别是分别是x1(n)和和x2(n)的的N/2的的DFT。另外，式中另外，式中k的取值范围是：的取值范围是：0，1，N/21。910因此，因此，只能计算出只能计算出X(k)的前一半值。的前一半值。12()()()kNX kX kW Xk后一半后一半X(k)值，值，N/2，N/2 1，N-1？1()2NXk2 1(2)120()Nr NkNrx r W2 1120()NrkNrx r W1()X k同理可得同理可得22()()2NXkXk11(2)NkkNNWW 因此可得后半部分因此可得后半部分X(k)2()2()2(221NkXWNkXNkXNkN12()()()k

5、NX kX kW Xk由前半部分由前半部分X(k)()(21kXWkXkNk=0，1，N/21k=0，1，N/2112k=k=0 0，1 1，N/2N/21 1结论：结论：12()()()kNX kXkW Xk12()()()2kNNX kX kW Xk2 11120()()NrkNrXkx r W2 12220()()NrkNrXkx r W因此，只要因此，只要求出求出2 2个个N/2N/2点的点的DFTDFT，即，即X X1 1(k)(k)和和X X2 2(k)(k)，再经过这种再经过这种运算就可求运算就可求出全部出全部X X(k)(k)的值的值1312()()()kNX kX kW Xk

6、12()()()2kNNX kX kW Xk蝶形运算式蝶形运算式蝶形运算蝶形运算信号流图符号信号流图符号 X1(k)X2(k)蝶形运算的运算量：每次蝶形含一次复数乘和两蝶形运算的运算量：每次蝶形含一次复数乘和两 次复数加次复数加0NW1NW2NW3NW以以N=8为例，为例，分解为分解为2个个4点点的的DFT，然后，然后做做8/2=4次蝶形次蝶形运算即可求出运算即可求出所有所有8点点X(k)的的值。值。1415复数乘法次数：复数乘法次数：N2复数加法次数：复数加法次数：N(N1)复数乘法次数：复数乘法次数：2*(N/2)2+N/2=N2/2+N/2复数加法次数：复数加法次数：2*(N/2)(N/

7、21)+2*N/2=N2/2nN点 每次蝶形含一次复数每次蝶形含一次复数乘和两次复数加乘和两次复数加16 由于由于N2M，因而，因而N/2仍是偶数仍是偶数，可以进一步把每个，可以进一步把每个N/2点点子序列再按其奇偶部分分解为两个子序列再按其奇偶部分分解为两个N/4点的子序列。点的子序列。以以N/2点序列点序列x1(r)为例为例 1314(2)()0,1,1(21)()4xlx lNlxlx l则有则有 rkNNrWrxkX212011)()(klNNllkNNlWlxWlx)12(21401221401)12()2(lkNNlkNlkNNlWlxWWlx41404241403)()()()(

8、42/3kXWkXkNk=0,1,14N17且且13/24()()4kNNXkXkWXkk=0,1,14N由此可见，一个由此可见，一个N/2点点DFT可分解成两个可分解成两个N/4点点DFT。同理，也可对同理，也可对x2(n)进行同样的分解，求出进行同样的分解，求出X2(k)。18概念：信号流图概念：信号流图100(0)()(0)(1)(0)(1)NnkNNnXx n WxxxW x10()()0,1NnkNnX kx n Wk 对此例对此例N=8，最后剩下的是，最后剩下的是4个个N/4=2点的点的DFT，2点点DFT也可以由蝶形运算来完成。也可以由蝶形运算来完成。N=2这说明，这说明，N=2

9、N=2M M的的DFTDFT可全部由蝶形运算来完成。可全部由蝶形运算来完成。因此：因此：N N必须是必须是2 2的整数次幂的整数次幂.“.“基基2”2”100(1)()(0)(1)(0)(1)NnkNNnXx n WxxxW x蝶形运算蝶形运算1920N=8按时间抽取法按时间抽取法FFT信号流图信号流图 x(n)抽取抽取X(k)顺序顺序DIT-FFT(Decimation-In-Time)FFT算法算法21由按时间抽取法FFT的信号流图可知，当N=2M时，共有 级蝶形运算；每级都由 个蝶形运算组成，而每个蝶形有 次复乘、次复加，因此每级运算都需 次复乘和 次复加。MN/2 N/2 12N22F

10、FT算法中，这样算法中，这样 级运算总共需要：级运算总共需要：M复数乘法：复数乘法：2log22NNMN复数加法：复数加法：2logN MNN直接直接DFT算法运算量算法运算量 复数乘法：复数乘法：复数加法：复数加法：N2N(N1)直接计算直接计算DFT与与FFT算法的计算量之比为算法的计算量之比为MNNNNNM222log2log2N=2M23NN2计算量之比M NN2计算量之比M 2414.012816 38444836.641644.025665 5361 02464.0864125.4512262 1442 304113.816256328.010241 048 5765 120204

11、.83210288012.820484 194 30411 264372.464404919221.4NN2log2NN2log22425N=1024;M=80;x=1:M,zeros(1,N-M);t=cputime;y1=fft(x,N);time_fft=cputime-t;t1=cputime;y2=dft(x,N);time_dft=cputime-t1;time_dft=6.0290time_fft=0.010026L=1L=2L=3倒倒序序27n同址运算（原位运算）n序列的逆序排列n蝶形运算两节点间的距离n相同旋转因子节点间的距离n (旋转因子）的确定pNW28 某一列任何两个节

12、点某一列任何两个节点k 和和j 的节点变量进行蝶形运算的节点变量进行蝶形运算后，得到结果为下一列后，得到结果为下一列k、j两节点的节点变量，而和其他两节点的节点变量，而和其他节点变量无关。这种原位运算结构可以节省存储单元，节点变量无关。这种原位运算结构可以节省存储单元，降低设备成本。降低设备成本。)(kX)2(NkX)(1kX)(2kXkNW运算前运算前运算后运算后()A j)(kA()A j)(kA例例n同址运算（原位运算）2930)(01221)()(BINMMDECnnnnnn 由于由于 x(n)被反复地按奇、偶分组，所以流图输被反复地按奇、偶分组，所以流图输入端的入端的排列不再是顺序的

13、，但仍有规律可循：排列不再是顺序的，但仍有规律可循：因为因为 N=2M，对于任意对于任意 n（0n N-1)，可以用可以用M个个二进制码表示为：二进制码表示为：10,01221nnnnnMM n 反复按奇、偶分解时，即按二进制码的反复按奇、偶分解时，即按二进制码的“0”“1”分解。分解。n序列的逆序排列3132自然顺序自然顺序 n二进制数二进制数倒位序二进制数倒位序二进制数倒位序顺序数倒位序顺序数0000000010011004201001023011110641000011510110156110011371111117 nIJ规律：最左端加规律：最左端加1，向右进位，向右进位倒序倒序33I

14、=J 不换不换IJ 不换不换IX=czt(x,M,W,A)n其中，其中，x是待变换的时域信号是待变换的时域信号x(n)，其长度为，其长度为N，M是是变换长度，变换长度，W确定变换的步长，确定变换的步长，A确定变换的起点。确定变换的起点。若若M=N，A=1，则，则CZT变成变成DFT。7610()()()()()Nmy nx nh nx m h nm设一离散线性移不变系统的冲激响应为设一离散线性移不变系统的冲激响应为 ,其输其输入信号为入信号为 .其输出为其输出为 .并且并且 的长度为的长度为N N点点,的长度为的长度为M M点点,计算其线性卷积。计算其线性卷积。)(nx)(nx)(nh)(nh

15、)(ny)(ny)(nx)(nh77用FFT算 的步骤：)(n ny y1.(),()1;x n h nLNM将补零点，至少为点2.()();H kFFT h nLFFT求,点3.()()X kFFT x n求,L点FFT；求)()()(.4kHkXkY。求)()(.5kYIFFTny78FFTFFTIFFTxx(n)h(n)y(n)X(k)H(k)Y(k)流程图流程图程序参考第三章幻灯片程序参考第三章幻灯片79例例4.8.2 设模拟信号设模拟信号 ，以，以 t=0.01n(n=0:N-1)进行取样，试用进行取样，试用fft函数对其做频谱分析。函数对其做频谱分析。N分别分别为：为：(1)N=1

16、0；(2)N=50；(3)N=100；(2)N=1024。()2sin(20)5cos(40)x ttt程序清单如下程序清单如下%计算计算N=10的的FFT并绘出其幅频曲线并绘出其幅频曲线N=10;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(20*pi*t)+5*cos(40*pi*t);y=fft(x,N);figure(1)subplot(2,2,1)plot(q,abs(y)title(FFT N=10)80%计算计算N=50的的FFT并绘出其幅频曲线并绘出其幅频曲线N=50;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(20*pi*

17、t)+5*cos(40*pi*t);y=fft(x,N);figure(1)subplot(2,2,2)plot(q,abs(y)title(FFT N=50)81%计算计算N=100的的FFT并绘出其幅频曲线并绘出其幅频曲线N=100;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(20*pi*t)+5*cos(40*pi*t);y=fft(x,N);figure(1)subplot(2,2,3)plot(q,abs(y)title(FFT N=100)82%计算计算N=1024的的FFT并绘出其幅频曲线并绘出其幅频曲线N=1024;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(20*pi*t)+5*cos(40*pi*t);y=fft(x,N);figure(1)subplot(2,2,4)plot(q,abs(y)title(FFT N=1024)Thank You世界触手可及世界触手可及携手共进，齐创精品工程携手共进，齐创精品工程