数字信号处理程佩青第三版课件_第四章_快速傅里叶变换(FFT)

第四章快速傅里叶变换 FFT 主要内容 DIT FFT算法DIF FFT算法IFFT算法Chirp FFT算法线性卷积的FFT算法 4 1引言 FFT FastFourierTransform1965年 Cooley Turky发表文章 机器计算傅里叶级数的一种算法 提出FFT算法 解决DFT运算量太大 在实际使用中受限制的问题 FFT的应用 频谱分析 滤波器实现 实时信号处理等 TI芯片TMS320系列 DSP芯片实现 TI公司的TMS320c30 10MHz时钟 基2 FFT1024点FFT时间15ms 典型应用 信号频谱计算 系统分析等 系统分析 频谱分析与功率谱计算 4 2直接计算DFT的问题及改进途径 1 DFT与IDFT 2 DFT与IDFT运算特点 同理 IDFT运算量与DFT相同 3 降低DFT运算量的考虑 FFT算法分类 时间抽选法DIT Decimation In Time频率抽选法DIF Decimation In Frequency 4 3按时间抽取 DIT 的FFT算法 DecimationInTime 1 算法原理设序列点数N 2L L为整数 若不满足 则补零 将序列x n 按n的奇偶分成两组 N为2的整数幂的FFT算法称基 2FFT算法 将N点DFT定义式分解为两个长度为N 2的DFT 记 1 再利用周期性求X k 的后半部分 将上式表达的运算用一个专用 蝶形 信流图表示 注 a 上支路为加法 下支路为减法 b 乘法运算的支路标箭头和系数 用 蝶形结 表示上面运算的分解 分解后的运算量 运算量减少了近一半 进一步分解 由于 仍为偶数 因此 两个点DFT又可同样进一步分解为4个点的DFT 蝶形 信流图表示 N点DFT分解为四个N 4点的DFT 类似的分解一直继续下去 直到分解为最后的两类蝶形运算为止 2点DFT 如上述N 8 23 N 4 2点中 类似进一步分解 进一步简化为蝶形流图 因此8点FFT时间抽取方法的信流图如下 FFT运算量与运算特点 1 N 2L时 共有L log2N级运算 每一级有N 2个蝶形结 2 每一级有N个数据中间数据 且每级只用到本级的转入中间数据 适合于迭代运算 3 计算量 每级N 2次复乘法 N次复加 每蝶形只乘一次 加减各一次 共有L N 2 N 2log2N次复乘法 复加法L N Nlog2N次与直接DFT定义式运算量相比 倍数 N2 Nlog2N 当N大时 此倍数很大 比较DFT 参考P150表4 1图4 6 可以直观看出 当点数N越大时 FFT的优点更突出 按时间抽取FFT蝶形运算特点 1 关于FFT运算的混序与顺序处理 位倒序处理 由于输入序列按时间序位的奇偶抽取 故输入序列是混序的 为此需要先进行混序处理 混序规律 x n 按n位置进行码位 二进制 倒置规律输入 而非自然排序 即得到混序排列 所以称为位倒序处理 位倒序实现 1 DSP实现采用位倒序寻址 2 通用计算机实现可以有两个方法 一是严格按照位倒序含义进行 二是倒进位的加N 2 倒位序 例计算 计算点FFT 用时间抽取输入倒序算法 问倒序前寄存器的数和倒序后的数据值 解 倒序前倒序倒序为倒序后 DITFFT中最主要的蝶形运算实现 1 参与蝶形运算的两类结点 信号 间 距离 码地址 与其所处的第几级蝶形有关 第m级的 结距离 为 即原位计算迭代 2 每级迭形结构为 蝶形运算两节点的第一个节点为k值 表示成L位二进制数 左移L m位 把右边空出的位置补零 结果为r的二进制数 3 的确定 第m级的r取值 四 FFT算法中一些概念 1 级 概念 将N点DFT先分成两个N 2点DFT 再是四个N 4点DFT 直至N 2个两点DFT 每分一次称为 一 级运算 因为N 2M所以N点DFT可分成M级如上图所示依次m 0 m 1 M 1共M级 2 组 概念 每一级都有N 2个蝶形单元 例如 N 8 则每级都有4个蝶形单元 每一级的N 2个蝶形单元可以分成若干组 每一组具有相同的结构 相同的因子分布 第m级的组数为 例 N 8 23 分3级 m 0级 分成四组 每组系数为m 1级 分成二组 每组系数为m 2级 分成一组 每组系数为 3 因子的分布 结论 每由后向前 m由M 1 0级 推进一级 则此系数为后级系数中偶数序号的那一半 DIT算法的其他形式流图 输入倒位序输出自然序输入自然序输出倒位序输入输出均自然序相同几何形状输入倒位序输出自然序输入自然序输出倒位序 参考P154 155 时间抽取 输入自然顺序 输出倒位序的FFT流图 例用FFT算法处理一幅N N点的二维图像 如用每秒可做10万次复数乘法的计算机 当N 1024时 问需要多少时间 不考虑加法运算时间 解当N 1024点时 FFT算法处理一幅二维图像所需复数乘法约为次 仅为直接计算DFT所需时间的10万分之一 即原需要3000小时 现在只需要2分钟 4 4按频率抽取 DIF 的FFT算法 与DIT FFT算法类似分解 但是抽取的是X k 即分解X k 成奇数与偶数序号的两个序列 设 N 2L L为整数 将X k 按k的奇偶分组前 先将输入x n 按n的顺序分成前后两半 DecimationInFrequency 一 算法原理 下面讨论 按k的奇偶将X k 分成两部分 显然 令 用蝶型结构图表示为 N 2仍为偶数 进一步分解 N 2 N 4 按照以上思路继续分解 即一个N 2的DFT分解成两个N 4点DFT 直到只计算2点的DFT 这就是DIF FFT算法 2个1点的DFT蝶形流图 进一步简化为蝶形流图 二 按频率抽取FFT蝶形运算特点 1 原位计算 L级蝶形运算 每级N 2个蝶形 每个蝶形结构 m表示第m级迭代 k j表示数据所在的行数 2 蝶形运算 对N 2L点FFT 输入自然序 输出倒位序 两节点距离 2L m N 2m 第m级运算 蝶形运算两节点的第一个节点为k值 表示成L位二进制数 左移m 1位 把右边空出的位置补零 结果为r的二进制数 存储单元 输入序列x n N个存储单元 系数 N 2个存储单元 三 DIT与DIF的异同 基本蝶形不同 DIT 先复乘后加减 DIF 先减后复乘 运算量相同 都可原位运算 DIT和DIF的基本蝶形互为转置 4 5IDFT的FFT算法 FFT应用一 一 从定义比较分析 与DFT的比较 1 旋转因子WN kn的不同 2 结果还要乘1 N 二 实现算法 直接使用FFT程序的算法 直接调用FFT子程序计算IFFT的方法 4 6线性调频Z变换 Chirp Z变换 算法 FFT应用二 单位圆与非单位圆采样 a 沿单位圆采样 b 沿AB弧采样 螺线采样 zk AW k k 0 1 M 1 Chirp Z变换的线性系统表示 由于系统的单位脉冲响应可以想象为频率随时间 n 呈线性增长的复指数序列 在雷达系统中 这种信号称为线性调频信号 ChirpSignal 因此 这里的变换称为线性调频Z变换 一 基本算法思路 4 7线性卷积的FFT算法 FFT应用三 若L点x n M点h n 则直接计算其线性卷积y n 需运算量 若系统满足线性相位 即 则需运算量 FFT法 以圆周卷积代替线性卷积 N 总运算量 次乘法 比较直接计算和FFT法计算的运算量 讨论 1 当 2 当 x n 长度很长时 将x n 分为L长的若干小的片段 L与M可比拟 1 重叠相加法 则 输出 其中 可以用圆周卷积计算 选 上面圆周卷积可用FFT计算 N 由于yi n 长度为N 而xi n 长度L 必有M 1点重叠 yi n 应相加才能构成最后y n 的 重叠相加法图形 和上面的讨论一样 用FFT法实现重叠相加法的步骤如下 计算N点FFT H k DFT h n 计算N点FFT Xi k DFT xi n 相乘 Yi k Xi k H k 计算N点IFFT yi n IDFT Yi k 将各段yi n 包括重叠部分 相加 重叠相加的名称是由于各输出段的重叠部分相加而得名的 例已知序列x n n 2 0 n 12 h n 1 2 1 试利用重叠相加法计算线性卷积 取L 5 y n 2 7 12 16 20 24 28 32 36 40 44 48 52 41 14 解 重叠相加法 x1 n 2 3 4 5 6 x2 n 7 8 9 10 11 x3 n 12 13 14 y1 n 2 7 12 16 20 17 6 y2 n 7 22 32 36 40 32 11 y3 n 12 37 52 41 14 2 重叠保存法 此方法与上述方法稍有不同 先将x n 分段 每段L N M 1个点 这是相同的 不同之处是 序列中补零处不补零 而在每一段的前边补上前一段保留下来的 M 1 个输入序列值 组成L M 1点序列xi n 如果L M 1 2m 则可在每段序列末端补零值点 补到长度为2m 这时如果用 DFT实现h n 和xi n 圆周卷积 则其每段圆周卷积结果的前 M 1 个点的值不等于线性卷积值 必须舍去 重叠保留法示意图 重叠保留法示意图 y k 2 7 12 16 20 24 28 32 36 40 44 48 52 41 14 x1 k 0 0 2 3 4 x2 k 3 4 5 6 7 x3 k 6 7 8 9 10 y1 k x1 k h k 11 4 2 7 12 x4 k 9