DSP第四章（长春工大李敏杰数字信号处理实验图像及书上完整版答案加考试范围）完整版共15篇-

1、第四章学习目标 理解按时间抽选的基-2FFT算法的算法原理、运 算流图、所需计算量和算法特点 理解按频率抽选的基-2FFT算法的算法原理、运 算流图、所需计算量和算法特点 理解IFFT算法 了解CZT算法 理解线性卷积的FFT算法及分段卷积方法 第四章 快速傅里叶变换 FFT: Fast Fourier Transform 1965年，Cooley, Tukey 机器计算傅里叶级数的一种算法 一、直接计算DFT的问题及改进途径 N 点有限长序列x(n DFT : nk X ( k = DFT x ( n = x ( n WN RN ( k n =0 N 1 IDFT : 1 x ( n = I

2、DFT X ( k = N X (k W k =0 N 1 nk N RN ( n 运算量 x(nW n =0 N 1 一个X(k nk N 复数乘法 复数加法 N N1 N2 N (N 1 N个X(k (N点DFT ( a + jb ( c + jd = ( ac bd + j ( ad + cb 一次复乘 一次复加 一个X (k N个X (k (N点DFT 4N 4N 2 实数乘法 4 实数加法 2 2 2N+2 (N 1=2 (2N 1 2N (2N 1 W 的特性 对称性 nk N nk WN = e j 2 nk N nk ( n (WN * = WN nk = WN N n k =

3、 WN ( N k nN nk Nk nk WN WN WN WN 周期性 可约性 W nk N =W ( N +n k N =W n ( N +k N nk mnk WN = WmN 2 j mnk mN nk nk / WN = WN / mm j 2 N N 2 e 0 ( k 特殊点： WN = 1 WNN / 2 = 1 WNk + N / 2 = WN e =e j = 1 FFT 算法的基本思想： 利用DFT 系数的特性，合并DFT 运算中的某些项， 把长序列DFT 短序列DFT，从而减少其运算量。 FFT算法分类: 按时间抽选法 DIT: Decimation-In-Time

4、按频率抽选法 DIF: Decimation-In-Frequency 二 、按时间抽选的基-2FFT算法 1、算法原理 设序列点数 N = 2L，L 为整数。 若不满足，则补零 N为2的整数幂的FFT算法称基-2FFT算法。 将序列x(n按n的奇偶分成两组： x ( 2r = x1 ( r r = 0,1,., N / 2 1 x ( 2r + 1 = x2 ( r 则x(n的DFT: nk nk nk X ( k = x ( n WN = x ( n WN + x ( n WN n =0 n =0 n =0 N 1 N 1 N 1 n为偶数 = = N / 21 r =0 n为奇数 ( 2

5、 r +1k N x ( 2r W r =0 1 2 rk N + rk N / 21 r =0 x ( 2r + 1W k N N / 2 1 r =0 2 N / 21 x ( r (W 2 N +W x ( r (W 2 N rk = N / 21 r =0 x1 ( r WNrk/ 2 + WNk N / 2 1 r =0 x2 ( r WNrk/ 2 k = X 1 ( k + WN X 2 ( k r , k = 0,1,.N / 2 1 再利用周期性求X(k的后半部分 Q X 1 ( k , X 2 ( k 是以N / 2为周期的 N X1 k + = X1 (k 2 又WN k

6、+ N 2 k k = WNN / 2WN = WN N X 2 k + = X 2 (k 2 k X ( k = X 1 ( k + WN X 2 ( k k = 0,1,., N / 2 1 N k X ( k + = X 1 ( k WN X 2 ( k 2 分解后的运算量： 复数乘法 一个N / 2点DFT (N / 22 两个N / 2点DFT N 2 / 2 一个蝶形 N / 2个蝶形 总计 1 N/2 N2 /2+ N /2 N2 /2 复数加法 N / 2 (N / 2 1 N (N / 2 1 2 N N ( N / 2 1 + N N2 /2 运算量减少了近一半 N / 2

7、仍为偶数，进一步分解：N / 2 x1 (2l = x3 (l x1 (2l + 1 = x4 (l N/4 l = 0,1,., N / 4 1 k X 1 (k = X 3 (k + WN / 2 X 4 (k N k = 0,1,., 1 N k 4 X 1 (k + 4 = X 3 (k WN / 2 X 4 (k 同理： k X 2 (k = X 5 (k + WN / 2 X 6 (k N k = 0,1,., 1 N 4 k X 2 (k + = X 5 (k WN / 2 X 6 (k 4 其中： X 5 (k = DFT x5 (l = DFT x2 (2l l = 0,1,

8、., N / 4 1 X 6 (k = DFT x6 (l = DFT x2 (2l + 1 k 2 统一系数：WN / 2 WN k 这样逐级分解，直到2点DFT 当N = 8时，即分解到X3(k，X4(k，X5(k， X6(k，k = 0, 1 X 3 (k = N / 4 1 l =0 x (l W 3 lk N /4 = x3 (l W l =0 1 lk N /4 k = 0,1 0 X 3 (0 = x3 (0W20 + W20 x3 (1 = x (0 + WN x (4 0 X 3 (1 = x3 (0W20 + W21 x3 (1 = x (0 WN x (4 X 4 (k

9、= N / 4 1 l =0 x4 (l W lk N /4 = x4 (l W l =0 1 lk N /4 k = 0,1 0 X 4 (0 = x4 (0W20 + W20 x4 (1 = x (2 + WN x (6 0 X 4 (1 = x4 (0W20 + W21 x4 (1 = x (2 WN x (6 2、运算量 当N = 2L时，共有L级蝶形，每级N / 2个蝶形，每 个蝶形有1次复数乘法2次复数加法。 N N 复数乘法： mF = L = log 2 N 2 2 复数加法： aF = NL = N log 2 N 比较DFT mF ( DFT N2 2N = = mF (

10、FFT N log N log 2 N 2 2 3、算法特点 1）原位运算 r X m (k = X m1 (k + X m1 ( j WN r X m ( j = X m1 (k X m1 ( j WN m表示第m级迭代，k，j表示数据所在的行数 2）倒位序 n0 0 n1 0 1 1 0 1 n2 0 1 0 1 0 1 0 1 x(n n = (n2 n1n0 2 倒位序 $ n = (n0 n1n2 2 000 100 010 110 001 101 011 111 0 4 2 6 1 5 3 7 自然序 n = (n2 n1n0 2 0 1 2 3 4 5 6 7 000 001 0

11、10 011 100 101 110 111 3）蝶形运算 对N = 2L点FFT，输入倒位序，输出自然序， 第m级运算每个蝶形的两节点距离为 2m1 第m级运算： r X m (k = X m1 (k + X m1 (k + 2m1 WN r X m (k + 2m1 = X m1 (k X m1 (k + 2m1 WN W r N 的确定 蝶形运算两节点的第一个节点为k值，表 示成L位二进制数，左移L m位，把右边 空出的位置补零，结果为r的二进制数。 r = (k 2 2 Lm 4）存储单元 输入序列x(n : N个存储单元 r WN ：N / 2个存储单元 系数 4、DIT算法的其他形

12、式流图 输入倒位序输出自然序 输入自然序输出倒位序 输入输出均自然序 相同几何形状 输入倒位序输出自然序 输入自然序输出倒位序 三 、按频率抽选的基-2FFT算法 1、算法原理 设序列点数N=2L，L为整数。 将X(k按k的奇偶分组前，先将输入x(n按n的 顺序分成前后两半： nk X (k = x(nWN = n =0 N 1 N / 2 1 n =0 nk x(nWN + n= N / 2 N 1 nk x(nWN = = N / 2 1 n =0 n= 0 x ( n W nk N + N / 2 1 n =0 N x n + W 2 N n + k 2 N N / 2 1 N Nk /

13、 2 nk x ( n + x n + 2 WN WN WNN / 2 = 1 = N / 2 1 n =0 N nk k x ( n + ( 1 x n + 2 WN k = 0,1,., N 1 按k的奇偶将X(k分成两部分： k = 2r k = 2r + 1 X (2r = N / 2 1 r = 0,1,., N / 2 1 N 2 nr x ( n + x n + 2 WN n =0 = N / 2 1 n =0 N nr x ( n + x n + 2 WN / 2 N n (2 r +1 x ( n x n + 2 WN N n nr x ( n x n + WN WN / 2

14、 2 X (2r + 1 = = N / 2 1 n =0 n =0 N / 2 1 令 N x1 ( n = x ( n + x n + 2 x (n = x(n x n + N W n N 2 2 N n = 0,1,., 1 2 则X(2r和X(2r+1分别是x1(n和x2(n的 N / 2 点DFT，记为X1(k和X2(k x(0 x(1 x(2 x(3 x(4 x(5 x(6 x(7 x1(0 x1(1 x1(2 x1(3 N/2点 DFT X1(0=X(0 X1(1=X(2 X1(2=X(4 X1(3=X(6 X2(0=X(1 N/2点 DFT X2(1=X(3 X2(2=X(5

15、X2(3=X(7 W -1 0 N x2(0 x2(1 x2(2 x2(3 1 WN -1 -1 -1 W 2 N 3 WN N /2仍为偶数，进一步分解：N /2 N /4 N x3 ( n = x1 ( n + x1 ( n + N / 4 n = 0,1,., 1 n 4 x4 ( n = x1 ( n x1 ( n + N / 4WN / 2 X 3 (k = X 1 (2k = DFT x3 (n X 4 (k = X 1 (2k + 1 = DFT x4 (n N k = 0,1,., 1 4 x1(0 x3(0 N/4点 DFT X3(0=X1(0=X(0 x1(1 x3(1 X

16、3(1=X1(2=X(4 x1(2 W -1 0 N /2 x4(0 N/4点 DFT X4(0=X1(1=X(2 x1(3 W -1 1 N /2 x4(1 X4(1=X1(3=X(6 同理： X 5 (k = X 2 (2k = DFT x5 (n X 6 (k = X 2 (2k + 1 = DFT x6 (n 其中： N k = 0,1,., 1 4 x5 ( n = x2 ( n + x2 ( n + N / 4 n x6 ( n = x2 ( n x2 ( n + N / 4WN / 2 N n = 0,1,., 1 4 逐级分解，直到2点DFT 当N=8时，即分解到x3(n，x4

17、(n，x5(n， x6(n，n=0,1 X 3 (k = N / 4 1 l =0 x (l W 3 lk N /4 = x3 (l W l =0 1 lk N /4 k = 0,1 X (0 = X 3 (0 = x3 (0W20 + W20 x3 (1 = x3 (0 + x3 (1 0 X (4 = X 3 (1 = x3 (0W20 + W21 x3 (1 = x3 (0 x3 (1WN X 4 (k = N / 4 1 l =0 x4 (l W lk N /4 = x4 (l W l =0 1 lk N /4 k = 0,1 X (2 = X 4 (0 = x4 (0W20 + W2

18、0 x4 (1 = x4 (0 + x4 (1 0 X (6 = X 4 (1 = x4 (0W20 + W21 x4 (1 = x4 (0 x4 (1WN 2、算法特点 1）原位运算 L级蝶形运算，每级N/2个蝶形，每个蝶形结构： X m (k = X m1 (k + X m1 ( j r X m ( j = X m1 (k X m1 ( j WN m表示第m级迭代，k，j表示数据所在的行数 X m1 (k X m 1 ( j r WN X m (k X m ( j -1 2）蝶形运算 对N=2L点FFT，输入自然序，输出倒位序， 两节点距离：2L-m=N / 2m 第m级运算： N X m

19、 ( k = X m 1 ( k + X m 1 k + m 2 X k + N = X (k X k + N W r m 1 m 1 N m 2m 2m r WN 的确定 蝶形运算两节点的第一个节点为k值，表示成L 位二进制数，左移m-1位，把右边空出的位置 补零，结果为r的二进制数。 r = (k 2 2 m 1 3、DIT与DIF的异同 基本蝶形不同 DIT: 先复乘后加减 DIF: 先减后复乘 运算量相同 N mF = log 2 N 2 aF = N log 2 N 都可原位运算 DIT和DIF的基本蝶形互为转置 四 、IFFT算法 比较： IDFT: DFT: 1 x(n = N

20、X ( k WN nk k =0 N 1 X (k = x(nW n =0 nk N N 1 nk N FFT : W W nk N 1 × N IFFT L 1 1 = N 2 1 N 1 nk x(n = X (k WN N n =0 1 N 1 * nk x* (n = X (k WN N k =0 * 1 1 * nk * x ( n = X ( k WN = DFT X ( k N k =0 N N 1 * 直接调用FFT子程序计算IFFT的方法： X (k 共轭 X * (k FFT 共轭 乘1/ N x ( n 例：设x (n 是2 N 点实数序列，试用一次N 点DFT

21、 来计算x (n 的2 N 点DFT: X (k 解：将x ( n 按奇偶分组，令 x1 ( n = x (2n x2 ( n = x (2n + 1 n = 0,1,., N 1 n = 0,1,., N 1 构成一个复序列 w( n = x1 ( n + jx2 ( n 对w( n 进行一次N 点DFT运算 W ( k = DFT w( n = X 1 ( k + jX 2 ( k 得 X 1 ( k = Wep ( k 1 X 2 ( k = Wop ( k j 均为N 点DFT 而X ( k 是2 N 点DFT ? 八 、线性调频 z变换（CZT算法 FFT不适用于： 只研究信号的某一

22、频段，要求对该频段抽样密 集，提高分辨率； 研究非单位圆上的抽样值； 需要准确计算N点DFT，且N为大的素数；等等。 CZT算法：对z变换采用螺线抽样，chirp-z变换 线性调频 z变换 1、算法原理 N点有限长序列，其z变换： X ( z = x(n z n =0 N 1 n 沿z平面上的一段螺线作等分角抽样，抽样点zk： zk = AW k 其中： z0 = A0e j0 k = 0,1,., M 1 W = W0 e j0 M为要分析的复频谱点数 则 zk = A0W0 k e j (0 + k0 抽样点： zk = A0W0 k e j (0 + k0 A0 ：起始抽样点的矢量半径长

23、度 0 ：起始抽样点的相角 0 ：相邻抽样点的角度差 0 > 0 ： 逆时针 W0 ：螺线的伸展率 W0>1：螺线内缩 W0<1： 螺线外伸 0 < 0 ：顺时针 当W0=1，则表示半径为A0的一段圆弧 若又有A0=1，则表示单位圆上的一段圆弧 若又有 0 = 0, 0 = 2 / N ，M=N ，即为序列的DFT。 求抽样点处的z变换： X ( zk = x ( n z n =0 N 1 n k = x ( n A W n n =0 N 1 nk k = 0,1,., M 1 NM次复乘 (N-1 M次 复加 由 nk = 1/ 2 n 2 + k 2 ( k n 2

24、 得 X ( zk = x ( n A nW W n =0 k 2 N 1 2 N 1 n2 2 ( k n 2 2 W k2 2 n = W x(n A W n =0 n2 2 W ( k n 2 2 n X ( zk = W x ( n A nW 2 n =0 k 2 N 1 2 2 ( k n W 2 2 令 g (n = x(n A W h( n = W n2 2 n n2 2 n = 0,1,., N 1 k2 2 则 X ( zk = W k 2 N 1 2 g ( n h ( k n = W g ( k * h ( k n =0 k = 0,1,., M 1 2、CZT的实现步骤

25、及运算量的估算 n X ( zk = W x ( n A W n =0 k 2 N 1 2 n2 2 W ( k n 2 2 其中 g ( n = x ( n A nW h( n = W n2 2 n 2 2 k = 0,1,., M 1 n = 0,1,., N 1 h( n : ( N 1 ( M 1 N + M 1点 g (n : 0 N 1 h( n * g (n : 2 N + M 2 Q X ( zk k = 0 M 1 L N + M 1 且L = 2 m 1 选择 L N + M 1 ，且 L = 2m 2 形成L点序列g(n： (3N n A W x ( n 0 n N 1

26、g ( n = 0 N n L 1 n2 2 系数 Cn = A W n n2 / 2 Cn +1 = A ( n +1 W ( n +12 / 2 =A W n n2 / 2 (W W n 1/ 2 A = Cn Dn 1 Dn = W nW 1/ 2 A1 = W W n 1W 1/ 2 A1 = WDn 1 D0 = W 1/ 2 A1 C0 = 1 求其L点FFT： L 1 j 2 rn L ( L/2*log2L 0 r L 1 G (r = FFT g (n = g (ne n =0 3）形成L点序列h(n: W h ( n = 0 ( L n 2 W 2 n2 2 (2N 0 n

27、 M 1 M n LN L N +1 n L 1 求其L点FFT: H (r = FFT h(n = h(ne n =0 L 1 (L/2*log2L j 2 rn L 0 r L 1 4）求乘积 Q(r = H (r G (r (L (L/2*log2L 5）求L点IFFT的 q (k 2 j rn 1 L1 q (k = IFFT Q (r = H (r G ( r e L L n =0 取 q( k = q( k RM ( k 6）求得抽样点的z变换: X ( zk = W q ( k k2 2 (M 0 k M 1 总运算量: 3 mF = L log 2 L + 5 N + L +

28、M 2 3、CZT算法的优点 N，M可为任意数，可不等 0 任意，易调整频率分辨率 周线可以是螺线，而不一定是圆弧 z0任意，从任意频率开始，便于窄带高分辨率分析 当 A = 1，M = N ，W = e j 2 N 时，CZT=DFT，解 决了N为素数的快速算法问题 九、线性卷积和线性相关的FFT算法 1、线性卷积的FFT算法 若L点x(n，M点h(n， 则直接计算其线性卷积y(n y (n = M 1 m =0 h(m x(n m 需运算量： md = LM 若系统满足线性相位，即： h ( n = ± h( M 1 n 则需运算量： md = LM / 2 FFT法：以圆周卷积

29、代替线性卷积 令 N = 2 M + L 1 x ( n 0 n L 1 x ( n = L n N 1 0 m h( n 0 n M 1 h ( n = M n N 1 0 则 y (n = x (n * h( n = x ( n N h(n 1 H(k = FFT h(n 2 X(k =FFT x(n 3 Y(k = H(kX(k 4 y(n = IFFT Y(k N /2*log2N N /2*log2N N N /2*log2N mF = N (1 + 3/ 2 * log 2 N 比较直接计算和FFT法计算的运算量 md ML Km = = mF 2 N (1 + 3/ 2* log 2 N