丁玉美_数字信号处理_第4章_快速傅里叶变换.ppt

第4章快速傅里叶变换 FFT 4 1引言4 2基2FFT算法4 3进一步减少运算量的措施4 4分裂基FFT算法4 5离散哈特莱变换 DHT 4 1引言 DFT是信号分析与处理中的一种重要变换 因直接计算DFT的计算量与变换区间长度N的平方成正比 当N较大时 计算量太大 所以在快速傅里叶变换 简称FFT 出现以前 直接用DFT算法进行谱分析和信号的实时处理是不切实际的 直到1965年发现了DFT的一种快速算法以后 情况才发生了根本的变化 4 2基2FFT算法 4 2 1直接计算DFT的特点及减少运算量的基本途径长度为N的有限长序列x n 的DFT为考虑x n 为复数序列的一般情况 对某一个k值 直接按 4 2 1 式计算X k 值需要N次复数乘法 N 1 次复数加法 4 2 1 如前所述 N点DFT的复乘次数等于N2 显然 把N点DFT分解为几个较短的DFT 可使乘法次数大大减少 另外 旋转因子WmN具有明显的周期性和对称性 其周期性表现为 4 2 2 其对称性表现为 或者 4 2 2时域抽取法基2FFT基本原理FFT算法基本上分为两大类 时域抽取法FFT DecimationInTimeFFT 简称DIT FFT 和频域抽取法FFT DecimationInFrequencyFFT 简称DIF FFT 下面先介绍DIF FFT算法 设序列x n 的长度为N 且满足 为自然数 按n的奇偶把x n 分解为两个N 2点的子序列 则x n 的DFT为 由于 所以 其中X1 k 和X2 k 分别为x1 r 和x2 r 的N 2点DFT 即 4 2 5 4 2 6 由于X1 k 和X2 k 均以N 2为周期 且 所以X k 又可表示为 4 2 7 4 2 8 图4 2 1蝶形运算符号 图4 2 2N点DFT的一次时域抽取分解图 N 8 与第一次分解相同 将x1 r 按奇偶分解成两个N 4长的子序列x3 l 和x4 l 即 那么 X1 k 又可表示为 4 2 9 式中 同理 由X3 k 和X4 k 的周期性和WmN 2的对称性Wk N 4N 2 WkN 2最后得到 4 2 10 用同样的方法可计算出 4 2 11 其中 图4 2 3N点DFT的第二次时域抽取分解图 N 8 图4 2 4N点DIT FFT运算流图 N 8 4 2 3DIT FFT算法与直接计算DFT运算量的比较每一级运算都需要N 2次复数乘和N次复数加 每个蝶形需要两次复数加法 所以 M级运算总共需要的复数乘次数为 复数加次数为 例如 N 210 1024时 图4 2 5FFT算法与直接计算DFT所需乘法次数的比较曲线 4 2 4DIT FFT的运算规律及编程思想1 原位计算由图4 2 4可以看出 DIT FFT的运算过程很有规律 N 2M点的FFT共进行M级运算 每级由N 2个蝶形运算组成 2 旋转因子的变化规律如上所述 N点DIT FFT运算流图中 每级都有N 2个蝶形 每个蝶形都要乘以因子WpN 称其为旋转因子 p称为旋转因子的指数 观察图4 2 4不难发现 第L级共有2L 1个不同的旋转因子 N 23 8时的各级旋转因子表示如下 L 1时 WpN WJN 4 WJ2L J 0L 2时 WpN WJN 2 WJ2L J 0 1L 3时 WpN WJN WJ2L J 0 1 2 3对N 2M的一般情况 第L级的旋转因子为 4 2 12 4 2 13 3 蝶形运算规律设序列x n 经时域抽选 倒序 后 存入数组X中 如果蝶形运算的两个输入数据相距B个点 应用原位计算 则蝶形运算可表示成如下形式 X J XL 1 J XL 1 J B WpNXL J B XL 1 J XL 1 J B WpN式中p J 2M L J 0 1 2L 1 1 L 1 2 M 下标L表示第L级运算 XL J 则表示第L级运算后数组元素X J 的值 如果要用实数运算完成上述蝶形运算 可按下面的算法进行 设T XL 1 J B WpN TR jTIXL 1 J X R J jX I J 式中下标R表示取实部 I表示取虚部 则 4 编程思想及程序框图 图4 2 6DIT FFT运算和程序框图 5 序列的倒序DIT FFT算法的输入序列的排序看起来似乎很乱 但仔细分析就会发现这种倒序是很有规律的 由于N 2M 所以顺序数可用M位二进制数 nM 1nM 2 n1n0 表示 图4 2 7形成倒序的树状图 N 23 表4 2 1顺序和倒序二进制数对照表 图4 2 8倒序规律 图4 2 9倒序程序框图 4 2 5频域抽取法FFT DIF FFT 在基2快速算法中 频域抽取法FFT也是一种常用的快速算法 简称DIF FFT 设序列x n 长度为N 2M 首先将x n 前后对半分开 得到两个子序列 其DFT可表示为如下形式 偶数 奇数 将X k 分解成偶数组与奇数组 当k取偶数 k 2r r 0 1 N 2 1 时 4 2 14 当k取奇数 k 2r 1 r 0 1 N 2 1 时 4 2 15 将x1 n 和x2 n 分别代入 4 2 14 和 4 2 15 式 可得 4 2 16 图4 2 10DIF FFT蝶形运算流图符号 图4 2 11DIF FFT一次分解运算流图 N 8 图4 2 12DIF FFT二次分解运算流图 N 8 图4 2 13DIF FFT运算流图 N 8 图4 2 14DIT FFT的一种变形运算流图 图4 2 15DIT FFT的一种变形运算流图 4 2 6IDFT的高效算法上述FFT算法流图也可以用于离散傅里叶逆变换 InverseDiscreteFourierTransform 简称IDFT 比较DFT和IDFT的运算公式 图4 2 16DIT IFFT运算流图 图4 2 17DIT IFFT运算流图 防止溢出 如果希望直接调用FFT子程序计算IFFT 则可用下面的方法 由于 对上式两边同时取共轭 得 4 3进一步减少运算量的措施 4 3 1多类蝶形单元运算由DIT FFT运算流图已得出结论 N 2M点FFT共需要MN 2次复数乘法 由 4 2 12 式 当L 1时 只有一种旋转因子W0N 1 所以 第一级不需要乘法运算 综上所述 先除去第一 二两级后 所需复数乘法次数应是从L 3至L M共减少复数乘法次数为 4 3 1 4 3 2 因此 DIT FFT的复乘次数降至 4 3 3 从实数运算考虑 计算N 2M点DIT FFT所需实数乘法次数为 4 3 4 4 3 2旋转因子的生成在FFT运算中 旋转因子WmN cos 2 m N jsin 2 m N 求正弦和余弦函数值的计算量是很大的 4 3 3实序列的FFT算法设x n 为N点实序列 取x n 的偶数点和奇数点分别作为新构造序列y n 的实部和虚部 即 对y n 进行N 2点FFT 输出Y k 则 根据DIT FFT的思想及式 4 2 7 和 4 2 8 可得到 由于x n 为实序列 所以X k 具有共轭对称性 X k 的另外N 2点的值为 4 4分裂基FFT算法 4 4 1分裂基FFT算法原理当n pq 且p N 4 q 4时 n可表示为 并有 再将上式中的k表示为 可得 对k0 0 1 2 3 并用k表示k1 用n表示n0 可以写出 4 4 1 4 4 2 4 4 3 令 则 4 4 2 式可写成如下更简明的形式 4 4 4 图4 4 1分裂基第一次分解L形流图 例如 N 16 第一次抽选分解时 由式 4 4 3 得x2 n x n x n 8 0 n 7x14 n x n x n 8 j x n 4 x n 12 Wn16 0 n 3x24 n x n x n 8 j x n 4 x n 12 W3n16 0 n 3把上式代入式 4 4 4 可得X 2k DFT x2 n 0 k 7X 4k 1 DFT x14 n 0 k 3X 4k 3 DFT x24 n 0 k 3 图4 4 2分裂基FFT算法L形排列示意图与结构示意图 a 分裂基FFT算法L形排列示意图 b 分裂基FFT算法运算流图结构示意图 图4 4 316点分裂基第一次分解L形流图 图中省去箭头 第二次分解 先对图4 4 3中N 2点DFT进行分解 令X1 l X 2l 则有X1 2l DFT y2 n 0 l 3X1 4l 1 DFT y14 n 0 l 1X1 4l 3 DFT y24 n 0 l 1 其中y2 n x2 n x2 n 4 0 n 3y14 n x2 n x2 n 4 x2 n 2 x n 6 Wn8 n 0 1y24 n x2 n x2 n 4 j x2 n 2 x2 n 6 W3n8 n 0 1 图4 4 4图4 4 4中N 2点DFT的分解L形流图 图4 4 54点分裂基L形运算流图 图4 4 616点分裂基FFT运算流图 4 4 2分裂基FFT算法的运算量设第j级有lj个L形 j 1 2 M 1 M log2N 则有l1 N 4 由图4 4 2 b 可见 第j 1列中的L形包含了第j列中的一部分结点的计算 即空白部分 所占结点数刚好等于第j 1列中所有L形对应结点的一半 所以第j列L形个数就减少lj 1 2个 即 由于每个L形有两次复 数 乘运算 所以全部复乘次数为 4 4 5 4 5离散哈特莱变换 DHT 4 5 1离散哈特莱变换定义设x n n 0 1 N 1 为一实序列 其DHT定义为 式中 cas cos sin 其逆变换 IDHT 为 4 5 3 逆变换证明如下 4 5 4 将式 4 5 2 代入式 4 5 3 得 4 5 2DHT与DFT之间的关系为了便于比较 重写DFT如下 4 5 5 4 5 6 容易看出 DHT的核函数DFT的核函数的实部与虚部之和 将XH k 分解为奇对称分量XHo k 与偶对称分量XHe k 之和 4 5 7 4 5 8 4 5 9 由DHT定义有 4 5 10a 4 5 10b 所以 x n 的DFT可表示为同理 x n 的DHT可表示为因此 已知x n 的DHT 则DFT可用下式求得 4 5 11 4 5 12 4 5 13 4 5 3DHT的主要优点 1 DHT是实值变换 在对实信号或实数据进行谱分析时避免了复数运算 从而提高了运算效率 相应的硬件也更简单 更经济 2 DHT的正 反变换 除因子1 N外 具有相同的形式 因而 实现DHT的硬件或软件既能进行DHT 也能进行IDHT 3 DHT与DFT间的关系简单 容易实现两种谱之间的相互转换 4 5 4DHT的性质1 线性性质 4 5 14 2 x N n 的DHT 4 5 15 其中 当k 0时 XH N k XH N XH 0 证明由DHT定义 而 3 循环移位性质 4 5 16 4 5 17 证明由DHT定义有 4 奇偶性奇对称序列和偶对称序列的DHT仍然是奇对称序列或偶对称序列 即DHT不改变序列的奇偶性 5 循环卷积定理 4 5 18 4 5 19 证明下面利用DFT的循环卷积定理和DHT与DFT之间的关系来证明其中 X1 k DFT x1 n X2 k DFT x2 n 根据DHT与DFT之间的关系 则有 将上面两式代入式 4 5 20 并整理后 得 所以式 4 5 18 成立 同理可证明式 4 5 19 亦成立 当x1 n 或x2 n 是偶对称序列时 则由DHT的奇偶性有 4 5 21 4 5 5DHT的快速算法 FHT 1 基2DIT FHT算法及运算流图仿照快速DFT的分解方法 可通过时域抽取或频域抽取的方式实现快速DHT x n 的N 2M点DHT如下式 对x n 进行奇偶抽取 4 5 22 4 5 23 与DFT的时域抽取分解比较 不是一个指数函数 所以处理要比W 2r 1 kN麻烦一些 根据三角公式有 4 5 24 4 5 25 令X0H k DHT x0 n X1H k DHT x1 n 并考虑DHT的周期性 4 5 25 式可写成 为了使以下推导中公式简明 记C k cos 2 N k S k sin 2 N k 将式 4 5 26 中的k分别取k N 2 k N 2 k和N k四个值 并考虑X0H k 和X1H k 以N 2为周期 得到 4 5 27 当k 0时 式 4 5 27 中有重复 可单独写成 4 5 28 同理 在k N 4时有 4 5 29 图4 5 1基2DIT FHT原理和哈特莱碟形 图4 5 24点DFHT蝶形图 图4 5 316点基2DIT FHT流图 2 基2D