第2章：离散傅里叶变换及其快速算法-2

2 3快速傅里叶变换 FFT 快速傅里叶变换 FFT 是计算DFT的一种快速有效方法 从前面的讨论中看到 有限长序列在数字技术中占有很重要的地位 有限长序列的一个重要特点是其频域也可以离散化 即离散傅里叶变换 DFT 虽然频谱分析和DFT运算很重要 但在很长一段时间里 由于DFT运算复杂 并没有得到真正的运用 而频谱分析仍大多采用模拟信号滤波的方法解决 直到1965年首次提出DFT运算的一种快速算法以后 情况才发生了根本变化 人们开始认识到DFT运算的一些内在规律 从而很快地发展和完善了一套高速有效的运算方法 快速付里变换 FFT 算法 FFT的出现 使DFT的运算大大简化 运算时间缩短一 二个数量级 使DFT的运算在实际中得到广泛应用 1 DFT运算的特点 首先分析有限长序列x n 进行一次DFT运算所需的运算量 一般 x n 和wnkN都是复数 因此 每计算一个X k 值 要进行N次复数相乘 和N 1次复数相加 X k 一共有N个点 故完成全部DFT运算 需要N2次复数相乘和N N 1 次复数相加 在这些运算中 乘法比加法复杂 需要的运算时间多 尤其是复数相乘 每个复数相乘包括4个实数相乘和2个实数相加 例又每个复数相加包括2个实数相加 所以 每计算一个X k 要进行4N次实数相乘和2N 2 N 1 2 2N 1 次实数相加 因此 整个DFT运算需要4N2实数相乘和N 4N 2 次实数相加 从上面的分析看到 在DFT计算中 不论是乘法和加法 运算量均与N2成正比 因此 N较大时 运算量十分可观 例 计算N 10点的DFT 需要100次复数相乘 而N 1024点时 需要1048576 一百多万 次复数乘法 如果要求实时处理 则要求有很高的计算速度才能完成上述计算量 反变换IDFT与DFT的运算结构相同 只是多乘一个常数1 N 所以二者的计算量相同 FFT算法的基本思想 考察DFT与IDFT的运算发现 利用以下两个特性可减少运算量 1 系数是一个周期函数 它的周期性和对称性可用来改进运算 提高计算效率 对称性 又如因此周期性 利用这些周期性和对称性 使DFT运算中有些项可合并 2 利用的周期性和对称性 把长度为N点的大点数的DFT运算依次分解为若干个小点数的DFT 因为DFT的计算量正比于N2 N小 计算量也就小 FFT算法正是基于这样的基本思想发展起来的 它有多种形式 但基本上可分为两类 时间抽取法和频率抽取法 2 按时间抽取的FFT N点DFT运算的分解 先从一个特殊情况开始 假定N是2的整数次方 N 2M M 正整数首先将序列x n 分解为两组 一组为偶数项 一组为奇数项 r 0 1 N 2 1 将DFT运算也相应分为两组 因为故 注意到 H k G k 有N 2个点 即k 0 1 N 2 1 还必须应用系数wkN的周期性和对称性表示 k 的N 2 N 1点 由得 依此类推 G k 和H k 可以继续分下去 这种按时间抽取算法是在输入序列分成越来越小的子序列上执行DFT运算 最后再合成为 点的DFT 可见 一个N点的DFT被分解为两个N 2点的DFT 这两个N 2点的DFT再合成为一个N点DFT 蝶形信号流图 将G k 和H k 合成X k 运算可归结为 蝶形运算的简化 图 a 为实现这一运算的一般方法 它需要两次乘法 两次加减法 考虑到 bW和bW两个乘法仅相差一负号 可将图 a 简化成图2 7 b 此时仅需一次乘法 两次加减法 图 b 的运算结构像一蝴蝶通常称作蝶形运算结构简称蝶形结 采用这种表示法 就可以将以上所讨论的分解过程用流图表示 N 23 8的例子 按照这个办法 继续把N 2用 除 由于N 2M 仍然是偶数 可以被 整除 因此可以对两个N 2点的DFT再分别作进一步的分解 即对 G k 和 H k 的计算 又可以分别通过计算两个长度为N 4 2点的DFT 进一步节省计算量 见图 这样 一个 点的DFT就可以分解为四个 点的DFT N 4点 N 4点 N 4点 N 4点 由四个2点DFT组成8点DFT 最后剩下的是2点DFT 它可以用一个蝶形结表示 从而得到一个8点的完整的按时间抽取运算的流图 由于这种方法每一步分解都是按输入时间序列是属于偶数还是奇数来抽取的 所以称为 按时间抽取法 或 时间抽取法 按时间抽取的8点FFT 时间抽取法FFT的运算特点 1 蝶形运算 2 原位计算 4 序数重排 3 蝶形类型随迭代次数成倍增加 1 蝶形运算对于N 2M 总是可以通过M次分解最后成为2点的DFT运算 这样构成从x n 到X k 的M级运算过程 从上面的流图可看到 每一级运算都由N 2个蝶形运算构成 因此每一级运算都需要N 2次复乘和N次复加 这样 经过时间抽取后M级运算总共需要的运算 复乘复加而直接运算时则与N2成正比 例N 2048 N2 4194304 N 2 log2N 11264 N2 N 2 log2N 392 4 FFT显然要比直接法快得多 2 原位计算当数据输入到存储器中以后 每一级运算的结果仍然储存在同一组存储器中 直到最后输出 中间无需其它存储器 这叫原位计算 每一级运算均可在原位进行 这种原位运算结构可节省存储单元 降低设备成本 还可节省寻址的时间 3 蝶形类型随迭代次数成倍增加观察8点FFT的三次迭代运算 第一级迭代 有一种类型的蝶形运算系数W08 两个数据点间隔为1第二级迭代 有二种类型的蝶形运算系数W08 W28 参加运算的两个数据点间隔为2 第三级迭代 有四类蝶形运算系数W08 W18 W28 W38 参加运算的两个数据点间隔为4 结论 每迭代一次 蝶形类型增加一倍 数据点间隔也增大一倍 每一级的取数间隔和蝶形类型种类均为2i 1 i 1 2 M 4 序数重排对按时间抽取FFT的原位运算结构 当运算完毕时 正好顺序存放着X 0 X 1 X 2 X 7 因此可直接按顺序输出 但这种原位运算的输入x n 却不能按这种自然顺序存入存储单元中 而是按x 0 x 4 x 2 x 6 x 7 的顺序存入存储单元 这种顺序看起来相当杂乱 然而它也是有规律的 当用二进制表示这个顺序时 它正好是 码位倒置 的顺序 例如 原来的自然顺序应是x 1 的地方 现在放着x 4 用二进制码表示这一规律时 则是在x 001 处放着x 100 x 011 处放着x 110 表2 1自然顺序与二进制按位翻转顺序 在实际运算中 一般直接将输入数据x n 按码位倒置的顺序排好输入很不方便 总是先按自然顺序输入存储单元 然后再通过变址运算将自然顺序的存储转换成码位倒置顺序的存储 然后进行FFT的原位计算 目前有许多通用DSP芯片支持这种码位倒置的寻址功能 第一次分偶 奇 根据最低位n0的0 1状态来分 若n0 0 则为偶序列 n0 1则为奇序列 得到两组序列 000010100110 001011101111第二次对这两个偶 奇序列再分一次偶 奇序列 这就要根据n1的 状态 若n1 0 则为偶序列 n1 1则为奇序列 得到四组序列 000100 010110 001101 011111同理 再根据n2的 状态来分偶 奇序列 直到不能再分偶 奇时为止 对于N 8 n2已是最高位 最后一次分得结果为000 100 010 110 001 101 011 111 3 按频率抽取的FFT 按输出X k 在频域的顺序上属于偶数还是奇数分解为两组 对于N 2M情况下的另外一种普遍使用的FFT结构是频率抽取法 对于频率抽取法 输入序列不是按偶数奇数 而是按前后对半分开 这样便将N点DFT写成前后两部分 又把X k 进一步分解为偶数组和奇数组 令a n x n x n N 2 b n x n x n N 2 wnN这两个序列都是N 2点的序列 将其代入上两式 得这正是两个N 2点的DFT运算 即将一个N点的DFT分解为两个N 2点的DFT 上式的运算关系可用下图表示 1 以N 8的频率抽取为例 先蝶形运算再DFT 按频率抽取将8点DFT分解成四个2点DFT 与时间抽取法一样 由于N 2M N 2仍是一个偶数 因此可以将N 2点DFT再分解为偶部和奇部 这样就进一步将N 2 点DFT分解为两个N 4点的DFT 进一步分解 一个N 2M点的DFT通过M次分解后 最后只剩下全部是2点的DFT 2点DFT实际上只有加减运算 但为了比较 也为了统一运算的结构 仍用一个系数为W0N的蝶形运算来表示 下图是N 8的频率抽取法FFT流图 频率抽取法的流图是时域抽取法流图的左右翻转 N 8的按频率抽取FFT运算流图 频率抽取法FFT的运算特点 1 蝶形运算对于任何一个2的整数幂N 2M 总是可以通过M次分解最后完全成为2点的DFT运算 这样的M次分解 就构成从x n 到X k 的M级运算过程 将频率抽取法的流图反转 并将输入变输出 输出变输入 得到的便是时间抽取法流图 反映了时域 频域的对称法 频率抽取法也共有M级运算 N 2M 其运算量与时间抽取法相同 2 原位计算类似于时间抽取法 当数据输入到存储器中以后 每一级运算的结果仍然储存在同一组存储器中 直到最后输出 中间无需其它存储器 所以频域抽取法也可进行原位计算 3 序数重排它的输入正好是自然顺序 但它的输出却是码位倒置的顺序 因此运算完毕后 要通过变址运算将码位倒置的顺序转换为自然顺序 然后输出 变址方法同时间抽取法 4 蝶形类型随迭代次数成倍减少 与时间抽取法相反 第一级迭代中有N 2种蝶形运算系数 参加蝶形运算的两个数据相隔N 2 随后每次迭代 蝶形类形比前一级减少一倍 间距也减少一倍 最后一级迭代 蝶形类形只有一种W0N 数据间隔为1 由这几点规律可以看出 频率抽取法与时间抽到法是两种等价的FFT运算 4 N为组合数的FFT 任意基数的FFT算法 以上讨论的都是以2为基数的FFT算法 即N 2M 这种情况实际上使用得最多 优点 程序简单 效率高 使用方便 实际应用时 有限长序列的长度N很大程度上由人为因素确定 因此多数场合可取N 2M 从而直接使用以2为基数的FFT算法 如N不能人为确定 N的数值也不是以2为基数的整数次方 处理方法有两种 补零 将x n 补零 使N 2M 例如N 30 补上x 30 x 31 0两点 使N 32 25 这样可直接采用以2为基数M 5的FFT程序 有限长度序列补零后并不影响其频谱X ejw 只是频谱的采样点数增加了 上例中由30点增加到32点 所以在许多场合这种处理是可接受的 如要求准确的N点DFT值 可采用任意数为基数的FFT算法 其计算效率低于以2为基数FFT算法 如N为复合数 可分解为两个整数p与q的乘积 像前面以2为基数时一样 FFT的基本思想是将DFT的运算尽量分小 因此 在N pq情况下 也希望将N点的DFT分解为p个q点DFT或q个p点DFT 以减少计算量 步骤 将DFT的时间顺序n和频率顺序k分别表示成二维形式 分别为0 1 Q 1 分别为0 1 P 1 N点DFT可以重新写成为 考虑到 再令 令 n0是参变量 n0为参变量 n1与k0为变量的P点DFT 以P 3 Q 4 N 12为例 1 先将x n 通过x n1Q n0 改写成x n1 n0 因为Q 4 n1 0 1 2 n0 0 1 2 3 故输入是按自然顺序的 即 x 0 0 x 0 x 0 1 x 1 x 0 2 x 2 x 0 3 x 3 x 1 0 x 4 x 1 1 x 5 x 1 2 x 6 x 1 3 x 7 x 2 0 x 8 x 2 1 x 9 x 2 2 x 10 x 2 3 x 11 2 求 个 点的DFT X1 k0 n0 乘以旋转因子得到X1 k0 n0 求P个 点的DFT 参变量是k0 将X2 k0 k1 通过X k0 k1P 恢复为X k N 12为组合数时的FFT 求 个P点DFT需要 P2次复数乘法和 P P 1 次复数加法 乘 个W因子需要 次复数乘法 求P个 点DFT需要PQ2次复数乘法和P 1 次复数加法 总的复数乘法量 QP2 N PQ2 N P Q 1 总的复数加法量 Q P P 1 P Q Q 1 N P Q 2 上述分解原则可推广至任意基数的更加复杂的情况 例如 如果N可分解为m个质数因子p1 p2 pm 即N p1p2p3 pm则第一步 可把N先分解为两个因子N p1q1 其中q1 p2p3 pm 并用上述讨论的方法将DFT分解为p1个q1点DFT 第二步 将q1分解为q1 p2q2 q2 p3p4 pm 然后将每个q1点DFT再分解为p2个q2点DFT 依此类推 通过m次分解 一直分到最少点数的DFT运算 从而获得最高的运算效率 其运算量近似为N p1 p2 pm 次复数乘法和复数加法 但计算效率的提高是要以编程的复杂性为代价的 一般较少应用 p1 p2 pm 2 为基2FFT算法 当组合数N P1P2P3 Pm中所有的Pi均为4时 就是基四FFT算法 以N 43为例 第一级运算的一般形式为 基四FFT只需要级间旋转因子的乘法运算每一级需要N次乘法 提出原因 特点 公式表达快速算法 卷积形式理解 运算量比较 5 Chirp z变换 采用FFT可以算出全部N点DFT值 即z变换X z 在z平面单位圆上的等间隔取样值 但要求N为复合数 问题的提出 不需要计算整个单位圆上z变换的取样 如对于窄带信号 只需要对信号所在的一段频带进行分析 这时 希望频谱的采样集中在这一频带内 以获得较高的分辨率 而频带以外的部分可不考虑 对其它围线上的z变换取样感兴趣 例如语音信号处理中 需要知道z变换的极点所在频率 如极点位置离单位圆较远 则其单位圆上的频谱就很平滑 如果采样不是沿单位圆而是沿一条接近这些极点的弧线进行 则在极点所在频率上的将出现明显的尖峰 由此可较准确地测定极点频率 要求能有效地计算当N是素数时序列的DFT 算法原理 螺旋线采样是一种适合于这种需要的变换 且可以采用FFT来快速计算 这种变换也称作Chirp z变换 已知x n 0 n N 1令zk AW k k 0 M 1 M 采样点数 A W 任意复数其中 A0表示起始取样点的半径长度 通常A0 1 0表示起始取样点z0的相角 0表两相邻点之间的等分角W0螺旋线的伸展率 W01则线内缩 反时针 W0 1则表示半径为A0的一段圆弧 若A0 1 这段圆弧则是单位圆的一部分 图螺旋线采样 计算z变换在采样点zk的值k 0 1 M 1显然 按照以上公式计算出全部M点采样值需要NM次复乘和 N 1 M次复加 当N及M较大时 计算量迅速增加 以上运算可转换为卷积形式 从而可采用FFT进行 这样可大大提高计算速度 利用zk的表示式代入 则得 nk可以用以下表示式来替换 则 定义 则 上式说明 如对信号x n 先进行一次加权处理 加权系数为 然后通过一个单位脉冲响应为h n 的线性系统 最后 对该系统的前M点输出再作一次的加权 就可得到全部M点螺旋线采样值 系统的单位脉冲响应与频率随时间成线性增加的线性调频信号相似 因此称为Chirp z变换 由于输入信号g n 是有限长的 长为N 但序列是无限长的 而计算0 M 1点卷积g k h k 所需要的h n 是取值在n N 1 M 1 那一部分的值 即L点 那么g n 和h n 两个序列线性卷积时需要延拓的周期为2N M 2 但是 由于我们只需要M点X k 对以后的卷积结果是否有混叠失真并不感兴趣 这样循环卷积的点数 即延拓周期可以仍是L Chirp z变换的计算步骤 对h n 以L为周期作周期延拓 再取其主值序列得到对长为N点的x n 求M点z变换 用FFT时两序列长度应扩展为L N M 1 1 求h n 的主值序列 2 用FFT求的傅里叶变换 H r FFT L点 3 对x n 加权并补零 得 4 G r FFT g n L点 5 Y r G r H r L点 6 y k IFFT Y r L点 7 0 k M 1 利用FFT计算Chirp z变换 计算量估算 乘法数估计 1 2 步可以事先计算 不必实时计算 3 7 两步两次加权共计N M次复乘 形成Y k 需L次复乘 一个FFT与一个IFFT需Llog2L次乘 所以总乘法数 L N M Llog2L 直接计算乘法数NM 当N及M较大时 用FFT实现Chirp Z变换 速度上有很大的改进 Chirp z变换的特点 1 输入序列的长度N与输出序列的长度M不需要相等 2 N及M不必是复合数 二者均可为素数 3 相邻采样点zk之间的角间隔 0是任意的 即频率分辨率是任意的 4 围线是任意的 不必是Z平面上的圆 5 起始点z0可任意选定 即可从任意频率上开始对输入数据进行窄带高分辨率分析 6 若A 1 M N 可用Chirp z变换计算DFT 即使N为素数 2 4FFT应用中的几个问题 1 IDFT的运算方法 以上所讨论的FFT算法可用于IDFT运算 简称为IFFT比较IDFT的定义式 IDFT DFT a IDFT与DFT的差别 1 把DFT中的每一个系数改为 2 再乘以常数1 N 则以上所讨论的时间抽取或频率抽取的FFT运算均可直接用于IDFT运算 当然 蝶形中的系数应改为 b 第二种方法 完全不需要改动FFT程序 而是直接利用它作IFFT 考虑到故IFFT计算分三步 将X k 取共轭 虚部乘以 1 对直接作FFT 对FFT的结果取共轭并乘以1 N 得x n 2 实数序列的FFT以上讨论的FFT算法都是复数运算 包括序列x n 也认为是复数 但大多数场合 信号是实数序列 任何实数都可看成虚部为零的复数 例如 求某实信号x n 的复谱 可认为是将实信号加上数值为零的虚部变成复信号 x n j0 再用FFT求其离散傅里叶变换 这种作法很不经济 因为把实序列变成复序列 存储器要增加一倍 且计算机运行时 即使虚部为零 也要进行涉及虚部的运算 浪费了运算量 合理的解决方法是利用复数据FFT对实数据进行有效计算 下面介绍两种方法 a 用一个N点FFT同时计算两个N点实序列的DFT设x n y n 是彼此独立的两个N点实序列 且X k DFT x n Y k DFT y n 则X k Y k 可通过一次FFT运算同时获得 首先将x n y n 分别当作一复序列的实部及虚部 令g n x n jy n 通过FFT运算可获得g n 的DFT值 利用离散傅里叶变换的共轭对称性 通过g n 的FFT运算结果G k 由上式可得到X k 的值 同理 可得 因此 作一次 点复序列的FFT 再通过加 减法运算就可以将X k 与Y k 分离出来 显然 这将使运算效率提高一倍 b 用一个N点的FFT运算获得一个2N点实序列的DFT设x n 是2N点的实序列 现人为地将x n 分为偶数组x1 n 和奇数组x2 n x1 n x 2n x2 n x 2n 1 n 0 1 N 1然后将x1 n 及x2 n 组成一个复序列 y n x1 n jx2 n 通过N点FFT运算可得到 Y k X1 k jX2 k 根据前面的讨论 得到 为求2N点x n 所对应X k 需求出X k 与X1 k X2 k 的关系 1 由x1 n 及x2 n 组成复序列 经FFT运算求得Y k 2 利用共轭对称性求出X1 k X2 k 3 最后利用上式求出X k 达到用一个N点的FFT计算一个2N点的实序列的DFT的目的 3 线性卷积的FFT算法线性卷积是求离散系统响应的主要方法之一 许多重要应用都建立在这一理论基础上 如卷积滤波等 以前曾讨论了用圆周卷积计算线性卷积的方法归纳如下 将长为N2的序列x n 延长到L 补L N2个零 将长为N1的序列h n 延长到L 补L N1个零 如果L N1 N2 1 则圆周卷积与线性卷积相等 此时 可用FFT计算线性卷积 方法如下 a 计算X k FFT x n b 求H k FFT h n c 求Y k H k X k k 0 L 1d 求y n IFFT Y k n 0 L 1可见 只要进行二次FFT 一次IFFT就可完成线性卷积计算 计算表明 L 32时 上述计算线性卷积的方法比直接计算线卷积有明显的优越 因此 也称上述循环卷积方法为快速卷积法 例x1 n 1 2 3 x2 n 5 4 3 2 循环卷积和线性卷积 L 4 3 1时 线性卷积等于循环卷积若L 6 则由于循环卷绕 循环卷积和线性卷积结果存在误差线性卷积重叠相加得到循环卷积结果 5 13 8 14 6 8 x1 1 2 3 x2 5 4 3 2 N 8 xk1 fft x1 N xk2 fft x2 N yk xk1 xk2 y ifft yk N N 4 y 8 8 20 4 N 8 y 5 14 20 4 13 6 0 0 上述结论适用于x n h n 两序列长度比较接近或相等的情况 如果 1 x n h n 长度相差较多 例如 h n 为某滤波器的单位脉冲响应 长度有限 用来处理一个很长的输入信号x n 按上述方法 h n 要补许多零再进行计算 计算量有很大的浪费 2 若x n 为外部采集数据 很长的x n 则意味着占用大量的存储单元和处理时间的延误 对系统的性能均带来不利的影响 为了保持快速卷积法的优越性 可将x n 分为许多段 每段的长度与h n 接近 处理方法有两种 h n x n 假定xi n 表示x n 序列的第i段 则输入序列可表为 于是输出可分解为 其中 a 重叠相加法 1 先对h n 及xi n 补零 具有L点长度 L N1 N2 1 一般选L 2M 2 用基2FFT计算yi n xi n h n 3 重叠部分相加构成最后的输出序列 由于yi n 的长度为N 而xi n 的长度为N2 因此相邻两段yi n 序列必然有L N2 N1 1点发生重叠 N1 h n 长度N2 分段长度L 循环卷积长度 补零相加 计算步骤 a 事先准备好滤波器参数H k DFT h n N点b 用L点FFT计算Xi k DFT xi n c Yi k Xi k H k d 用N点IFFT求yi n IDFT Yi k e 将重叠部分相加函数fftfilt 由分段卷积的各段相加构成总的卷积输出 重叠 卷积结果末端相加 卷积结果 b 重叠保留法 这种方法分段序列不是补零 除了第一段 补N1 1个0 而是保留原来的前一段留下来的输入序列值 此时只是简单地将x n 分成长度为L的段 每一段的前段有N1 1个点与前一段相同 如果L 2M 则可在每段序列末端补零 h n 的线性移位 长度为N1 和循环移位 长度为L 在0 m N1 2处结果不同 h n 的循环移位在主值序列末端出现非零值 而此处分段是有数值存在的 所以在0 n N1 2这一部分的循环卷积结果将混入上述乘积值 从而不同于线性卷积的结果 但是从n N1 1开始 二者结果完全相同 所以需要舍去循环卷积开头的N1 1点 各输出段留下来的序列衔接起来 就构成最后的输出结果 重叠 卷积输入末端保留 输入舍去 输出前段 重叠保留法与重叠相加法的计算量差不多 但省去了重叠相加法最后的相加运算 4 用FFT计算相关函数相关的概念很重要 互相关运算广泛应用于信号分析与统计分析 如通过相关函数峰值的检