数字信号处理第四章

2020 4 7 信息与通信工程学院 1 第四章快速傅立叶变换 FFT FastFourierTransform 4 1引言 4 2按时间抽取的基 2FFT算法 4 3按频率抽取的基 2FFT算法 4 4离散傅立叶反变换的快速计算方法 IFFT 本章主要内容 2020 4 7 2 4 1引言 FFT不是新的变换形式 只是DFT的一种快速算法 且根据对序列分解与选取方法的不同而产生了FFT的多种算法 应用很广 一DFT直接计算存在的问题设x n 为N点有限长序列 正反变换计算量相同 例N 4 信息与通信工程学院 2020 4 7 3 例N 4 N点DFT的直接计算量 1 乘法次数 对每一个k N次复数乘法 4N次实乘和2N次实加 N个k 次复数乘法 信息与通信工程学院 2020 4 7 4 加法次数 对每一个k N 1次复加 2 N 1 次实加 N个k N N 1 次复加即和成正比 例N 1024 则有1048576次复乘 约400万次实乘 假定运算器的指令速度为100MIPS 则计算时间大约为4秒 信息与通信工程学院 2020 4 7 5 旋转因子的对称性和周期性 例N 4 旋转因子W的矩阵形式为 二 改进途径 信息与通信工程学院 2020 4 7 6 旋转因子W的性质 得 可约性 信息与通信工程学院 2020 4 7 7 如前所述 N点DFT的复乘次数与N的平方成比例 显然N较小时乘法次数大大减少 利用上述旋转因子的特性 可以将有些项合并 并将DFT分解为短序列 从而降低运算次数 提高运算速度 1965年 库利 cooley 和图基 Tukey 首先提出FFT算法 对于N点DFT 仅需 N 2 log2N次复数乘法运算 例如N 1024 210时 需要 1024 2 log2210 512 10 5120次复乘运算 5120 1048576 4 88 速度提高20倍 分解有时域抽取 DIT 和频域抽取 DIF 两大类 改进方法概述 信息与通信工程学院 2020 4 7 8 4 2按时间抽取的基 2FFT算法 DecimationInTimeFFT 仍以N 4为例 运算流图 信息与通信工程学院 2020 4 7 9 一算法原理 先将x n 按n的奇偶分为两组作DFT 设N 2L L大于等于2 不足时 可补些零 这样一个N点的DFT分解为两个N 2点的DFT 一 N 2点DFT 信息与通信工程学院 2020 4 7 10 DIT FFT分解说明 信息与通信工程学院 2020 4 7 11 X k 后一半的确定 由旋转因子的周期性 这就是说 E k 和F k 的后一半分别等于其前一半的值 同理 可见 X k 的后一半 也完全由E k 和F k 的前一半确定 即N点的DFT可由两个N 2点的DFT来计算 信息与通信工程学院 2020 4 7 12 蝶形运算 流图表示 蝶形运算流图 N 2个蝶形 信息与通信工程学院 2020 4 7 13 二 N 4点DFT 由于N 2L 所以N 2仍为偶数 可以进一步把每个N 2点的序列再按其奇偶部分分解为两个N 4的子序列 其中 蝶形图同上类似 信息与通信工程学院 2020 4 7 14 三 2点DFT 按上述办法不断划分下去 直到最后剩下2点DFT 2点DFT实际上只有加减运算 DIT FFT算法 是在时间上对输入序列的次序是属于偶数还是属于奇数来进行分解的 所以称作按时间抽取的算法 例 作N 8时FFT的运算流图由上推证 令与 信息与通信工程学院 2020 4 7 15 例续 N 8时FFT的运算 展开得 信息与通信工程学院 2020 4 7 16 例续 N 8点FFT的运算流图 信息与通信工程学院 2020 4 7 17 二DIT FFT运算量分析 由上例分析可知 N 8 需三级蝶形运算 由此可推 N 2L共需L级蝶形运算 L log2N 每级都由N 2个蝶形运算组成 每个蝶形运算有一次复乘 两次复加 因此 N点的FFT的运算量为 复乘 mF N 2 L N 2 log2N复加 aF 2 N 2 L N log2N即计算量远比DFT的计算量小 以乘法作为基准 二者的比较结果如p 101图4 2 5所示 N越大 FFT的优点越明显 信息与通信工程学院 2020 4 7 18 原位运算 输入数据 中间结果和输出均用同一存储器 信息与通信工程学院 2020 4 7 19 通用蝶形图 可见 在某列进行蝶形运算的任意两个节点 行 k和j的节点变量就完全可以确定蝶形运算的结果 与其它行 节点 无关 这样 蝶形运算的两个输出值仍可放回蝶形运算的两个输入所在的存储器中 即实现所谓原位运算 每一组 列 有N 2个蝶形运算 所以只需N个存储单元 可以节省存储单元 信息与通信工程学院 2020 4 7 20 倒位序规律 由图知输出X k 按自然顺序排列在存储单元 而输入是按下列顺序排列 这种顺序称作倒位序 即二进制数倒位 造成倒位序的原因在于输入x n 按标号n的奇偶不断分组 倒位序实际是将序号n表示成L位的二进制数 把此二进制数位序颠倒 即得对应倒位序数 例 N 8 L 3 取x n 中的n 6 二进制表示为 110 倒位序得 011 即3 见p104表4 2 1 为避免倒位过程中重复变换位置 则设原序号为n 倒序后为 仅当时 与两值互换位置 否则不变换 信息与通信工程学院 2020 4 7 21 级 与 组 或群 N点FFT共有L级迭代 流图中的每一列为一级迭代 每一级中 由交叉蝶形组成一个群 第一级有N 2个群 第二级N 4个群 依次类推 最后一级1个群 蝶形运算两节点的 距离 指一个蝶形中 两输入或输出节点间间隔的点数 其值为2m 1 其中 m表示第m级 且m 1 L 第一级距离为1 第二级为2 第三级为4 第四级为8 依次类推 最后一级为N 2 信息与通信工程学院 2020 4 7 22 WNr因子的确定 由图可见 第一级 第二级 第三级 第i级 第L级 由旋转因子的可约性 个因子可表示为 信息与通信工程学院 2020 4 7 23 四DIT FFT其它形式流图 输入自然顺序 输出倒位序的DIT FFT运算流图 pp107图4 2 13 信息与通信工程学院 2020 4 7 24 输入输出均为自然顺序的DIT FFT运算流图 pp108图4 2 15 信息与通信工程学院 2020 4 7 25 五编程思想及程序框图 p104 105 倒序程序框图 DIT FFT运算程序框图 信息与通信工程学院 2020 4 7 26 3 3按频率抽取的基 2FFT算法 一算法原理N点DFT的另一种表达形式 DIF DecimationInFrequencyFFT 信息与通信工程学院 2020 4 7 27 N点DFT按k的奇偶分组分为两个N 2点的DFT 当k为偶数 即k 2r时 1 k 1 当k为奇数 即k 2r 1时 1 k 1 这时X k 可分为两部分 可见 上面两式均为N 2的DFT 信息与通信工程学院 2020 4 7 28 蝶形运算 令 则蝶形运算图如下 信息与通信工程学院 2020 4 7 29 再将N 2点DFT按k的奇偶各分解为两个N 4点的DFT 令r 2l与r 2l 1 信息与通信工程学院 2020 4 7 30 2点DFT 依上述方法 按k的奇偶不断分解L 1次 直至最终只有2点的DFT这种分解方法是先把输入序列分成前后两半 然后把输出序列按序号k的奇偶进行分解为短序列 即在频域进行抽取 简称DIF FFT 信息与通信工程学院 2020 4 7 31 二例如N 8时DIF FFT的运算规则 信息与通信工程学院 2020 4 7 32 例续N 8时DIF的运算规则 信息与通信工程学院 2020 4 7 33 例续N 8时DIF的运算规则 信息与通信工程学院 2020 4 7 34 例续N 8时DIF FFT的信号流图 信息与通信工程学院 2020 4 7 35 三DIF FFT的运算量分析 复乘 复加 即和DIT FFT的计算量相同 信息与通信工程学院 2020 4 7 36 四DIF FFT的运算特点 原位运算蝶形运算的通用表达式 通用蝶形图 信息与通信工程学院 2020 4 7 37 DIF FFT的运算特点 续 码位倒序 同DIT FFT 级 与 群 同DIT FFT 蝶形节点间距离 一般公式为2L m N 2m例如N 23 8 1 m 1时的距离为8 2 4 2 m 2时的距离为8 4 2 3 m 3时的距离为8 8 1 信息与通信工程学院 2020 4 7 38 DIF FFT的运算特点 续 旋转因子W的分布 由于DIF蝶形运算的两节点的距离为N 2m 所以蝶形运算可表示为 r的求法 k n2n1n0 左移m 1位 右边空出补零 得 r 2 k 22m 1 例如 N 8 1 m 1 k 2 k 010 2左移0位 r 2 010 2 2 2 m 2 k 1 k 001 2左移1位 r 2 010 2 2 3 m 2 k 5 k 101 2左移1位 r 2 010 2 2 信息与通信工程学院 2020 4 7 39 五DIF法与DIT法的异同 相同点 1 进行原位运算 2 运算量相同 复乘 N 2 Log2N次 复加NLog2N次 不同点 1 DIT输入为倒位序 输出为自然顺序 DIF正好与此相反 但DIT也有输入为自然顺序 输出为倒位序的情况 2 蝶形运算不同 见下页图 两种蝶形运算流图互为转置 将流图的所有支路方向都反向 交换输入和输出 节点变量值不变 基本蝶形的复乘 DIT在输入端 先乘后加 DIF在输出端 先加后乘 信息与通信工程学院 2020 4 7 40 DIT法与DIF法的蝶形运算比较 矩阵表示 信息与通信工程学院 2020 4 7 41 六DIF FFT的其它形式流图1 DIF FFT运算流图 N 8 信息与通信工程学院 2020 4 7 42 DIF FFT的其它形式流图2 DIT FFT的一种变形运算流图 N 8 信息与通信工程学院 2020 4 7 43 p127 1 设x n 是长度为2N的有限长实序列 X k 为x n 的2N点DFT 要求 试设计用一次N点FFT完成计算X k 的高效算法 解 本题的解题思路是DIT FFT思想 在时域分别抽取偶数点和奇数点x n 得到两个N点实序列x1 n 和x2 n x1 n x 2n n 0 1 N 1x2 n x 2n 1 n 0 1 N 1根据DIT FFT的思想 只要求得x1 n 和x2 n 的N点DFT 再经过简单的一级蝶形运算就可得到x n 的2N点DFT 因为x1 n 和x2 n 均为实序列 所以根据DFT的共轭对称性 可用一次N点FFT求得X1 k 和X2 k 做法如下 信息与通信工程学院 2020 4 7 44 1 解续 令y n x1 n jx2 n Y k DFT y n 则X1 k DFT x1 n Yep k Y k Y N k 2jX2 k DFT x2 n Yop k Y k Y N k 22N点DFT x n X k 可由X1 k 和X2 k 得到 X k X1 k W2NkX2 k X k N X1 k W2NkX2 k k 0 1 N 1这样 通过一次N点FFT计算完成了计算2N点DFT 信息与通信工程学院 2020 4 7 45 2 设x n 是长度为2N的有限长实序列 X k 为x n 的2N点DFT 若已知X k 要求 试设计用一次N点IFFT完成计算x n 的高效算法 解 设x1 n x 2n n 0 1 N 1x2 n x 2n 1 n 0 1 N 1X1 k DFT x1 n k 1 1 N 1X2 k DFT x2 n k 1 1 N 1由DIT FFT的算法 有以下关系式 X k X1 k W2NkX2 k k 1 1 N 1X k N X1 k W2NkX2 k 由以上两可解出X1 k 和X2 k 信息与通信工程学院 2020 4 7 46 2 解续 X1 k X k X k N 2X2 k W2N k X k X k N 2由上面分析可得出运算过程如下 1 由X k 计算出X1 k 和X2 k 2 由X1 k 和X2 k 构成N点频域序列Y k Y k X1 k jX2 k Yep k Yop k 其中Yep k X1 k Yop k jX2 k 进行N点IFFT运算 得到y n IFFT Y k Re y n jIm y n n 0 1 N 1由DFT的共轭对称性Re y n y n y n 2 IDFT Yep k x1 n jIm y n y n y n 2 IDFT Yop k jx2 n 由x1 n 和x2 n 合成x n 当n 偶数 n 0 1 2N 1 时x n x1 n 2 当n 奇数 n 0 1 2N 1 时x n x2 n 1 2 信息与通信工程学院 2020 4 7 47 4 4离散傅立叶反变换的快速计算方法 IFFT 一变动FFT程序和参数法实现IFFT 比较两式可知 只要DFT的每个系数换成 最后再乘以常数1 N就可以得到IDFT的快速算法 IFFT 可以将常数1 N分配到每级运算中 也就是每级蝶形运算均乘以1 2 信息与通信工程学院 2020 4 7 48 DIT IFFT的运算流图 防止溢出 信息与通信工程学院 2020 4 7 49 二利用FFT的程序直接实现IFFT 即先将X k 取共轭 将X k 的虚部乘 1 直接利用FFT程序计算DFT 然后再取一次共轭 最后再乘1 N 即为所求的 所以FFT IFFT可共用一个子程序 信息与通信工程学院 2020 4 7 50 三实序列的FFT算法 设x n 为N点实序列 取x n 的偶数点和奇数点分别作为新构造序列y n 的实部和虚部 即 对y n 进行N 2点FFT 输出Y k 则 由DIT FFT的思想得 由于x n 为实序列 所以X k 具有共轭对称性 X k 的另外N 2点的值为 信息与通信工程学院 2020 4 7 51 4 10FFT的应用 一线性卷积的长度设一离散线性移不变系统的冲激响应为h n 其输入信号为x n 其输出为y n 并且x n 的长度为L点 h n 的长度为M点 FIR 则 信息与通信工程学院 2020 4 7 52 线性卷积例子 信息与通信工程学院 2020 4 7 53 例续 线性卷积 信息与通信工程学院 2020 4 7 54 例续 线性卷积 信息与通信工程学院 2020 4 7 55 例续 线性卷积 结论 线性卷积结果的长度为L M 1 信息与通信工程学院 2020 4 7 56 二用FFT计算线性卷积的步骤 将x n h n 补零至长为N L M 1则 说明 当M L时 用圆周卷积计算线性卷积的速度快 点数越多 速度越快 N 64时 速度增加明显 L M时 速度增加不太明显 此时宜采用分段卷积 信息与通信工程学院 2020 4 7 57 三x n 为长序列时线性卷积的FFT算法 重叠相加法 设h n 点数为M 信号x n 为很长的输入序列 将x n 分解为很多段 每段长为L L和M的数量级相同 用xi n 表示x n 的第i段 长度为L M 1 信息与通信工程学院 2020 4 7 58 重叠相加法图示 由于yi n 长度为N L M 1 所以xi n 和h n 需要补零至相同长度N 因此相邻两段输出序列必然有N L个点重叠 最终结果y n 需要相加 故称重叠相加法 仅当N L M 1时才有M 1 N L 信息与通信工程学院 2020 4 7 59 用FFT实现重叠相加法的步骤 选L 使 将x n 分解成多个长为L的子序列xi n 之和 将xi n 和h n 补零至长为N的序列 分别求H k 与Xi k Yi k H k Xi k 并求yi k IDFT Yi k 将yi n 的后N L个样值与yi 1 n 的前N L个样值重叠相加 信息与通信工程学院 2020 4 7 60 重叠相加法的例子 信息与通信工程学院 2020 4 7 61 用表格法求解8点圆周卷积 信息与通信工程学院 2020 4 7 62 重叠保留法 重叠保留法的出发点是由于圆周卷积中发生混叠时只影响序列的一部分 如x1 n h n 的长度分别为L和M 则线性卷积的长度为L M 1 若L M 1 N 则N点圆周卷积的前L M 1 N个值不等于线性卷积 而后边余下的值相等 若通过 在序列前边 补数使序列长度增加 并把卷积后 在序列前边 多余的序列值舍去 则可得线性卷积结果 其它作法与重叠相加法相同 证明见教材p181 182 信息与通信工程学院 2020 4 7 63 重叠保留法的步骤 设h n 长为M x n 分段后每段长为L 选L 使 x0 n x0 n 是x n 前边补M 1个零值 凑N L M 1个点 x1 n 是由x0 n 后边取M 1个数值 再取其后x n 值凑够N点 如此反复取xi n 求 将yi n 的前M 1个值舍去 余下的值顺序排列为y n 信息与通信工程学院 2020 4 7 64 重叠保留法例1 x n 分为两段 每段长为L 6 M 3 设N 8 则 信息与通信工程学院 2020 4 7 65 重叠保留法例2 若选L 2 M 3 N M L 1 4 则有 信息与通信工程学院 2020 4 7 66 四线性相关的FFT算法 自相关与互相关在统计通信与数字信号处理中应用较广 同利用圆周卷积计算线性卷积一样 相关的计算也可以利用圆周相关来代替线性相关 称快速相关 线性相关 圆周相关 信息与通信工程学院 2020 4 7 67 快速相关的算法步骤 设x n 长为L y n 长为M 选 求 求 求 信息与通信工程学院 2020 4 7 68 4 9线性调频z变换 Chirp z变换 算法 一DFT算法及存在问题DFT算法 DFT是序列x n 的z变换X z 在z平面单位圆上N个等间隔点的抽样值 即整个单位圆上 0 2 的频谱 可用FFT计算 DFT算法存在的问题 FFT要求N 2m或高度复合数 N r1 r2 r3 无法求解到单位圆外的变换值 若信号为窄带信号 则不得不算许多冗余值 线性调频z变换是解决这一问题的适用方法 它是在Z平面上采用螺旋抽样 Chirp zTransform CZT 信息与通信工程学院 2020 4 7 69 二Chirp z变换的原理 设x n 为有限长 N点 序列 其z变换为 序列x n 长度为N 要分析z平面上M点频谱采样值 分析点为zk k 0 1 2 M 1 设zk AW k 0 k M 1式中A和W为复数 用极坐标形式表示为A A0ej 0W W0e j 0 zk A0W0 kej 0 k 0 式中 A0和W0为实数 当k 0 有z0 A0ej 0则zk表示z平面上起始半径为A0 角度为 0的螺旋线上角度间隔为 0的抽样点 如图所示 信息与通信工程学院 2020 4 7 70 Chirp Z变换分析频率点分布图 zM 1 当w0 1 螺线内缩当w00 逆时针转当