数字信号处理（西安理工）：DSP_5 (Chapter 4)

李洁《数字信号处理》2005? Digital Signal Processing__Jie Li 2005? 1 数字信号处理 第四章 快速 Fourier变换 (FFT) Digital Signal Processing 主讲教师：李洁 李洁李洁--《《数字信号处理数字信号处理》》--第四章第四章快速快速Fourier变换变换2 / 30 0 500 1000 1500 2000 0 10203040 Number of samples, N Number of Op era t ions DFT ∝ N 2 FFT ∝ N log 2 N 直直 接计算接计算 DFT 由于由于 W N kn 的的 周期性和对周期性和对 称性，会使计算过程出现许多冗余称性，会使计算过程出现许多冗余 直直接计算接计算DFT总共需要总共需要N 2 次次 复数乘法和复数乘法和N（（N-1）次复数）次复数 加法！加法！ V E R Y B A D ! 首先，首先， FFT并不是一种新的并不是一种新的 Fourier变换，只是变换，只是 DFT的一种快速算法。的一种快速算法。 1,,1,0)()]([)( 1 0 ?=== ∑ ? = NkWnxnxDFTkX N n kn N L W 8 0 W 8 4 W 8 2 W 8 5 W 8 6 W 8 3 W 8 7 W 8 1 512010485761024 44816384128 80102432 4164 Radix-2DFTN 李洁《数字信号处理》2005? Digital Signal Processing__Jie Li 2005? 2 李洁李洁--《《数字信号处理数字信号处理》》--第四章第四章快速快速Fourier变换变换3 / 30 例 李洁李洁--《《数字信号处理数字信号处理》》--第四章第四章快速快速Fourier变换变换4 / 30 § 4.2 基 2FFT算法 1. 的特性 m N W 周期性 22 ()jmlN jm mlN m NN NN We e W ππ ?+ ? + === 对称性 W 8 0 W 8 4 W 8 2 W 8 5 W 8 6 W 8 3 W 8 7 W 8 1 W 16 0 W 16 8 W 16 4 W 16 10 W 16 12 W 16 6 W 16 14 W 16 2 W 16 1 W 16 3 W 16 5 W 16 7 W 16 9 W 16 11 W 16 13 W 16 15 怎样减少怎样减少 DFT的计算量？的计算量？ WWWWWW m N N m N m N mN N mN N m N ? === +??? 2 )*( 可约性 WWWW mnk mN nk N mnk mN nk N / / == 当当 当当 当当 N=2 M 时，？？？时，？？？ 李洁《数字信号处理》2005? Digital Signal Processing__Jie Li 2005? 3 李洁李洁--《《数字信号处理数字信号处理》》--第四章第四章快速快速Fourier变换变换5 / 30 2. 时域抽取法基 2FFT基本原理 1,,1,0)()]([)( 1 0 ?=== ∑ ? = NkWnxnxDFTkX N n kn N L W 8 0 W 8 4 W 8 2 W 8 5 W 8 6 W 8 3 W 8 7 W 8 1 0 1 2 3 4 5 6 7 8 9 10 n 0 N x(n) 1 3,2,1,0,)2()( 1 == rrxrx 3,2,1,0,)12()( 2 =+= rrxrx 3,2,1,0,)()( 3 0 4 11 == ∑ = krxkX n kn W 3,2,1,0,)()( 3 0 4 22 == ∑ = krxkX n kn W WWWWWWWWW kkkkkkkk n kn xxxxxxxxnxkX 7 8 5 8 3 8 1 8 6 8 4 8 2 8 0 8 7 0 8 )7()5()3()1()6()4()2()0()()( +++++++== ∑ = WW kn knjknj kn ee 2 8 2 8 2 4 2 4 === ?? ππ ))(()( ))7()5()3()1(())6()4()2()0(( ))7()5()3()1(())6()4()2()0(( 2 8 1 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 kXkX xxxxxxxx xxxxxxxx W WWWWWWWWW WWWWWWWWW k kkkkkkkkk kkkkkkkkk += +++++++= +++++++= 3,2,1,0,))(()())4(()4()4( 3,2,1,0,))(()()( 2 8 12 4 8 1 2 8 1 =?=+++=+ =+= + kkXkXkXkXkX kkXkXkX WW W kk k 李洁李洁--《《数字信号处理数字信号处理》》--第四章第四章快速快速Fourier变换变换6 / 30 ? 蝶形 WN k X 1 (k) X 1(k)+WN k X 2(k) X 1 (k)-W N k X 2 (k) X2(k) Butterfly 李洁《数字信号处理》2005? Digital Signal Processing__Jie Li 2005? 4 李洁李洁--《《数字信号处理数字信号处理》》--第四章第四章快速快速Fourier变换变换7 / 30 李洁李洁--《《数字信号处理数字信号处理》》--第四章第四章快速快速Fourier变换变换8 / 30 李洁《数字信号处理》2005? Digital Signal Processing__Jie Li 2005? 5 李洁李洁--《《数字信号处理数字信号处理》》--第四章第四章快速快速Fourier变换变换9 / 30 李洁李洁--《《数字信号处理数字信号处理》》--第四章第四章快速快速Fourier变换变换10 / 30 4. DIT-FFT的运算规律及编程思想 Ⅰ）原位计算 李洁《数字信号处理》2005? Digital Signal Processing__Jie Li 2005? 6 李洁李洁--《《数字信号处理数字信号处理》》--第四章第四章快速快速Fourier变换变换11 / 30 Ⅱ）旋转因子的变化规律 p N W N点点DIT――FFT运算流图中，每级都有运算流图中，每级都有N/2个蝶形。个蝶形。 每个蝶形都要乘以因子每个蝶形都要乘以因子W N p ，称其为旋转因子，，称其为旋转因子，p称为旋转因子的指数。称为旋转因子的指数。 观察图不难发现，第观察图不难发现，第L级共有级共有2 L-1 个不同的旋转因子。个不同的旋转因子。 LM Jp ? ?= 2 12,...,2,1,0 1 ?= ?L J 李洁李洁--《《数字信号处理数字信号处理》》--第四章第四章快速快速Fourier变换变换12 / 30 Ⅲ）蝶形运算规律 1 2 ? = L B 同一旋转因子对应着间隔为同一旋转因子对应着间隔为 2 L 点的点的 2 M-L 个蝶形。个蝶形。 李洁《数字信号处理》2005? Digital Signal Processing__Jie Li 2005? 7 MATLAB Ⅳ）编程思想及程序流程 从输入端第1级开始，逐级进 行，共M级，所以L=1:M 第L级中每个蝶形 相距B=2 L-1 个点 第L级共有2 L-1 个不 同的旋转因子， J=[0:B-1]则旋转 因子指数分别为 P=2 M-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 = 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 李洁《数字信号处理》2005? Digital Signal Processing__Jie Li 2005? 8 李洁李洁--《《数字信号处理数字信号处理》》--第四章第四章快速快速Fourier变换变换15 / 30 (1)从输入端第从输入端第 1级开始，逐级进行，共级开始，逐级进行，共 M级，所以级，所以 L=1:M (2)第第 L级中每个蝶形相距级中每个蝶形相距 B=2 L-1 个点个点 (3)第第 L级共有级共有 2 L-1 个不同的旋转因子，个不同的旋转因子， J=[0:B-1]则旋则旋 转因子指数分别为转因子指数分别为 P=2 M-L *J (4)使用同一个旋转因子的蝶形相距使用同一个旋转因子的蝶形相距 2 L 个点，所以同个点，所以同 一组蝶形的起始下标为一组蝶形的起始下标为 k=J:2 L :N-1 流图解读流图解读 李洁李洁--《《数字信号处理数字信号处理》》--第四章第四章快速快速Fourier变换变换16 / 30 Ⅴ）序列的倒序 ? 倒序的规律 你能写出 N=16时的倒序数吗？ 规律规律 1：写成：写成 M位位 2进制数后，将其倒置，再转换成相应的十进制数进制数后，将其倒置，再转换成相应的十进制数 李洁《数字信号处理》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变换变换18 / 30 7 6 5 4 3 2 1 0 n J 7 J 6 J 5 J 4 J 3 J 2 J 1 J 0 J n N/2 N/4 N/8 李洁《数字信号处理》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<J) Temp=X(I+1); X(I+1)=X(J+1); X(J+1)=Temp; end K=LH; while(J>=K) J=J-K; K=K/2; end J=J+K; end J赋初值 J 1 =0+(N/2) J 0 和J N 不用 交换所以 I=[1：N-2] 为了保证只交换一 次，仅在I<J时交换 由J i 计算J i+1 李洁李洁--《《数字信号处理数字信号处理》》--第四章第四章快速快速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 李洁《数字信号处理》2005? Digital Signal Processing__Jie Li 2005? 11 李洁李洁--《《数字信号处理数字信号处理》》--第四章第四章快速快速Fourier变换变换21 / 30 5. 频率抽取法 FFT （ DIF-FFT） 李洁李洁--《《数字信号处理数字信号处理》》--第四章第四章快速快速Fourier变换变换22 / 30 李洁《数字信号处理》2005? Digital Signal Processing__Jie Li 2005? 12 李洁李洁--《《数字信号处理数字信号处理》》--第四章第四章快速快速Fourier变换变换23 / 30 李洁李洁--《《数字信号处理数字信号处理》》--第四章第四章快速快速Fourier变换变换24 / 30 6. FFT的变形运算流图 李洁《数字信号处理》2005? Digital Signal Processing__Jie Li 2005? 13 李洁李洁--《《数字信号处理数字信号处理》》--第四章第四章快速快速Fourier变换变换25 / 30 6. FFT的变形运算流图 李洁李洁--《《数字信号处理数字信号处理》》--第四章第四章快速快速Fourier变换变换26 / 30 7.快速逆 DFT ∑ ? = == 1 0 )()]([)( N n kn N W nxnxDFTkX ∑ ? = ? == 1 0 )( 1 )]([)( N n kn N W kX N kXIDFTnx ∑∑ ? = ? = ? == 1 0 1 0 )(* 1 )*)(* 1 )(* ( N n kn N N n kn N W W kX N kX N nx *)]}(*[{ 1 ]*)(*[ 1 )( 1 0 kXDFT N kX N nx N n kn N W == ∑ ? = 上式两边同时取共轭上式两边同时取共轭 书上有错！ 李洁《数字信号处理》2005? Digital Signal Processing__Jie Li 2005? 14 MATLAB 7. IFFT function X=ifftM(X) %利用 FFT实现 IDFT的子程序 [N,o]=size(X); if (~(rem(log2(N),1)==0)) error('序列长度不是 2的整数次幂！不能进行 IFFT！ ') end X=conj(fft(conj(X)))/N; 习题5 *)]}(*[{ 1 ]*)(*[ 1 )( 1 0 kXDFT N kX N nx N n kn N W == ∑ ? = 李洁李洁--《《数字信号处理数字信号处理》》--第四章第四章快速快速Fourier变换变换28 / 30 § 4.3 实序列的 FFT 李洁《数字信号处理》2005? Digital Signal Processing__Jie Li 2005? 15 习题 习题 4 Page 习题开讲 习题 习题 3 Page 已知已知X(k)和和Y(k)分别是两个分别是两个N点实序列点实序列x(n)和和 y(n)的的N点点DFT，若要求，若要求x(n)和和y(n)，为提高运，为提高运 算效率，试设计用一次算效率，试设计用一次N点点IFFT来完成。来完成。 )()()()()( kFkFkjYkXkF opep +=+= 习题开讲 )](Im[)](Re[)]([)( nfjnfkFIFFTnf +== )()](Re[)]([)]([ nxnfkFIDFTkXIDFT ep === )()](Im[)]([)]([ njynfjkFIDFTkjYIDFT op === )](*)([ 2 1 )( nfnfnx += )](*)([ 2 1 )( nfnf j ny ?=