CHP4快速傅立叶变换ppt课件

第二部分傅立叶变换及其快速算法之第四章快速傅里叶变换 赵发勇zfy 72 物电学院 目录 4 1概述4 2时间抽取 DIT 基2FFT算法4 3频率抽取 DIF 基2FFT算法4 5分裂基算法4 6线性调频Z变换4 7与本章有关节MATLAB文件 4 1概述 快速傅里叶变换 FFT 是求解离散傅里叶变换 DFT 的快速算法 问题的提出直接计算N点DFT需要的计算量是多少 计算一个X k 需要N次复数乘法和N一1次复数加法 算出全部N点X k 共需N2次复数乘法和N N一1 次复数加法 总运算量近似地正比于N2 当N值很大 如2 D图像处理 运算量将非常庞大 同时 对于实时性很强的信号处理来说 必将对计算速度有十分苛刻的要求 为此 需要改进对DFT的计算方法 以减少总的运算次数 4 1概述 在正交矩阵中 虽然有N2个元素 但只有N个不同的值 且有些取值特别简单 主要由于旋转因子具有如下的特点 对称性 周期性 下面以四点DFT为例来说明快速算法的思路 4 1概述 4 1概述 交换矩阵第二列和第三列得 从上面的结果可以看出 利用对称性和周期性 求四点DFT只需要一次复数乘法 称为Coolkey Tukey算法 4 1概述 算法分类 N为2的整次幂 按基数分为基 2FFT算法 基 4FFT算法 混合基FFT算法 分裂基FFT算法 当N不是2的整次幂 典型的有Winograd算法 按抽取方法分 时间抽取 Decimation in Time 简称DIT 频率抽取 Decimation in Frequency 简称DIF 4 1概述 4 2时间抽取 DIT 基2FFT算法 为了将大点数的DFT分解为小点数的DFT运算 要求序列的长度N为N 2M M为正整数 该情况下的变换称为基2FFT 核心思想是 问题是如何分最有效 可以对时间变量分 DIT 也可对频率变量分 DIF 4 2时间抽取 DIT 基2FFT算法 基本思路 从时域将N点序列x n 按奇偶项分解为两组 分别计算两组N 2点DFT 然后再合成一个N点DFT 按此方法继续下去 直到2点DFT 从而减少运算量 算法具体步骤 一 算法的推导 1 序列x n 按奇偶项分解为两组 将DFT运算也相应分为两组 则 4 2时间抽取 DIT 基2FFT算法 2 两个N 2点的DFT合成一个N点DFT 问题 A k B k 都只有N 2个点 怎样得到X k 的后N 2点 利用周期性和对称性得 4 2时间抽取 DIT 基2FFT算法 4 2时间抽取 DIT 基2FFT算法 4 2时间抽取 DIT 基2FFT算法 3 继续分解 一直分解到两点DFT变换 A K 和B K 仍是高复合数 N 2 的DFT 我们可按上述方法继续以分解 令r 2l r 2l十1 l 0 1 N 4 1 则A K 和B K 可分别表示为 4 2时间抽取 DIT 基2FFT算法 4 2时间抽取 DIT 基2FFT算法 4 2时间抽取 DIT 基2FFT算法 令 则 4 2时间抽取 DIT 基2FFT算法 4 2时间抽取 DIT 基2FFT算法 同理 令 则 按此方法一直分解下去直到2点DFT 当N 8时 如下 4 2时间抽取 DIT 基2FFT算法 4 2时间抽取 DIT 基2FFT算法 4 2时间抽取 DIT 基2FFT算法 下面通过讨论寻找FFT的一般规律 二 算法的讨论 1 级 的概念在分解过程中 每分一次 称为一级运算 因为M log2N 所以N点DFT可以分解为M级 按抽取算法的信号流图中来定义 从左到右分别称为0级 1级到M 1级 4 2时间抽取 DIT 基2FFT算法 2 蝶形单元在算法的信号流图中 第m级存在这种运算 这种结构几何形状像蝴蝶 称为蝶形单元 p q是参于本蝶形单元运算的上 下节点的序号 由于第m级序号的两点只参于这一个蝶形单元的运算 其输出在第m十l级 且这一蝶形单元也不再涉及别的点 由于这一特点 在计算机编程时 我们可将蝶形单元的输出仍放在输入数组中 这一特点称为 同址运算 4 2时间抽取 DIT 基2FFT算法 4 2时间抽取 DIT 基2FFT算法 4 2时间抽取 DIT 基2FFT算法 由于一级都含有N 2个蝶形单元 每个蝶形单元需要1次复数乘法和两次复数加法 因此完成log2N级共需要的复数乘法和加法分别为 直接计算DFT时所需的复乘数与复加数都是与N2成正比的 所以采用FFT算法使运算量大大减少 显然 N值愈大 节省的运算量愈多 4 2时间抽取 DIT 基2FFT算法 3 组 的概念在分解过程中 每一级的N 2个蝶形单元可以分成若干组 每一组具有相同的结构和W因子分布 第m级可分成N 2m 1组 例 N 8 23 分3级 第一级的分组及Wr因子如下 m 0级 分成四组 因子为m 1级 分成二组 因子为m 2级 分成一组 因子为 4 2时间抽取 DIT 基2FFT算法 4 Wr因子的分布由上分析可知 结论 每由后向前 m由M 1 0级 推进一级 则此系数为后级系数中偶数序号的那一半 4 2时间抽取 DIT 基2FFT算法 5 码位倒置在FFT算法中 输出的频谱依照正常次序排列 但输入的序列x n 是按奇偶分开的 分开的规律 以N 8为例 按如下方法进行排序 1 将x n 的序号写成二进制x 000 x 001 x 110 x 111 2 将二进制的码进行翻转 得x 000 x 100 x 011 x 111 3 将二进制的翻转码转换为对应的十进制x 0 x 4 x 3 x 7 这就是按奇偶抽取得到的顺序 4 2时间抽取 DIT 基2FFT算法 说明 在上述的基2FFT算法中 由于每一步分解都是按输入序列x n 在时域上的次序是属于偶数还是奇数来抽取的 所以称为 按时间抽取法 或 时间抽取 上述的基2FFT算法中 抽取也可在频域进行 引出频率抽取 DIF 基2FFT算法 4 3频率抽取 DIF 基2FFT算法 设输入序列长度为N 2M M为正整数 频率抽取法将输入序列不是按奇 偶分组 而是按前后对半分开 这样可将N点DFT写成前后两部分 将该序列的频域的输出序列X k 也是N点序列 按其频域顺序的奇偶分解为越来越短的子序列 称为基2按频率抽取的FFT算法 也称为Sander Tukey算法 4 3频率抽取 DIF 基2FFT算法 算法分析现将输入x n 按n的顺序分前后两部分 前半子序列x n 0 n N 2 1 后半子序列x n N 2 0 n N 2 1 例 N 8时 前半序列为 x 0 x 1 x 2 x 3 后半序列为 x 4 x 5 x 6 x 7 考虑N点的DFT 由DFT定义得 4 3频率抽取 DIF 基2FFT算法 4 3频率抽取 DIF 基2FFT算法 按k的奇偶将X k 分成奇偶两部分 k 2r和k 2r 1 考虑k为偶数情况 令 4 3频率抽取 DIF 基2FFT算法 考虑k为奇数情况 令 4 3频率抽取 DIF 基2FFT算法 结论一个N点的DFT被分解为两个N 2点 与时间抽取法的推演过程一样 由于N 2M 因此 N 2仍为偶数 所以可以将N 2点DFT的输出X k 再分为偶数组和奇数组 这样就将一个N 2点的DFT分成两个N 4点DFT的输入 也是将N 2点的DFT的输入上 下对半分后通过蝶形运算而形成 直至最后为2点DFT 8点DIF基2FFT算法流图 4 3频率抽取 DIF 基2FFT算法 4 3频率抽取 DIF 基2FFT算法 4 3频率抽取 DIF 基2FFT算法 DIT与DIF的相同之处 1 DIF与DIT两种算法均为原位运算 2 DIF与DIT运算量相同 DIT与DIF的不同之处 1 DIF与DIT两种算法结构倒过来 DIF为输入顺序 输出乱序 运算完毕再运行 二进制倒读 程序 DIT为输入乱序 输出顺序 先运行 二进制倒读 程序 再进行求DFT 2 DIF与DIT根本区别 在于蝶形结不同 DIT的复数相乘出现在减法之前 DIF的复数相乘出现在减法之后 4 5分裂基算法 自从基2快速算法出现以来 人们仍在不断寻求更快的算法 1984年 法国的杜梅尔 P Dohamel 和霍尔曼 H Hollmann 将基2分解和基4分解糅合在一起 提出了分裂基FFT算法 其运算量比前几种算法都有所减少 运算流图却与基2FFT很接近 运算程序也很短 它是目前一种实用的高效新算法 4 5分裂基算法 分裂基算法又称基2 4算法 算法的基本思路是 对偶号输出使用基2算法 对奇号输出使用基4算法 下面首先介绍基4算法 令N 4M 对N点DFT可按下面方法进行频率抽取 分别令k 4r k 4r 2 k 4r 1 k 4r 3 其中 r 0 1 2 N 4 1 有 4 5分裂基算法 4 5分裂基算法 4 5分裂基算法 4 5分裂基算法 算法分析对于N 4M点DFT 当输出项k为偶数时 采用基2算法 即 当输出项k为奇数时 采用基4算法 即 4 5分裂基算法 4 5分裂基算法 4 5分裂基算法 分析 一个N点DFT可以分解为一个N 2点DFT和两个N 4点DFT 由x n x n N 4 x n N 2 和x n 3N 4 求N 2点DFT和N 4的DFT只需要两次乘法 可以减少运算量 N 2点DFT可进一步分解为一个N 4点DFT和两个N 8的DFT N 4的点DFT进一步分解为一个N 8点DFT和两个N 16的DFT 这样一直下去 直到分解为两点或4点DFT为止 4 5分裂基算法 结论 分裂基FFT算法结构同基2FFT算法结构相似 适用于N 2M的场合 并由M级运算实现 运算流图输入为顺序 输出为倒序 分裂基FFT算法的计算量 以上提出FFT算法 可以很快地求出全部DFT值 即求出有限长序列x n 的z变换X z 在单位园上N个等间隔抽样点zk处的抽样值 它要求N为高度复合数 即N可以分解成一些因子的乘积 例N 2L实际上 1 也许对其它围线上z变换取样发生兴趣 如语音处理中 常常需要知道某一围线z变换的极点所处的复频率 2 只需要计算单位圆上某一段的频谱 即M不等于N 如窄带信号 希望在窄带频率内频率抽样能够非常密集 提高分辨率 带外则不考虑 3 若N是大素数时 不能加以分解 又如何有效计算这种序列DFT 例N 311 若用基2则须补N 28 512点 要补211个零点 4 6线性调频Z变换 4 6线性调频Z变换 问题提出为了提高DFT的灵活性 须用新的方法 线性调频z变换 CZT 就是适用这种更为一般情况下 由x n 求X zk 的快速变换 CZT来自于雷达专业的专用词汇 Z变换采用螺线抽样 可计算单位圆上任一段曲线的Z变换 适用于更一般情况下 M不等于N 由x n 求X zr 的快速算法 达到频域细化的目的 这种变换称为线性调频Z变换 简称CZT 为适应z可以沿平面内更一般的路径取值 我们沿z平面上的一段螺线作等分角的抽样 则z的取样点Zr可表示为 已知N点序列x n 0 n N 1 其z变换为 其中M 表示欲分析的复频谱的点数 M不一定等于N A W都为任意复数 令 4 6线性调频Z变换 一 CZT的定义 4 6线性调频Z变换 上式即为CZT的定义 现在讨论A0 W0 0 0的含义 为输出M点的变换域值 r 0时的A0ej 0是CZT的起点 随着r的变化 r0 r1 RM 1构成CZT的变化路径 对于M 1点其极坐标为 4 6线性调频Z变换 4 6线性调频Z变换 CZT在现Z平面上的变换路径是一条螺旋线 1 A为起始样点位置 起点半径 大于1时 表示螺旋线在单位圆外 反之 在单位圆内 起点半相角 2 当W0 1 螺旋线内旋 反之外旋 3 当A0 W0 1时 CZT的变换路径为单位圆上的一段弧 起点为P 终点为Q 且M不一定等于N 4 6线性调频Z变换 4 6线性调频Z变换 考虑A0 W0 1时 在单位圆上CZT 且M不一定等于N 4 6线性调频Z变换 4 6线性调频Z变换 CZT的线性滤波计算步骤 4 6线性调频Z变换 二 CZT的计算方法分析 从上面的推导过程可以看出 计算CZT关键是计算一个线性卷积 其中 g n 应为N点序列 h n 应为偶对称的无限长序列 考虑到输出M点序列 h n 的实际长度应为M点 因此 可用DFT来实现两者的卷积 具体步骤如下 4 6线性调频Z变换 1 计算并设置g n 4 6线性调频Z变换 2 计算并设置h n 将h n 设置成长度为L的序列 考虑到其为偶对称序列 且取M点 输出M点序列 取 如下图所示 4 6线性调频Z变换 4 6线性调频Z变换 3 计算h n 和g n 的DFT 得到L点序列H k 和G k 4 令Y k H k G k 乘积 后作Y k 的IDFT 反变换 得到时域输出序列y r 5 取y r 的前M点 并乘以W r2 2 则得最后的输出X zr 即 与标准FFT算法相比 CZT算法有以下特点 1 输入序列长度N及输出序列长度M不需要相等 且N及M不必是高度合成数 二者均可为素数 2 Zk的角间隔是任意的 说明其频率分辨率也是任意可控的 角间隔小 分辨率高 反之 分辨率低 3 周线不必是z平面上的圆 在语音分析中螺旋周线具有某些优点 4 由于起始点z0