李洁《数字信号处理》20058 数字信号处理 Digital Signal processing 第四章快速 Fourier变换(FFT 主讲教师:李洁 形影x去 首先,FFT并不是一种新的 Fourier变换,只是DFT的一种快速算法 X(k)=DFx(m=2∑xmHk=01…,N-1,we 直接计算D由W如的周期性和对 称性,会使计算过程出现许多冗余 直接计算DFT总共需要N次 复数乘法和N(N-1)次复数o 加法! FTRadix-2 FT∝N FFT o Nlog2N 1281638448 102410485765120 李洁一《数字信号处理 2|30 Digital Signal processing Jie Li 2005
李洁《数字信号处理》2005® Digital Signal Processing__Jie Li 2005® 1 数字信号处理 第四章 快速Fourier变换(FFT) Digital Signal Processing 主讲教师：李洁 李洁 -- 《数字信号处理》 -- 第四章 快速Fourier Fourier变换 2 / 30 0 500 1000 1500 2000 0 10 20 30 40 Number of samples, N Number of Operations DFT ∝ N2 FFT ∝ N log2N 直接计算DFT DFT 由于WN kn 的周期性和对 称性，会使计算过程出现许多冗余 直接计算DFT总共需要N2次 复数乘法和N（N-1）次复数 加法！ VERY BAD VERY BAD ! 首先，FFT并不是一种新的Fourier变换，只是DFT的一种快速算法。 ( ) [ ( )] ( ) 0,1, , 1 1 0 = = ∑ = − − = X k DFT x n x n W k N N n kn N L W8 0 W8 4 W8 2 W8 5 W8 6 W8 3 W8 7 W8 1 1024 1048576 5120 128 16384 448 32 1024 80 4 16 4 N DFT Radix-2
李洁《数字信号处理》20058 如果一台通用计算机的速度为平均每次复乘5s,每次复加0.5s,用它来计算 的DFTx(n)],问直接计算需要多少时间,用FFT运算需要多少时间 ①直接利用DFT计算:复乘次数为N2复加次数为N(N-1) ②利用FFT计算:复乘次数为 复加次数为Nog:N (1)直接计算 T=5×10×N=5x104×5122=1.31072s 复加所需时间 =0.5×10×N×(N-1)=0.5×10×512×(512 T-T1+T2=1.441536s 复乘所需时间 T=5×10×2N-5×10+×2912-0.012 复加所需时间 所T1=0.5×10×NogN=0.5×10×512512=0.00304 T=T+T:=0.013824s 李洁一《数字信号处理 3/ 影x去 4.2基2FF算法当当当N-2时,?? 1.H的特性怎神少D的计算量 周期性Wm+=e -(m+IN 对称性W=W"(W")=WW=-W果 可约性W=WmW=Wm 李洁一《数字信号处理 4l30 Digital Signal processing Jie Li 2005
李洁《数字信号处理》2005® Digital Signal Processing__Jie Li 2005® 2 李洁 -- 《数字信号处理》 -- 第四章 快速Fourier Fourier变换 3 / 30 例 李洁 -- 《数字信号处理》 -- 第四章 快速Fourier Fourier变换 4 / 30 §4.2 基2FFT算法 1. m 的特性 WN 周期性 2 2 j m lN j m ( ) m lN m N N We e W N N π π −+ − + = == 对称性 W8 0 W8 4 W8 2 W8 5 W8 6 W8 3 W8 7 W8 1 W160 W168 W164 W1610 W1612 W166 W1614 W162 W161 W16 W 3 165 W167 W169 W1611 W1613 W1615 怎样减少DFT的计算量？ W W W W W W m N N m N m N N m N N m N m N = = = − − − − + ( )* 2 可约性 W W W W nk m N m nk N mnk mN nk N / / = = 当 当 当N=2M时，？？？
李洁《数字信号处理》20058 2.时域抽取法基2FF基本原理 X(k)=DFx(m)=∑x(n)W如k=01…N-1 012345678910 x(r)=x(2r),r=0,2,3x1(k)=∑x)4,k=012,3 x4k)=∑x(r)H+, X(k)=∑xn)W=x0)W“+x(2)W+x(4W+x16)1+x(DW+x(3)W+x5)H1+x7W =(x(0)H“+x02)2+x(4)W“+x(6)W4)+W(x()W"+x(3)W+x5W+x7)W =(x(0)W“+x(2)W+x4)W“+x(6)W“)+“(x(D)W“+x3)W“+x(5H:+x7)W2 x(k)=x1(k)+W:(x3(k),k=0123 X(+4)=x1(k+4)+W(x2(k+4)=X1(k)-H,(x2(k),k=0123 5 影x去 蝶形 Butterfly 李洁一《数字信号处理 630 Digital Signal processing Jie Li 2005
李洁《数字信号处理》2005® Digital Signal Processing__Jie Li 2005® 3 李洁 -- 《数字信号处理》 -- 第四章 快速Fourier Fourier变换 5 / 30 2. 时域抽取法基2FFT基本原理 ( ) [ ( )] ( ) 0,1, , 1 1 0 = = ∑ = − − = X k DFT x n x n W k N N n kn N L W8 0 W8 4 W8 2 W8 5 W8 6 W8 3 W8 7 W8 1 0 1 2 3 4 5 6 7 8 9 10 n 0 N x(n) 1 ( ) (2 ), 0,1,2,3 x1 r = x r r = ( ) (2 1), 0,1,2,3 x2 r = x r + r = ( ) ( ) , 0,1,2,3 3 0 1 1 4 = ∑ = = X k x r k n kn W ( ) ( ) , 0,1,2,3 3 0 2 2 4 = ∑ = = X k x r k n kn W W W W W W W W W W k k k k k k k k n kn X k x n x x x x x x x x 7 8 5 8 3 8 1 8 6 8 4 8 2 8 0 8 7 0 8 ( ) = ∑ ( ) = (0) + (2) + (4) + (6) + (1) + (3) + (5) + (7) = W W kn j kn j kn kn e e 2 8 2 8 2 4 2 4 = = = − − π π ( ) ( ( )) ( (0) (2) (4) (6) ) ( (1) (3) (5) (7) ) ( (0) (2) (4) (6) ) ( (1) (3) (5) (7) ) 1 8 2 3 4 2 4 1 4 0 4 1 8 3 4 2 4 1 4 0 4 6 8 4 8 2 8 0 8 1 8 6 8 4 8 2 8 0 8 X k X k x x x x x x x x x x x x x x x x W W W W W W W W W W W W W W W W W W W k k k k k k k k k k k k k k k k k k k = + = + + + + + + + = + + + + + + + ( 4) ( 4) ( ( 4)) ( ) ( ( )), 0,1,2,3 ( ) ( ) ( ( )), 0,1,2,3 2 1 8 2 4 1 8 1 8 2 + = + + + = − = = + = + X k X k X k X k X k k X k X k X k k W W W k k k 李洁 -- 《数字信号处理》 -- 第四章 快速Fourier Fourier变换 6 / 30 蝶形 WN k X1(k) X1(k)+WN k X2(k) X1(k)-WN k X2(k) X2(k) Butterfly
李洁《数字信号处理》20058 X(0) N2点x1(1) X1(2) X(2) x(6) X1(3) X(3) x(1) 2点x2(1) X2(2) x(5) DFTX (3)wN 71 影x去 X X1(1) x(2) DFT X(① X2(0) X, (1 (3) 李洁一《数字信号处理 8|30 Digital Signal processing Jie Li 2005
李洁《数字信号处理》2005® Digital Signal Processing__Jie Li 2005® 4 李洁 -- 《数字信号处理》 -- 第四章 快速Fourier Fourier变换 7 / 30 李洁 -- 《数字信号处理》 -- 第四章 快速Fourier Fourier变换 8 / 30
李洁《数字信号处理》20058 李洁一《数字信号处理 9 影x去 4.DIT-FFT的运算规律及编程思想 I)原位计算 X(3) 李洁一《数字信号处理 10/30 Digital Signal processing Jie Li 2005 5
李洁《数字信号处理》2005® Digital Signal Processing__Jie Li 2005® 5 李洁 -- 《数字信号处理》 -- 第四章 快速Fourier Fourier变换 9 / 30 李洁 -- 《数字信号处理》 -- 第四章 快速Fourier Fourier变换 10 / 30 4. DIT-FFT的运算规律及编程思想 Ⅰ）原位计算
李洁《数字信号处理》20058 垂1x去 Ⅱ)旋转因子的变化规律 观察图不难发现,第L级共有21个不同的旋转因子,,P=小24∠ J=0,1,2, L-1-1 A0)x(0) 43)、4 A(7) A(7) (7) N点DT一FFT运算流图中,每级都有N2个蝶形 每个蝶形都要乘以因子WP,称其为旋转因子,p称为旋转因子的指数 李洁一《数字信号处理 130 影x去 Ⅲ)蝶形运算规律 x1()←Xz-1(J)+X1-1(+B)WR X,(+ B)eXi-1J)-XL-(+ B)W& 式中p=J·2M-;J=0,1,…,22-1-1;L=1,2,…,M x(0-4(0 B 0) X(1) (2)4(2) (3)4(6) 同一旋转因子对应间隔为2点的2M个蝶形。 李洁一《数字信号处理 12/30 Digital Signal processing Jie Li 2005
李洁《数字信号处理》2005® Digital Signal Processing__Jie Li 2005® 6 李洁 -- 《数字信号处理》 -- 第四章 快速Fourier Fourier变换 11 / 30 Ⅱ）旋转因子的变化规律 p WN N点DIT―FFT运算流图中，每级都有N/2个蝶形。 每个蝶形都要乘以因子WN p，称其为旋转因子，p称为旋转因子的指数。 观察图不难发现，第L级共有2 L-1个不同的旋转因子。 M L p J − = ⋅ 2 0,1,2,...,2 1 1 = − L− J 李洁 -- 《数字信号处理》 -- 第四章 快速Fourier Fourier变换 12 / 30 Ⅲ）蝶形运算规律 1 2 − = L B 同一旋转因子对应着间隔为2L点的2M-L个蝶形
李洁《数字信号处理》20058 ATLAB Ⅳ)编程思想及程序流程 function x= Myfft(x) [Xold, M, N]=Pinv(x)i 输入端第1级开始,逐级进 for I=l: M 行,共M级,所以L=1:M for J=0: B-1 第L级中每个蝶形 相距B=21个点 Np=exp(-1*2*pi*P/N); for kJJ: 2.L: N-1 x(k+1)=Xold(k+1)+Xold(k+B+l).*wNp; X(k+B+1)=Xo1d(k+1)-Xo1d(k+B+1) end k=1,81, 第L级共有21个不 Xold=X 同的旋转因子 end J=[0:B-们则旋转 因子指数分别为 个旋 P=2M-L*J 转因子的蝶 形相距2个 MATLAB function X= Myfft(x) TXold, M,N=Pinv(x); for l=l: M 八 for J=0: B-1 P=(2(ML)J; WNp=exp(-]*2*pi"P/N); for k=J: 2.L: N-1 X(k+l)=Xold(k+1)+Xold(k+B+l). * WNp; X(k+B+l)=Xold(k+1)-Xold(k+B+1). WNp; end XoldX d Digital Signal processing Jie Li 2005 7
李洁《数字信号处理》2005® Digital Signal Processing__Jie Li 2005® 7 MATLAB Ⅳ）编程思想及程序流程 从输入端第1级开始，逐级进 行，共M级，所以L=1:M 第L级中每个蝶形 相距B=2L-1个点 第L级共有2L-1个不 同的旋转因子， J=[0:B-1]则旋转 因子指数分别为 P=2M-L*J 使用同一个旋 转因子的蝶 形相距2 L 个 点 function X = Myfft(x) [Xold,M,N]=Pinv(x); for L=1:M B=2.^(L-1); for J=0:B-1 P=(2.^(M-L)).*J; WNp=exp(-j*2*pi*P/N); for k=J:2.^L:N-1 X(k+1)=Xold(k+1)+Xold(k+B+1).*WNp; X(k+B+1)=Xold(k+1)-Xold(k+B+1).*WNp; end end Xold=X; end MATLAB function X = function X = Myfft(x Myfft(x) [Xold,M,N]=Pinv(x); for L=1:M B=2.^(L-1); for J=0:B-1 P=(2.^(M-L)).*J; L)).*J; WNp=exp(-j*2*pi*P/N); j*2*pi*P/N); for k=J:2.^L:N-1 X(k+1)=Xold(k+1)+Xold(k+B+1).* X(k+1)=Xold(k+1)+Xold(k+B+1).*WNp; X(k+B+1)=Xold(k+1) X(k+B+1)=Xold(k+1)-Xold(k+B+1).* Xold(k+B+1).*WNp; end end Xold=X; end
李洁《数字信号处理》20058 (1)从输入端第1级开始,逐级进行,共级,所以 L=1:M (2)第L级中每个蝶形相距B=2-个点 (3)第L级共有2-个不同的旋转因子,J[0:1则旋 转因子指数分别为P=2J (4)使用同一个旋转因子的蝶形相距2个点,所以同 一组蝶形的起始下标为k=:24:N1 流图解读 k=,N-1,2 )=X(k)+X X(+B)e X()-X(+ B)WA 李洁一《数字信号处理 1530 影x去 V)序列的倒序 你能写出N=16时的倒序数吗? ■倒序的规律 规律1:写成M位2进制数后,将其倒置,再转换成相应的十进制数 0 0000 0 0102 1106 0 (n2 n, no)2 0011 0113 1117 李洁一《数字信号处理 1630 Digital Signal processing Jie Li 2005
李洁《数字信号处理》2005® Digital Signal Processing__Jie Li 2005® 8 李洁 -- 《数字信号处理》 -- 第四章 快速Fourier Fourier变换 15 / 30 (1)从输入端第1级开始，逐级进行，共M级，所以 L=1:M (2)第L级中每个蝶形相距B=2L-1个点 (3)第L级共有2L-1个不同的旋转因子，J=[0:B-1]则旋 转因子指数分别为P=2M-L*J (4)使用同一个旋转因子的蝶形相距2L个点，所以同 一组蝶形的起始下标为k=J:2L:N-1 流图解读 李洁 -- 《数字信号处理》 -- 第四章 快速Fourier Fourier变换 16 / 30 Ⅴ）序列的倒序 倒序的规律 你能写出N=16时的倒序数吗？ 规律1：写成M位2进制数后，将其倒置，再转换成相应的十进制数
李洁《数字信号处理》20058 MATLAB function [x, M, N]= Pinv(x) [N, o]=size(x) M=1og2(N); fox工=1:N-2 serl=sscanf (dec2bin(I),'Bc,inf while (M=Ns) serl=strcat('0', serl) [ms, Ns]=size(serl); ser3=' for K-M: -1:1 ser3=strcat(ser3, serl(l,k)) end ser=bin2dec(ser3) x(工+1)=temg 影x去 顺序和倒序二进制数对照表 顺序 倒序 n十进制数1二进制数二进制数十进制数了L 0 000 0000-1 N/2N/4 N/8 04261 JJJ 010 110 6 111 111 JJJJ 11 李洁一《数字信号处理 16 8/30 Digital Signal processing Jie Li 2005
李洁《数字信号处理》2005® Digital Signal Processing__Jie Li 2005® 9 MATLAB function [x,M,N]= Pinv(x) [N,o]=size(x); M=log2(N); for I=1:N-2 ser1=sscanf(dec2bin(I),'%c',inf); [ms,Ns]=size(ser1); while (M~=Ns) ser1=strcat('0',ser1); [ms,Ns]=size(ser1); end ser3=''; for K=M:-1:1 ser3=strcat(ser3,ser1(1,K)); end ser=bin2dec(ser3); if (I<ser) temp=x(ser+1); x(ser+1)=x(I+1); x(I+1)=temp; end end 李洁 -- 《数字信号处理》 -- 第四章 快速Fourier Fourier变换 18 / 30 7 6 5 4 3 2 1 0 n J7 J6 J5 J4 J3 J2 J1 J0 Jn N/2 N/4 N/8
李洁《数字信号处理》20058 MATLAB V)序列的倒序 function [X,M, N]= invM (X) IN,o]=size(x) J赋初值 LH=N/2 J=0+(N2) J=lH: 交换所以 N=N-2 for I=1: N-2 [1:N2] i￡f(工<J) 为了保证只交换 次,仅在<J时交换 x(工+1)=X(J+1); X(J+1)=Tem 由计算J end K=工H J←J= K=K/2 JEJ+K 5.频率抽取法FFT(DIF-FFT) x()=2(x)+x(n+ N X()=∑x1(n)W M/2点 x(r)-r(n)-rl x(2r+1)= x2(n)wnra 李洁一《数字信号处理 0/30 Digital Signal processing Jie Li 2005 10
李洁《数字信号处理》2005® Digital Signal Processing__Jie Li 2005® 10 MATLAB Ⅴ）序列的倒序 function [X,M,N]= invM(X) [N,o]=size(X); LH=N/2; J=LH; for I=1:N-2 if(I=K) J=J-K; K=K/2; end J=J+K; end J赋初值 J1=0+(N/2) J0和JN不用 交换所以 I=[1：N-2] 为了保证只交换一 次，仅在I<J时交换 由Ji 计算Ji+1 李洁 -- 《数字信号处理》 -- 第四章 快速Fourier Fourier变换 20 / 30 5. 频率抽取法FFT（DIF-FFT） 0 1 2 3 4 5 6 7 8 9 10 n 0 N x(n) 1 N/2点DFT