快速傅里叶变换

1、2021/3/231n自己能编写自己能编写FFT算法算法n理论上理解理论上理解FFT算法算法2021/3/232n直接计算直接计算DFT的问题及改进的途径的问题及改进的途径 2021/3/233 2021/3/234 设复序列设复序列x(n) 长度为长度为N点点,其其DFT为为10( )( )NnkNnX kx n Wk=0,N-1 （1）计算一个）计算一个X(k) 值的运算量值的运算量复数乘法次数复数乘法次数: N复数加法次数复数加法次数: N1（2）计算全部）计算全部N个个X(k) 值的运算量值的运算量复数乘法次数复数乘法次数: N2复数加法次数复数加法次数: N(N1)2021/3/23

2、5N点点DFT的复数乘法次数举例的复数乘法次数举例NN2NN22464404941612816384864256 65 536 16256512 262 144 3210241024 1 048 576 结论结论:当当N很大时很大时,其运算量很大其运算量很大,对实时性很强的信号处理对实时性很强的信号处理来说来说,要求计算速度快要求计算速度快,因此需要改进因此需要改进DFT的计算方法的计算方法,以大以大大减少运算次数。大减少运算次数。 2021/3/236 nkNW主要原理是利用系数主要原理是利用系数 的以下特性对的以下特性对DFT进行分解进行分解: nkNW（2）对称性）对称性 ()nkNW(

3、)k N nNW（1）周期性）周期性 ()()n N kn kNnkNNNWWW（3）可约性）可约性 mnknkmNNWW/nknk mNN mWW另外另外,12/NNWkNNkNWW)2/(2021/3/237DIT-FFT(Decimation-In-Time)2021/3/2381(2 )( )xrx r设设N2M,将将x(n)按按 n 的奇偶分为两组的奇偶分为两组: 2(21)( )xrx rr =0,1, 12N10( ) ( )( )NnkNnX kDFT x nx n W则则1010)()(NnnnkNNnnnkNWnxWnx为奇数为偶数2021/3/239120)12(1202

4、) 12()2(NrkrNNrrkNWrxWrx1202212021)()(NrrkNkNNrrkNWrxWWrx)()(21kXWkXkN1100( )( )( )NNnknkNNnnnnX kx n Wx n W为偶数为奇数式中式中,X1(k)和和X2(k)分别是分别是x1(n)和和x2(n)的的N/2的的DFT。另外另外,式中式中k的取值范围是的取值范围是:0,1, ,N/21 。2021/3/2310因此因此, 只能计算出只能计算出X(k)的前一半值。的前一半值。12( )( )( )kNX kX kW Xk后一半后一半X(k) 值值, N/2 , N/2 1, ,N-1 ?1()2N

5、Xk2 1(2)120( )Nr NkNrx r W2 1120( )NrkNrx r W1( )X k同理可得同理可得22()( )2NXkXk2021/3/2311(2)NkkNNWW 因此可得后半部分因此可得后半部分X(k) )2()2()2(221NkXWNkXNkXNkN12( )( )( )kNX kX kW Xk由前半部分由前半部分X(k) )()(21kXWkXkNk=0,1, ,N/21k=0,1, ,N/212021/3/2312k=k=0,10,1, , , ,N/2N/21 1结论结论:12( )( )( )kNX kXkW Xk12()( )( )2kNNX kX k

6、W 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)的值的值2021/3/231312( )( )( )kNX kX kW Xk12()( )( )2kNNX kX kW Xk蝶形运算式蝶形运算式蝶形运算蝶形运算信号流图符号信号流图符号 X1(k)X2(k)蝶形运算的运算量蝶形运算的运算量:每次蝶形含一次复数乘和两每次蝶形含一次复数

7、乘和两 次复数加次复数加2021/3/23140NW1NW2NW3NW以以N=8为例为例,分分解为解为2个个4点的点的DFT,然后做然后做8/2=4次蝶形运次蝶形运算即可求出所算即可求出所有有8点点X(k)的值。的值。2021/3/2315复数乘法次数复数乘法次数: N2复数加法次数复数加法次数: N(N1)复数乘法次数复数乘法次数: 2*(N/2)2+N/2=N2/2+N/2复数加法次数复数加法次数: 2*(N/2)(N/21)+2*N/2=N2/2nN点 每次蝶形含一次复数每次蝶形含一次复数乘和两次复数加乘和两次复数加2021/3/2316 由于由于N2M,因而因而N/2仍是偶数仍是偶数

8、,可以进一步把每个可以进一步把每个N/2点点子序列再按其奇偶部分分解为两个子序列再按其奇偶部分分解为两个N/4点的子序列。点的子序列。 以以N/2点序列点序列x1(r)为例为例 1314(2 )( )0,1,1(21)( )4xlx lNlxlx l则有则有 rkNNrWrxkX212011)()(klNNllkNNlWlxWlx)12(21401221401) 12()2(lkNNlkNlkNNlWlxWWlx41404241403)()()()(42/3kXWkXkNk=0,1, 14N2021/3/2317且且13/24( )( )4kNNXkXkWXkk=0,1, 14N由此可见由此可

9、见,一个一个N/2点点DFT可分解成两个可分解成两个N/4点点DFT。 同理同理,也可对也可对x2(n)进行同样的分解进行同样的分解,求出求出X2(k)。 2021/3/2318概念概念:信号流图信号流图2021/3/2319100(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,N=2M M的的DFTDFT可全部由蝶形运算来完成。可全部由蝶形运算来完成。因

10、此因此:N:N必须是必须是2 2的整数次幂的整数次幂. . “基基2 2”100(1)( )(0)(1)(0)(1)NnkNNnXx n WxxxW x蝶形运算蝶形运算2021/3/2320N=8按时间抽取法按时间抽取法FFT信号流图信号流图 x(n)抽取抽取X(k)顺序顺序DIT-FFT(Decimation-In-Time)FFT算法算法2021/3/2321由按时间抽取法FFT的信号流图可知,当N=2M时,共有 级蝶形运算;每级都由 个蝶形运算组成,而每个蝶形有 次复乘、 次复加,因此每级运算都需 次复乘和 次复加。 MN/2 N/2 12N2021/3/2322FFT算法中算法中,这样

11、这样 级运算总共需要级运算总共需要: M复数乘法复数乘法: 2log22NNMN复数加法复数加法: 2logN MNN直接直接DFT算法运算量算法运算量 复数乘法复数乘法: 复数加法复数加法: N2N(N1)直接计算直接计算DFT与与FFT算法的计算量之比为算法的计算量之比为MNNNNNM222log2log2N=2M2021/3/2323NN2计算量之比M NN2计算量之比M 2414.012816 38444836.641644.025665 5361 02464.0864125.4512262 1442 304113.816256328.010241 048 5765 120204.83

12、210288012.820484 194 30411 264372.464404919221.4NN2log2NN2log22021/3/23242021/3/2325N=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.01002021/3/2326L=1L=2L=3倒倒序序2021/3/2327n同址运算（原位运算）n序列的逆序排列n蝶形运算两节点间的距离n相同旋

13、转因子节点间的距离n (旋转因子）的确定pNW2021/3/2328 某一列任何两个节点某一列任何两个节点k 和和j 的节点变量进行蝶形运算的节点变量进行蝶形运算后后,得到结果为下一列得到结果为下一列k、j两节点的节点变量两节点的节点变量,而和其他节而和其他节点变量无关。这种原位运算结构可以节省存储单元点变量无关。这种原位运算结构可以节省存储单元,降低降低设备成本。设备成本。)(kX)2(NkX)(1kX)(2kXkNW运算前运算前运算后运算后( )A j)(kA( )A j )(kA例例n同址运算（原位运算）2021/3/23292021/3/2330)(01221)()(BINMMDECn

14、nnnnn 由于由于 x(n) 被反复地按奇、偶分组被反复地按奇、偶分组,所以流图输所以流图输入端的入端的排列不再是顺序的排列不再是顺序的,但仍有规律可循但仍有规律可循: 因为因为 N=2M , 对于任意对于任意 n（0n N-1),可以用可以用M个个二进制码表示为二进制码表示为:10,01221nnnnnMM n 反复按奇、偶分解时反复按奇、偶分解时,即按二进制码的即按二进制码的“0” “1” 分解。分解。n序列的逆序排列2021/3/23312021/3/2332自然顺序自然顺序 n二进制数二进制数倒位序二进制数倒位序二进制数倒位序顺序数倒位序顺序数00000000100110042010

15、01023011110641000011510110156110011371111117 nIJ规律规律:最左端加最左端加1,向右进位向右进位倒序倒序2021/3/2333I=J 不换不换IJ 不换不换IX= czt(x, M, W, A)n其中其中,x是待变换的时域信号是待变换的时域信号x(n),其长度为其长度为N,M是变换是变换长度长度,W确定变换的步长确定变换的步长,A确定变换的起点。若确定变换的起点。若M= N,A= 1,则则CZT变成变成DFT。2021/3/237610( )( )( )( ) ()Nmy nx nh nx m h nm设一离散线性移不变系统的冲激响应为设一离散线性

16、移不变系统的冲激响应为 ,其输其输入信号为入信号为 .其输出为其输出为 .并且并且 的长度为的长度为N N点点, 的长度为的长度为M M点点,计算其线性卷积。计算其线性卷积。)(nx)(nx)(nh)(nh)(ny)(ny)(nx)(nh2021/3/2377用FFT算 的步骤: )(n ny y1.( ), ( )1 ;x n h nLNM将补零点，至少为点2.( )( );H kFFT h nLFFT求, 点3.( )( )X kFFT x n求,L点FFT；求)()()(.4kHkXkY。求)()(.5kYIFFTny2021/3/2378FFTFFTIFFTxx(n)h(n)y(n)X

17、(k)H(k)Y(k)流程图流程图程序参考第三章幻灯片程序参考第三章幻灯片2021/3/2379例例4.8.2 设模拟信号设模拟信号 ,以以 t= 0.01n (n=0: N-1) 进行取样进行取样,试用试用fft函数对其做频谱分析。函数对其做频谱分析。N分别为分别为:(1) N=10;(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)2021/3/2380%计算计算N=50的的FFT并绘出其幅频曲线并绘出其幅频曲线N=50;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,2)plot(q,abs(y)title(FFT N=50)2021/3/2381%计算计算N=100的的FFT并绘出其幅频曲线并绘出其幅频曲线N=100;n=0:N-1;t=0.0