第5章快速傅里叶变换

第4章快速傅里叶变换 FFT 本章将介绍快速傅立叶变换 FFT 的具体实现方法 本章内容包括 直接计算DFT的问题及改进的途径 按时间抽取 DIT 的FFT算法 库利 图基算法 频率抽取 DIF 的FFT算法 桑德 图基算法 离散傅里叶反变换 IDFT 的快速计算方法 混合基算法 分裂基FFT算法 2020 4 7 2 目录 4 1直接计算DFT的问题及改进途径4 2按时间抽取的FFT算法4 3按频率抽取的FFT算法4 4离散傅里叶反变换的快速算法4 5N为复合数的FFT算法 混合基算法4 6分裂基FFT算法 2020 4 7 3 DFT是信号分析与处理中的一种重要变换 因直接计算DFT的计算量与变换区间长度N的平方成正比 当N较大时 计算量太大 所以在快速傅里叶变换 简称FFT 出现以前 直接用DFT算法进行谱分析和信号的实时处理是不切实际的 直到1965年发现了DFT的一种快速算法以后 情况才发生了根本的变化 2020 4 7 4 4 1直接计算DFT的问题及改进的途径 设x n 为N点有限长序列 其DFT为 k 0 1 2 3 4 N 1 4 1 反变换 IDFT 为 n 0 1 2 3 4 N 1 4 2 二者的差别只在于WN的指数符号不同 以及差一个常数乘因子1 N 因而下面只讨论DFT正变换 4 1 式的运算量 4 2 式运算量是完全相同的 2020 4 7 5 一般来说 x n 和都是复数 X k 也是复数 因此每计算一个X k 值 需要N次复数乘法 x n 与相乘 以及 N 1 次复数加法 而X k 一共有N个点 k从0取到N 1 所以完成整个DFT运算总共需要N2次复数乘法及N N 1 次复数加法 已知复数运算实际上是由实数运算来完成的 4 1 式可写成 4 3 2020 4 7 6 式 4 3 运算量计算 由上式可见 一次复数乘法需用四次实数乘法和二次实数加法 一次复数加法则需二次实数加法 因而每运算一个X k 需要4N次复数乘法及2N 2 N 1 2 2N 1 次实数加法 所以整个DFT运算共需要4N2次实数乘法和N 2 2N 1 2N 2N 1 次实数加法 因而 直接计算DFT 乘法次数和加法次数都是个N2成正比的 但N很大时 运算量是很可观的 例如 当N 8时 DFT需64次复乘 而当N 1024时 DFT所需复乘为1 048 576次 既一百多万次复乘运算 这对实时性很强的信号处理来说 对计算速度要求太高了 因而需要改进DFT的计算方法 以大大减少运算次数 2020 4 7 7 2 减少运算工作量的途径 仔细观察DFT的运算就可看出 利用系数的以下固有特性 就可减少DFT的运算量 1 的对称性 2 的周期性 由此可得出 2020 4 7 8 利用这些特性 使DFT运算中有些项可以合并 利用的对称性和周期性 可以将长序列的DFT分解为短序列的DFT 而前面已经说到 DFT的运算量是与N2成正比的 所以N越小越有利 因而小点数序列的DFT比大点数序列的DFT的运算量要小 快速傅立叶变换 FFT 算法正是基于这样的基本思路而发展起来的 它的算法基本上可以分成两大类 即按时间抽取 Decimation In Time 简称DIT 法和按频率抽取 Decimation In Frequency 简称DIF 法 2020 4 7 9 4 2按时间抽取的FFT算法 4 2 1算法原理先设序列长度为N 2L L为整数 如果不满足这个条件 可以人为地加上若干零值点 使得达到这一要求 这种N为2的整数幂的FFT 也称基 2FFT 将N 2L的序列x n n 0 1 N 1 先按n的奇偶分成两组 r 0 1 2 1 4 4 则可将DFT化为 2020 4 7 10 由于 则上式可表示成 4 5 2020 4 7 11 式中X1 k 及X2 k 分别是x1 r 及x2 r 的N 2点DFT 由 4 5 式看出 一个N点DFT已分解两个N 2点的DFT 它们按 4 5 式又组合成一新的DFT 但是 x1 r x2 r 以及X1 k X2 k 都是N 2点的序列 即r k满足r k 0 1 2 N 2 1 而X k 却有N点 而用 4 5 式计算得到的只是X k 的前一半项数的结果 要用X1 k X2 k 来表达全部X k 值 还必须应用系数的周期性 即 4 6 4 7 2020 4 7 12 这样可得到 同理可得 4 8 4 9 4 8 式 4 9 式说明了后半部分k值 N 2 k N 1 所对应的X1 k X2 k 分别等于前半部分k值 0 k N 2 1 所对应的X1 k X2 k 2020 4 7 13 再考虑到的对称性 把 4 8 式 4 9 式 4 10 式代入到 4 5 式 就可将X k 表达为前后两部分前半部分X k k 0 1 1 4 11 后半部分X k k 0 1 1 4 12 4 10 2020 4 7 14 蝶形信号流图 这样 只要求出0 N 2 1 区间的所有X1 k 和X2 k 值 即可求出0 N 1 区间内的所有X k 值 这就大大节省了运算 图4 1的蝶形信号流图符号表示式 4 11 和 4 12 支路上没有标出系数时 则该支路的传输系数为1 图4 1时间抽取法蝶形运算流图符号 2020 4 7 15 图4 2表示N 23 8分解过程 其中输出值X 0 X 3 由公式 4 11 给出的 而输出值X 4 X 7 由公式 4 12 给出 图4 2按时间抽取 将N点DFT分解为两个N 2点DFT 2020 4 7 16 与第一次分解相同 将x1 r 按奇偶分解成两个N 4长的子序列x3 l 和x4 l 即 那么 X1 k 又可表示为 4 13 2020 4 7 17 式中 且 k 0 1 1 4 14 4 15 图4 3给出N 8时 将一个N 2点DFT分解成两个N 4点DFT 由此两个N 4点DFT组合成一个N 2点DFT的流图 2020 4 7 18 图4 3由两个N 4点DFT组合成一个N 2点DFT 2020 4 7 19 X2 k 也可进行同样的分解 式中 k 0 1 1 k 0 1 1 4 16 4 17 2020 4 7 20 N 8点DFT 将系数统一为 则一个N 8点DFT就可以分解为4个N 4 2点的DFT 这样可得如图4 4的流图 根据上面同样的分析知道利用四个N 4点的DFT及两级蝶形组合运算来计算N点DFT 比只用一次分解蝶形组合方式的计算量又减少了大约一半 完整的8点DFT流图 如图4 5所示 这种方法的每一步分解都是按输入序列在时间上的次序是属于偶数还是属于奇数来分解为两个更短的子序列 所以称为 按时间抽取法 DIT 2020 4 7 21 图4 4按时间抽取 将一个N点DFT分解为四个N 4点DFT N 8 2020 4 7 22 图4 5N 8按时间抽取法FFT运算图 这种方法的每一步分解都是按输入序列在时间上的次序是属于偶数还是属于奇数来分解为两个更短的子序列 所以称为 按时间抽取法 DIT FFT 2020 4 7 23 4 2 2运算量 由按时间抽取法FFT的流图可见 当N 2L共有L级蝶形 每级都由N 2个蝶形运算组成 每个蝶形有一次复乘 二次复加 因而每级运算都需N 2次复乘和N次复加 这样L级运算总共需要复乘数 4 18 4 19 复加数 直接计算DFT与FFT算法的计算量之比为 4 20 2020 4 7 24 在图4 6中 表示出直接计算DFT与FFT算法所需运算量分别与点数N的关系曲线 可以更加直观地看出FFT算法的优越性 尤其是当点数N越大时 FFT的优点就更突出 图4 6FFT算法与直接计算DFT所需乘法次数的比较曲线 2020 4 7 25 4 2 3按时间抽取的FFT算法的特点 为求得任何N 2L点的时间抽取基 2FFT信号流图或设计出硬件实现电路 下面来考察这种按时间抽取法在运算方式上的特点 1 原位运算原位运算从图4 5可以看出这种运算是很有规律的 其每级 每列 计算都是由N 2个蝶形运算构成 每一个蝶形结构完成下述基本迭代运算 4 21 式中m表示第m列迭代 k j为数据所在行数 式 4 21 的蝶形运算结构由一次复乘和两次复加 减 组成 如图4 7所示 2020 4 7 26 图4 7时间抽取法蝶形运算结构 由图4 5看出 某一列的任何两个节点k和j的节点变量进行蝶形运算后 得到结果为下一列k j两节点的节点变量 而与其他节点变量无关 因而可以采用原位运算 这样存储数据只需N个存储单元 下一级的运算仍采用这种原位方式 只不过进入蝶形结的组合关系有所不同 这种原位运算结构可以节省存储单元 降低设备成本 2020 4 7 27 2 倒位序规律 由图4 5看出 按原位计算时 FFT的输出X k 是按正常顺序排列在存储单元中 即按X 0 X 1 X 7 的顺序排列 但是这时输入x n 却不是按自然顺序存储的 而是按x 0 x 4 x 7 的顺序存入存储单元 看起来好象是 混乱无序 的 实际上是有规律的 我们称之为倒位序 造成倒位序的原因是输入x n 按标号n的偶奇的不断分组而造成 如果n用二进制数表示为 n2n1n0 2 当N 8时 二进制为三位 第一次分组 由图4 2看出 n为偶数在上半部分 n为奇数在下半部分 这可以观察n的二进制数的最低位n1 n0 0则序列值对应于偶数抽样 n0 1 则序列值对应于奇数抽样 2020 4 7 28 图4 8描述倒位序的树状图 下一次则根据次最低位n1的0 1来分偶奇 而不管原来的子序列是偶序列还是奇序列 这种不断分成偶数子序列和齐数子序列的过程可用图4 8的二进制树状图来描述 这就是DIT的FFT算法输入序列的序数成为倒位序的原因 2020 4 7 29 3 倒位序的实现 一般实际运算中 总是先按自然顺序将输入序列存入存储单元 为了得到倒位序的排列 可以通过变址运算来完成 如果输入序列的序号n用二进制数 例如n2n1n0 表示 则其倒位序二进制数 就是 n2n1n0 这样 在原来自然顺序时应该放x n 的单元 现在倒位序后应存放x 例如N 8时 x 3 的标号是n 3 它的二进制数是011 倒位序的二进制数是110 即 6 所以存放x 011 的单元现在应该存放入x 110 表4 2列出了N 8时的自然顺序二进制数以及相应的倒位序二进制数 2020 4 7 30 表4 2码位的倒位序 N 8 2020 4 7 31 FFT原位运算所要求的倒位序的变址功能如图4 9所示 当n 时 不必调换 当n 时 必须将原来存放数据x n 的存信单元内调入数据x 而将存放x 的存储单元内调入x n 为了避免把已调换过的数据再次调换 保证只调换一次 否则又回到原状 只需看 是否比n小 若 比n小 则意味着此x n 在前边已和x 互相调换过 不必再调换了 只有当 n是地 才将原存放的x n 及存放x 的存储单元内的内容互换 这样就得到输入所需的倒位序的顺序 可以看出 其结果与图4 5的要求是一致的 2020 4 7 32 图4 9倒相序的变址处理 4 蝶形运算两结点的 距离 以图4 5的8点FFT为例 其输入是倒位序的 输出是自然顺序的 其第一级 第一列 每个蝶形的两节点 距离 为1 第二级 每一蝶形的两节点 距离 为2 第三级 每个蝶形的两节点 距离 为4 由此类推得 对N 2L点FFT 当输入为倒位序 输出为正常顺序时 其第m级运算 每个蝶形的两节点间 距离 为2m 1 2020 4 7 33 5 的确定 由于对第m级运算 一个DIT蝶形运算的两节点间 距离 为为2m 1 因而 4 21 式可写成 4 22 r的求解方法为 把 4 22 式中 蝶形运算两节点中的第一个节点标号值 即k值 表示成L位 注意为2L 二进制数 把此二进制数乘上为2L m 即将此L位二进制数左移L m位 注意m是第m级的运算 把右边空出的位置补零 此数即为所求r的二进制数 2020 4 7 34 6 存储单元 由于是原位运算 只需有输入序列x n n 0 1 2 N 1 的N个存储单元 加上系数 r 0 1 N 2 1 的个存储单元 从图4 5看出 因子最后一列有N 2种 顺序为 其余可类推 2020 4 7 35 4 2 4按时间抽取的FFT算法的其他形式流图 显然 对于任何流图 只要保持各节点所连的支路及其传输系数不变 则不论节点位置怎么排列所得流图总是等效的 所得最后结果都是x n 的DFT的正确结果 只是数据的提取和存放的次序不同而已 这样就可得到按时间抽取的FFT算法的若干其他形式流图 将图4 5中跟x 4 水平相连的所有节点与跟x 1 水平相连的节点位置对调 再将跟x 6 水平相连的所有节点与跟x 3 水平相连的所有节点对调 其余诸节点保持不变 则可得图4 10的流图 2020 4 7 36 图4 10时间抽取 输入自然顺序 输出倒位序的FFT流图 2020 4 7 37 经过简单变换 可得图4 11的流图 它的输入与输出都是按自然顺序排列的 不需要倒位序重排数据 但是它却不能进行原位运算 这样N个输入数据 至少要2N个复数存储单元 可以在专用硬件实现中利用这种流程 2020 4 7 38 图4 11时间抽取 输入输出皆为自然顺序的FFT流图 2020 4 7 39 当没有随机存储器时 采用图4 12的流图特别有用 此流图输入是倒位序的 而输出则是自然顺序的 流图各级的几何形状完全一样 只是级与级之间的支路传输比是改变的 这就有可能按顺序存取数据 图4 13的流图也是各级几何形状相同的流图 不过它的输入是自然顺序排列而输出是倒位序排列的 2020 4 7 40 图4 12时间抽取 各级具有相同几何形状 输入倒位序 输出自然顺序的FFT流图 2020 4 7 41 图4 13时间抽取 各级具有相同几何形状 输入倒位序 输出倒位序的FFT流图 2020 4 7 42 4 3按频率抽取的FFT算法 本节讨论按频率抽取 DIF 的FFT算法 它是把输出序列X k 也是N点序列 按其顺序的偶 奇来分解为越来越短的序列 4 3 1算法原理仍设序列长度为N 2L L为整数 在把输出X k 按k的偶奇分组之前 先把输入按n的顺序分成前后两半 注意 这不是频率抽取 2020 4 7 43 由于 故 式中用的是 而不是 因而这并不是N 2点DFT 则 k 0 1 N 1 k 0 1 N 1 4 23 2020 4 7 44 当k为偶数时 1 k 1 k为奇数时 1 k 1 为此按k的奇偶可将X k 分为两部分 令 r 0 1 1 则 4 24 4 25 2020 4 7 45 式 4 24 为前一半与后一半输入之和的N 2点DFT 4 25 为前一半输入与后一半输入之差再与之积的N 2点DFT 令 n 0 1 1 4 26 4 27 r 0 1 1 2020 4 7 46 图4 14频率抽取法碟形运算流图符号 式 4 26 所表示的运算关系可以用图4 14所示的蝶形运算来表示 这样就把一个N点DFT按频率k的奇偶分解为两个N 2点的DFT了 如式 4 27 所示 N 8时 上述分解过程示于图4 15 2020 4 7 47 图4 15按频率抽取 将N点DFT分解为两个N 2点DFT的组合 N 8 2020 4 7 48 与时间抽取法的推导过程一样 由于N 2L N 2仍是一个偶数 因而可以将每个N 2点DFT的输出再分解为偶数组与奇数组 这就将N 2点DFT进一步分解为2个N 4点DFT 这两个N 4点DFT的输入是将N 2点DFT的输入上下对半分开后通过蝶形运算而形成 图4 16示出了这一步分解的过程 2020 4 7 49 图4 16按频率抽取一个N点DFT分解4个N 4点 N 8 2020 4 7 50 这样的分解可以一直进行到第L次 N 2L 第L次实际上是做2点DFT 它只有加减运算 但是为了比较并为了统一运算结构 仍然采用系数为的蝶形运算来表示 这N 2个2点DFT的N个输出就是x n 的N点的DFT的结果X k 图4 17 a 表示一个N 8的完整的按频率抽取的FFT结构 从图4 17 a 的流图看出 频率抽取法FFT的这一流图仍是原位运算的 这里就不再讨论了 2020 4 7 51 图4 17 a 按频率抽取法FFT流图 N 8 2020 4 7 52 图4 17 b 频率抽取法蝶形运算结构 4 3 4频率抽取法与时间抽取法的异同由图4 17 a 与图4 5相比较看出 DIF法与DIT法的区别是 图4 17 a 的DIF输入是自然顺序 输出是倒位序的 这与图4 5的DIT法正好相反 但这不是实质性的区别 因为DIF法与DIT法一样 都可将输入或输出进行重排 使二者的输入输出顺序变成自然顺序 2020 4 7 53 由于输入输出是按自然数顺序还是按倒位序排列的不同 有四种原位运算的FFT流图 时间抽取法两种 频率抽取法两种 这里把它们的蝶形结对偶节点距离 第m级运算中蝶形结的计算公式以及中r的求法 列在表4 3中 供读者参考 2020 4 7 54 4 4离散傅里叶反变换的快速计算方法 前面讨论的FFT算法 同样可以适用于离散傅里叶反变换 IDFT 运算 并可简称为IFFT 即快速傅里叶反变换 由IDFT公式 4 29 并与DFT公式 4 30 比较得 若把DFT运算中的每一个系数换成 并且最后再乘以常数1 N 则以上所有时间抽取或频率抽取的FFT算法都可以拿来运算IDFT 2020 4 7 55 图4 18时间抽取IFFT流图 N 8 这样可得到图4 18的IDFT流图 2020 4 7 56 当把频率抽取的FFT流图用于IDFT时 由于输出变量变成x n 此时是按x n 的奇偶分组了 因而应改称时间抽取IDFT流图 同样 如果将时间抽取FFT流图用于IDFT时 由于输入变成X k 则是按输入X k 的奇偶分组了 因而应改称为频率抽取IDFT流图 以上讨论的IDFT算法 排程序虽很方便 但是要稍稍改动FFT的程序和参数才能实现 下面讨论一种完全不用改变FFT的程序就能计算IFFT的方法 将IDFT公式 4 29 取共轭 2020 4 7 57 因而 这说明 只要先将X k 取共轭 就可以直接利用FFT子程序 最后再将运算结果取一次共轭 并乘以1 N 即得到x n 值 这样FFT和IFFT运算都可共用一个子程序块 就很方便了 4 31 2020 4 7 58 4 5N为复合数的FFT算法 混合基算法 以上讨论的是序列的点数为N 2L 即以2为基数的FFT算法 这种算法在实际中使用得最多 因为它的程序简单 效率高 使用方便 若不满足为N 2L 则可有以下几种办法 1 将x n 补一些零点值的办法 以使N增长到最临近的一个2L数值 由DFT的性质知道 有限长序列补零之后 并不影响其频谱X ej 只不过其频谱的抽样点数增加了 所造成的结果是增加了计算量而已 但是 有时计算量增加太多 造成很大的浪费 例如x n 的点数N 300 则须补到N 28 512 要补212个零值点 因而研究N 2L时的FFT算法十分必要 2020 4 7 59 2 如果要求准确的N点DFT 而N又是素数 则只能采用直接DFT方法 或者用CZT Chirpz变换 方法 3 若N是一个复合数 即它可以分解为一些因子的乘积 则可以用FFT的一般算法 即混合基FFT算法 而基2算法只是这种一般算法的特例 4 5 1整数的多基多进制表示形式 1 对于二进制 N 2L 则任一个n 2L的正整数n可以用2为基数表示成二进制形式 nL 1 nL 2 n1 n0 其中ni为0或1 i 0 1 L 1 这个二进制数所表示的数值为 2020 4 7 60 将此二进制数倒位序后成为 n 2 nL 1 nL 2 n1 n0 2 其所代表的数值为 2 对于r进制 多进制 设N rL r和L皆为大于1的整数 r 2时 即为二进制 则任一个n 2L的正整数n 可以用r为基数表示成r进制形式 nL 1 nL 2 n1 n0 其中 ni 0 1 r 1 i 0 1 L 1这个r进制数所表示的数值为 4 32 4 33 4 34 2020 4 7 61 将此r进制数倒位序后成为 n nL 1 nL 2 n1 n0 r 其所代表的数值为 3 对于多基多进制或称混合基 这是最一般的情况 可以包括上面两种单基的情况 此时N可表示成复合数N r1 r2 rL 则对于n r1 r2 rL的任何一个正整数n 可以按L个基r1 r2 rL表示为多基多进制形式 nL 1 nL 2 n1 n0 r1 r2 rL 这一多基多进制数所代表的数值为 4 35 4 36 4 37 其倒位序形式为 它所代表的数值为 2020 4 7 62 4 5 2N r1r2的快速算法 1 算法原理要计算的N点DFT为 k 0 1 N 1 4 39 设N是一个复合数N r1r2 按上面讨论 可将n n N 的数用下面的公式表达 4 40 同样 若令N r2r1则可将频率变量k k N 用下式表达 4 41 2020 4 7 63 将式 4 40 与 4 41 代入式 4 39 可得 上面推导中应用了的结果 这里是用n1和n0表示的 所以要对n1和n0的所有位求和 所以n的单求和号变成了n1和n0的两个求和号 4 42 2020 4 7 64 4 42 式可进一步表示为 注意 表达式中x n1 n0 表示x n 的n为r2进制顺序排列 X k0 k1 表示X k 的k为r1进制顺序排列 而x n1 n0 则应表示x n 的n为r2进制倒位顺序排列 X k0 k1 则应表示X k 的k为r1进制倒位序排列 式中 4 43 2020 4 7 65 4 44 4 45 4 46 4 47 2020 4 7 66 式 4 44 表示n0为参量时 n0 0 1 r2 1 输入变量n1与输出变量k0之间n1点DFT 总共有r2个 n0的数目 r1点DFT 而X1 k0 n0 的序列值为r1 r1 N 式 4 45 表示式 4 44 的X1 k0 n0 乘以因子所组成的新序列 成为旋转因子 式 4 46 表示k0为参量时 k0 0 1 r1 1 输入变量n0与输出变量k1之间的点DFT 总共有r1个 n0的数目 r2点DFT 而X2 k0 k1 的序列值为r2 r1 N 同时可以看出 X2 k0 k1 中的变量是按进位制位序排列的 式 4 47 则表示 最后要利用k r1k1 k0进行整序 以恢复出X k1 k0 X k 2020 4 7 67 因而可将N为复合数N r1r2的DFT算法的步骤归纳如下 1 将x n 改写为x n1 n0 利用 2 利用 4 44 式 把r2个r1点DFT 得X k1 k0 3 利用 4 45 式 把N个X1 k0 n0 乘以相应的旋转因子 组成 4 利用 4 46 式 做r1个r2点DFT 得X2 k0 k1 5 利用 4 47 式 进行整序 得到X k1 k0 X k 其中k r1k1 k0 2020 4 7 68 4 6分裂基FFT算法 从基 2按时间和按频率抽选的推导过程中看出 每级抽选时 每一组的偶序号部分 时间抽选看输入序号 频率抽选看输出序号 都不乘