离散傅里叶变换及其快速算法.ppt

第二章离散傅里叶变换及其快速算法 1753年 Bernoulli就推断一振动的弦可以表示成正弦加权和的形式 但是他未能给出所需的加权系数 Jean Baptiste JosephFourier于1768年3月出生在法国的Auxerre 当他8岁时不幸成了一名孤儿 被收养在一个宗教界主办的军事学校中 在此期间 Fourier对数学产生了浓厚的兴趣 21岁那年 Fourier在巴黎学术界论述了有关数值方程解的著名论作 这一工作使他在巴黎的数学界出名 Fourier不仅是公认的大数学家 而且他还是一位杰出的教师 他灵活运用历史典故使得他的讲座非常生动 实际上 Fourier所研究的主要领域是数学史 Fourier是最早以应用的眼光来解释抽象数学概念的研究者之一 1798年 拿破仑侵略埃及 在侵略队伍中一些有名的数学家和科学家 Fourier就是其中的一位 他负责组织修建第一条从格勒诺布尔到都灵的道路 Fourier也是一个拥有独特想法的一个怪才 例如 他认为酷热是理想的环境 因此 他喜欢居住在严热的小屋里 并穿上厚厚的衣服 1801 法国决心召回自己的军队 于是Fourier才得以重返家园 回国后 Fourier被任命为格勒诺布尔伊泽尔省的长官 就是在此期间 Fourier完成了其经典之作Theorieanalytiquedelachaleur 热能数学原理 在该著作中 他证明了任一周期函数都可以表示成正弦函数和的形式 其中正弦函数的频率为频率的整数倍 Fourier 离散傅里叶变换不仅具有明确的物理意义 相对于DTFT他更便于用计算机处理 但是 直至上个世纪六十年代 由于数字计算机的处理速度较低以及离散傅里叶变换的计算量较大 离散傅里叶变换长期得不到真正的应用 快速离散傅里叶变换算法的提出 才得以显现出离散傅里叶变换的强大功能 并被广泛地应用于各种数字信号处理系统中 近年来 计算机的处理速率有了惊人的发展 同时在数字信号处理领域出现了许多新的方法 但在许多应用中始终无法替代离散傅里叶变换及其快速算法 二 DFS离散付里级数的推导意义 用数字计算机对信号进行频谱分析时 要求信号必须以离散值作为输入 而且上面讨论可知 只有第四种形式 DFS 对数字信号处理有实用价值 但如果将前三种形式要么在时域上采样 要么在频域上采样 变成离散函数 就可以在计算机上应用 所以我们要先了解如何从以上三种形式推出DFS 1 由非周期连续时间信号推出DFS X t 经过抽样为x nT 对离散的时间信号进行DTFT得到周期连续频谱密度函数 再经过抽样 得到周期性离散频谱密度函数即为DFS x t t 取样 x t t DTFT X ej T 采样 X ejw w 2 周期性连续时间信号函数 周期性连续时间信号函数经采样后 得到周期性的离散时间函数 DFS x t X ejw t w 采样 3 非周期离散时间信号 非周期离散时间信号经过序列付里时变换 即单位园上的Z变换 DTFT 得到周期连续谱密度函数 再经采样为周期离散频谱密度函数 DFS x t t X ej T w X ejw DTFT 采样 三 推导DFS正变换 以下由第三种付里叶级数形式为例推导出离散付里叶级数变换 非周期信号x n 其DTFT 单位园上Z变换 为其为周期连续频谱密度函数 对其进行采样 使其成为周期性离散频谱函数 设在一周期内采样N个点 则两采样点间距为 得到频间距为 代入DTFT式子中得 四 DFS的反变换 即证 证明 已知两边同乘以 并对一个周期求和用n置换r得 根据正交定理 回顾DFS 设x n 为周期为N的周期序列 则其离散傅里叶级数 DFS 变换对为 正变换反变换其中 五 离散付里叶级数性质 可以由抽样Z变换来解析DFS 它的许多性质与Z变换性质类似 它们与Z变换主要区别为 1 与两者具有周期性 与Z变换不同 2 DFS在时域和频域之间具有严格的对偶关系 它们主要性质分为 线性 序列移位 循环移位 调制性 周期卷积和 1 线性 2 序列移位 循环 移位 时域频域 3 调制性 4 时域卷积 周期卷积和与以前卷积不同 它的卷积过程限在一个周期内称为周期卷积 时域卷积等于频域相乘 频域 5 频域卷积 时域 证 同理 对于复序列其共轭序列满足 3 共轭对称性 已知序列x1 n R4 n x2 n n 1 R5 n 分别将序列以周期为6周期延拓成周期序列 本节涉及内容周期序列的离散傅立叶级数 DFS 正变换 反变换 2 1 2离散傅里叶变换 DFT 我们知道周期序列实际上只有有限个序列值有意义 因此它的许多特性可推广到有限长序列上 一个有限长序列x n 长为N 为了引用周期序列的概念 假定一个周期序列 它由长度为N的有限长序列x n 延拓而成 它们的关系 周期序列的主值区间与主值序列 对于周期序列 定义其第一个周期n 0 N 1 为的 主值区间 主值区间上的序列为主值序列x n x n 与的关系可描述为 数学表示 RN n 为矩形序列 符号 n N是余数运算表达式 表示n对N求余数 例 是周期为N 8的序列 求n 11和n 2对N的余数 因此 频域上的主值区间与主值序列 周期序列的离散傅氏级数也是一个周期序列 也可给它定义一个主值区间 以及主值序列X k 数学表示 长度为N的有限长序列x n 其离散傅里叶变换X k 仍是一个长度为N的有限长序列 它们的关系为 x n 与X k 是一个有限长序列离散傅里叶变换对 已知x n 就能唯一地确定X k 同样已知X k 也就唯一地确定x n 实际上x n 与X k 都是长度为N的序列 复序列 都有N个独立值 因而具有等量的信息 有限长序列隐含着周期性 DFT的矩阵方程表示 DFT特性 以下讨论DFT的一些主要特性 这些特性都与周期序列的DFS有关 假定x n 与y n 是长度为N的有限长序列 其各自的离散傅里叶变换分别为 X k DFT x n Y k DFT y n 1 线性DFT ax n by n aX k bY k a b为任意常数 2 循环移位有限长序列x n 的循环移位定义为 f n x n m NRN n 含义 1 x n m N表示x n 的周期延拓序列的移位 2 x n m NRN n 表示对移位的周期序列x n m N取主值序列 所以f n 仍然是一个长度为N的有限长序列 f n 实际上可看作序列x n 排列在一个N等分圆周上 并顺时针旋转m位 循环移位 循环移位 移位前 左移两位后 证 利用周期序列的移位特性 实际上 利用WN mk的周期性 将f n x n m NRN n 代入DFT定义式 同样很容易证明 序列循环移位后的DFT为F k DFT f n x k 同样 对于频域有限长序列X k 的循环移位 有如下反变换特性 IDFT X k l NRN k x n 3 循环卷积若F k X k Y k 则或 证 这个卷积可看作是周期序列卷积后再取其主值序列 将F k 周期延拓 得 则根据DFS的周期卷积公式 因0 m N 1时 x m N x m 因此经过简单的换元可证明 这一卷积过程与周期卷积比较 过程是一样的 只是这里只取结果的主值序列 由于卷积过程只在主值区间0 m N 1内进行 所以实际上就是y m 的圆周移位 称为 循环卷积 习惯上常用符号 表示循环卷积 以区别于线性卷积 同样 若f n x n y n 则 4 有限长序列的线性卷积与循环卷积 循环卷积的应用 实际问题的大多数是求解线性卷积 如信号x n 通过系统h n 其输出就是线性卷积y n x n h n 而循环卷积比起线性卷积 在运算速度上有很大的优越性 它可以采用快速傅里叶变换 FFT 技术 若能利用循环卷积求线性卷积 会带来很大的方便 现在我们来讨论上述x n 与h n 的线性卷积 如果x n h n 为有限长序列 则在什么条件下能用循环卷积代替而不产生失真 有限长序列的线性卷积 假定x n 为有限长序列 长度为N y n 为有限长序列 长度为M 它们的线性卷积f n x n y n 也应是有限长序列 因x m 的非零区间 0 m N 1 y n m 的非零区间 0 n m M 1 这两个不等式相加 得 0 n N M 2 在这区间以外不是x m 0 就是y n m 0 因而f n 0 因此 f n 是一个长度为N M 1的有限长序列 循环卷积 重新构造两个有限长序列x n y n 长度均为L max N M 序列x n 只有前N个是非零值 后L N个为补充的零值 序列y n 只有前M个是非零值 后L M个为补充的零值 为了分析x n 与y n 的循环卷积 先看x n y n 的周期延拓 根据前面的分析 f n 具有N M 1个非零序列值 因此 如果周期卷积的周期L N M 1 那么f n 周期延拓后 必然有一部分非零序列值要重叠 出现混淆现象 只有L N M 1时 才不会产生交叠 这时f n 的周期延拓中每一个周期L内 前N M 1个序列值是f n 的全部非零序列值 而剩下的L N M 1 点的序列则是补充的零值 循环卷积正是周期卷积取主值序列 所以使圆周卷积等于线性卷积而不产生混淆的必要条件是 L N M 1 三 DFT涉及的基本概念 1 主值 主值区间 主值序列 2 移位 线性移位 圆周移位 3 卷积 线性卷积 圆周卷积 4 对称 序列的对称性 序列的对称分量 5 相关 线性相关 圆周相关 1 主值 主值区间 主值序列 主值区间 设有限长序列x n 0 n N 1 将其延拓为周期序列 周期序列长度为N 则的第一个周期n 0到n N 1的区间称为主值区间 主值序列 设有限长序列x n 0 n N 1 将其延拓为周期序列 周期为N 则主值区间内的序列x n 0 n N 1 即为主值序列 2 移位 线性移位 序列沿坐标轴的平移 圆周移位 将有限长序列x n 以长度N为周期 延拓为周期序列 并加以线性移位后 再取它的主值区间上的序列值 m点圆周移位记作 其中 N表示N点周期延拓 1 有限长序列圆周移位的实现步骤 2 例子1 2 1 3 1 0 5 1 周期延拓 N 5时 n x n 2 1 3 1 x n 0 5 2 1 3 1 0 5 1 1 2 0 5 n 2 周期延拓 N 6时 补零加长 2 1 3 1 x n 0 5 2 1 3 1 0 5 1 1 2 3 n 3 2 例子 2 1 3 1 0 5 n x n 3 M 1时 左移 取主值 1 3 1 x n 0 5 2 4 M 2时 右移 取主值 2 1 3 1 n x n 0 5 n 3 卷积 卷积在此我们主要介绍 1 线性卷积 2 圆周卷积 3 圆周卷积与线性卷积的性质对比 1 线性卷积 线性卷积定义 有限长序列x1 n 0 n N1 1 x2 n 0 n N1 1则线性卷积为注意 线性卷积结果长度变为N1 N2 1 2 圆周卷积 令则圆周卷积结果长度不变 为N 圆周卷积的实现步骤 例子线性卷积与圆周卷积步骤比较1 2 3 1 x n 5 4 n 0 N1 5 2 1 3 h n n 0 N2 3 线性卷积 圆周卷积 N 7 补零加长 2 3 1 x k 5 4 k 0 N1 5 2 3 1 x k 5 4 0 N 7 k 例子线性卷积与圆周卷积步骤比较2 2 3 1 h k 0 k 2 线卷积无需周期延拓 而圆周卷积需进行周期延拓 线卷积的反折 圆卷积的反折 并取主值区间 2 3 1 2 3 1 2 3 1 h k k 0 2 3 1 h k k 0 例子线性卷积与圆周卷积步骤比较3 3 平移 2 3 1 h 1 k k 0 2 3 1 h 1 k k 0 4 相乘x k h k 5 1 5x k h 1 k 5 2 4 1 14x k h 2 k 5 3 4 2 3 1 26 2 3 1 x k 5 4 k 0 2 3 1 x k 5 4 0 N 7 k x k h 3 k 4 3 3 2 2 1 20 x k h 4 k 3 3 2 2 1 1 14x k h 5 k 2 3 1 2 8x k h 6 k 1 3 3 例子线性卷积与圆周卷积步骤比较4 5 相加得到线性卷积的示意图 相加得到圆周卷积的示意图 14 26 5 n y n 20 14 8 3 0 14 26 5 n y n 20 14 8 3 0 可见 线性卷积与圆周卷积相同 当N N1 5 N2 3 1 7时 例子线性卷积与圆周卷积步骤比较5 若圆周卷积取长度为N 5 则求圆周卷积 2 3 1 x k 5 4 0 N 5 k 2 3 1 h k k 0 求得圆周卷积x k h k 5 1 2 3 1 2 13x k h 1 k 5 2 4 1 1 3 17x k h 2 k 5 3 4 2 3 1 26 x k h 3 k 4 3 3 2 2 1 20 x k h 4 k 3 3 2 2 1 1 14看出圆卷积与线卷积不同 17 13 26 y n n 0 20 14 作业 求 1 x n x n 的线卷积 2 N 4 不加长 3 N 6 补零加长 4 N 7 补零加长 5 N 8 补零加长 1 2 x n 1 2 0 n 3 圆周卷积与线性卷积的性质对比 4 对称 分为 1 序列的对称性 2 序列的对称分量 1 序列的对称性 a 奇对称 序列 和偶对称 序列 b 圆周奇对称 序列 和圆周偶对称 序列 c 共轭对称 序列 和共轭反对称 序列 d 圆周共轭对称 序列 和圆周共轭反对称 序列 a 奇对称 序列 和偶对称 序列 称x n 与 x n 互为奇对称 满足x0 n x0 n 的序列x0 n 称为奇对称序列 称x n 与x n 互为偶对称 满足xe n xe n 的序列xe n 称为偶对称序列 例子 0 xe n n 0 x n n 0 x n n 互为偶对称 为偶对称序列 0 x n n 0 x n n 互为奇对称 0 xo n n 为奇对称序列 b 圆周奇对称 序列 和圆周偶对称 序列 长度为N的有限长序列x n 与 x n NRN n 互为圆周奇对称 长度为N的有限长序列x0 n 若满足x0 n x0 n NRN n 则x0 n 是圆周奇对称序列 长度为N的有限长序列x n 与x n NRN n 互为圆周偶对称 长度为N的有限长序列xe n 若满足xe n xe n NRN n 则是圆周偶对称序列 c 共轭对称 序列 和共轭反对称 序列 序列x n 与x n 互为共轭对称 共轭对称序列是满足xe n x e n 的序列xe n 对于实序列来说 这一条件变成xe n xe n 即为偶对称序列 序列x n 与 x n 互为共轭反对称 共轭反对称序列是满足xe n x o n 的序列 对于实序列来说 即为xo n xo n 奇对称序列 d 圆周共轭对称 序列 和圆周共轭反对称 序列 N点有限长序列x n 与x n NRN n 互为圆周共轭对称 圆周共轭对称序列是满足xep n xep n NRN n 的序列即xep n 的模是圆周偶对称 辐角是圆周奇对称 或说实部圆周偶对称 虚部圆周奇对称 即把xep n 看成分布在N等分的圆上 在n 0的左半圆与右半圆上 序列是共轭对称的 圆周共轭对称 序列 的例子 虚部 实部 实部圆周偶对称 虚部圆周奇对称 圆周共轭反对称 序列 圆周共轭反对称序列是满足xop n xop n NRN n 的序列即xop n 的模是圆周奇对称 辐角是圆周偶对称 或说实部圆周奇对称 虚部圆周偶对称 即把xop n 看成分布在N等分的圆上 在n 0的左半圆与右半圆上 序列是共轭反对称的 圆周共轭反对称 序列 例子 实部 虚部 实部圆周偶对称 虚部圆周奇对称 2 序列的对称分量 a 奇对称分量和偶对称分量 b 圆周奇对称分量和圆周偶对称分量 c 共轭对称分量和共轭反对称分量 d 圆周共轭对称分量和圆周共轭反对称分量 a 奇对称分量和偶对称分量 x n 为任一序列 实或纯虚序列 x n 总能表示成一个奇对称序列xo n 和一个偶对称序列xe n 之和 即x n xo n xe n 其中 xo n 奇对称序列称为x n 的奇对称分量 xe n 偶对称序列称为x n 的偶对称分量 看出这样得到的xo n 和xe n 分别满足奇对称和偶对称的条件 且二者之和为x n 说明 若x n 为有限长序列且0 n N 1 则与点数均为 2N 1 区别于奇对称 序列 和偶对称 序列 b 圆周奇对称分量和圆周偶对称分量 设x n 是一长度为N的有限长序列 总能表示成一个圆周奇对称序列xop n 和一个圆周偶对称序列xep n 之和 即x n xep n xop n 其中xop n 称为x n 的圆周奇对称分量 xep n 称为x n 的圆周偶对称分量 看出满足圆周奇对称和圆周偶对称的条件 且二者之和为x n c 共轭对称分量和共轭反对称分量 任一序列x n 总能表示成一个共轭对称序列xo n 和一个共轭反对称序列xe n 之和 即x n xo n xe n 其中 xo n 共轭反对称序列称为x n 的共轭反对称分量 xe n 共轭对称序列称为x n 的共轭对称分量 看出xo n 和xe n 分别满足奇对称和偶对称的条件 且二者之和为x n d 圆周共轭对称分量和圆周共轭反对称分量 设x n 是一长度为N的有限长序列 总能表示成一个圆周共轭反对称序列xop n 和一个圆周共轭对称序列xep n 之和 即x n xep n xop n 其中xop n 称为x n 的圆周共轭反对称分量 xep n 称为x n 的圆周共轭对称分量 看出满足圆周奇对称和圆周偶对称的条件 且二者之和为x n f 选频性 对 0有限制 对复指数函数进行采样得复序列x n 0 n N 1其中q为整数 当 0 2 N时 x n ej2 nq N 其离散傅里叶变换为写成闭解形式可见 当输入频率为q 0时 变换X K 的N个值中只有X q N 其余皆为零 如果输入信号为若干个不同频率的信号的组合 经离散傅里叶变换后 不同的k上 X k 将有一一对应的输出 因此 离散傅里叶变换算法实质上对频率具有选择性 例 求余弦序列 的DFT g DFT与Z变换有限长序列可以进行z变换比较z变换与DFT变换 可见 当z w kN时 即 图DFT与z变换 o o o o o o o o o o o X ej X k o Re z jIm z o 变量 周期 分辨率 是z平面单位圆上幅角为的点 即将z平面上的单位圆N等分后的第k点 1 X k 也就是z变换在单位圆上等间隔的采样值 2 X k 也可看作是对序列傅氏变换X ej 的采样 采样间隔为 N 2 N 即 结论 采样定律告诉我们 一个频带有限的信号 可以对它进行时域采样而不丢失任何信息 DFT变换进一步告诉我们 对于时间有限的信号 有限长序列 也可以对其进行频域采样 而不丢失任何信息 这正反应了傅立叶变换中时域 频域的对称关系 它有十分重要的意义 由于时域上的采样 使我们能够采用数字技术来处理这些时域上的信号 序列 而DFT的理论不仅在时域 而且在频域也离散化 因此使得在频域采用数字技术处理成为可能 FFT就是频域数字处理中最有成效的一例 8 DFT形式下的Parseval定理 令k 0 得到 2 2利用DFT做连续信号的频谱分析 利用DFT计算连续信号的频谱 1 混迭对连续信号x t 进行数字处理前 要进行采样采样序列的频谱是连续信号频谱的周期延拓 周期为fs 如采样率过低 不满足采样定理 fs 2fh 则导致频谱混迭 使一个周期内的谱对原信号谱产生失真 无法恢复原信号 进一步的数字处理失去依据 2 泄漏处理实际信号序列x n 时 一般总要将它截断为一有限长序列 长为N点 相当于乘以一个矩形窗w n RN n 矩形窗函数 其频谱有主瓣 也有许多副瓣 窗口越大 主瓣越窄 当窗口趋于无穷大时 就是一个冲击函数 我们知道 时域的乘积对应频域的卷积 所以 加窗后的频谱实际是原信号频谱与矩形窗函数频谱的卷积 卷积的结果使频谱延伸到了主瓣以外 且一直延伸到无穷 当窗口无穷大时 与冲激函数的卷积才是其本身 这时无畸变 否则就有畸变 例如 信号为 是一单线谱 但当加窗后 线谱与抽样函数进行卷积 原来在 0处的一根谱线变成了以 0为中心的 形状为抽样函数的谱线序列 原来在一个周期 s 内只有一个频率上有非零值 而现在一个周期内几乎所有频率上都有非零值 即的频率成份从 0处 泄漏 到其它频率处去了 考虑各采样频率周期间频谱 泄漏 后的互相串漏 卷积后还有频谱混迭现象产生 3 栅栏效应N点DFT是在频率区间 0 2 上对信号频谱进行N点等间隔采样 得到的是若干个离散的频谱点X k 且它们限制在基频的整数倍上 这就好像在栅栏的一边通过缝隙看另一边的景象一样 只能在离散点处看到真实的景象 其余部分频谱成分被遮挡 所以称之为栅栏效应 减小栅栏效应方法 尾部补零 使谱线变密 增加频域采样点数 原来漏掉的某些频谱分量就可能被检测出来 4 DFT的分辨率 填补零值可以改变对DTFT的采样密度 人们常常有一种误解 认为补零可以提高DFT的频率分辨率 事实上我们通常规定DFT的频率分辨率为 这里的N是指信号x n 的有效长度 而不是补零的长度 不同长度的x n 其DTFT的结果是不同的 而相同长度的x n 尽管补零的长度不同其DTFT的结果应是相同的 他们的DFT只是反映了对相同的DTFT采用了不同的采样密度 参数选择的一般原则 若已知信号的最高频率 为防止混叠 选定采样频率 根据频率分辩率 确定所需DFT的长度 3 和N确定以后 即可确定相应模拟信号的时间长度这里T是采样周期 5 周期信号的谱分析 对于连续的单一频率周期信号 为信号的频率 可以得到单一谱线的DFT结果 但这是和作DFT时数据的截取长度选得是否恰当有关 截取长度N选得合理 可完全等于的采样 6 8 10 k X k a b c d 不同截取长度的正弦信号及其DFT结果 第二章离散傅里叶变换及其快速算法 2 3快速付里叶变换 FFT FastFourietTransformer 一 快速付里叶变换FFT 有限长序列通过离散傅里叶变换 DFT 将其频域离散化成有限长序列 但其计算量太大 很难实时地处理问题 因此引出了快速傅里叶变换 FFT FFT并不是一种新的变换形式 它只是DFT的一种快速算法 并且根据对序列分解与选取方法的不同而产生了FFT的多种算法 FFT在离散傅里叶反变换 线性卷积和线性相关等方面也有重要应用 二 FFT产生故事 当时加文 Garwin 在自已的研究中极需要一个计算付里叶变换的快速方法 他注意到图基 J W Turkey 正在写有关付里叶变换的文章 因此详细询问了图基关于计算付里叶变换的技术知识 图基概括地对加文介绍了一种方法 它实质上就是后来的著名的库利 CooleyJ W 图基算法 在加文的迫切要求下 库利很快设计出一个计算机程序 1965年库利 图基在 MathematicofComputation杂志上发表了著名的 机器计算付里级数的一种算法 文章 提出一种快速计算DFT的方法和计算机程序 揭开了FFT发展史上的第一页 促使FFT算法产生原因还有1967年至1968年间FFT的数字硬件制成 电子数字计算机的条件 使DFT的运算大简化了 三 本节主要内容 1 直接计算DFT算法存在的问题及改进途径 2 多种DFT算法 时间抽取算法DIT算法 频率抽取算法DIF算法 线性调频Z变换即CZT法 3 FFT的应用 直接计算DFT算法存在的问题及改进途径 一 直接计算DFT计算量 问题提出 设有限长序列x n 非零值长度为N 计算对x n 进行一次DFT运算 共需多大的运算工作量 1 比较DFT与IDFT之间的运算量 其中x n 为复数 也为复数所以DFT与IDFT二者计算量相同 2 以DFT为例 计算DFT复数运算量 计算一个X k 一个频率成分 值 运算量为例k 1则要进行N次复数乘法 N 1 次复数加法所以 要完成整个DFT运算 其计算量为 N N次复数相乘 N N 1 次复数加法 3 一次复数乘法换算成实数运算量 复数运算要比加法运算复杂 需要的运算时间长 一个复数乘法包括4个实数乘法和2个实数相法 a jb c jd ac bd j bc ad 4次复数乘法 2次实数加法 4 计算DFT需要的实数运算量 每运算一个X k 的值 需要进行4N次实数相乘和2N 2 N 1 2 2N 1 次实数相加 整个DFT运算量为 4N2次实数相乘和2N 2N 1 次实数相加 由此看出 直接计算DFT时 乘法次数与加法次数都是和N2成比例的 当N很大时 所需工作量非常可观 例子 例1 当N 1024点时 直接计算DFT需要 N2 220 1048576次 即一百多万次的复乘运算这对实时性很强的信号处理 如雷达信号处理 来讲 它对计算速度有十分苛刻的要求 迫切需要改进DFT的计算方法 以减少总的运算次数 例2 石油勘探 24道记录 每道波形记录长度5秒 若每秒抽样500点 秒 每道总抽样点数 500 5 2500点24道总抽样点数 24 2500 6万点DFT运算时间 N2 60000 2 36 108次 二 改善DFT运算效率的基本途径 利用DFT运算的系数的固有对称性和周期性 改善DFT的运算效率 1 合并法 合并DFT运算中的某些项 3 分解法 将长序列DFT利用对称性和周期性 分解为短序列DFT 利用DFT运算的系数的固有对称性和周期性 改善DFT的运算效率 的对称性 的周期性 因为 由此可得出 例子 例 利用以上特性 得到改进DFT直接算法的方法 1 合并法 步骤1分解成虚实部 合并DFT运算中的有些项对虚实部而言所以带入DFT中 1 合并法 步骤2代入DFT中 展开 1 合并法 步骤3合并有些项 根据 有 1 合并法 步骤4结论 由此找出其它各项的类似归并方法 乘法次数可以减少一半 例 2 将长序更DFT利用对称性和周期性分解为短序列DFT 思路 因为DFT的运算量与N2成正比的如果一个大点数N的DFT能分解为若干小点数DFT的组合 则显然可以达到减少运算工作量的效果 2 将长序更DFT利用对称性和周期性分解为短序列DFT 方法 把N点数据分成二半 其运算量为 再分二半 这样一直分下去 剩下两点的变换 2 将长序更DFT利用对称性和周期性分解为短序列DFT 结论 快速付里时变换 FFT 就是在此特性基础上发展起来的 并产生了多种FFT算法 其基本上可分成两大类 按抽取方法分 时间抽取法 DIT 频率抽取法 DIF 按 基数 分 基 2FFT算法 基 4FFT算法 混合基FFT算法 分裂基FFT算法其它方法 线性调频Z变换 CZT法 2 3 1基 2按时间抽取的FFT算法Decimation in Time DIT Coolkey Tukey 一 算法原理 设输入序列长度为N 2M M为正整数 将该序列按时间顺序的奇偶分解为越来越短的子序列 称为基2按时间抽取的FFT算法 也称为Coolkey Tukey算法 其中基数2 N 2M M为整数 若不满足这个条件 可以人为地加上若干零值 加零补长 使其达到N 2M 例子 设一序列x n 的长度为L 9 应加零补长为N 24 16应补7个零值 0123456789101112131415n x n 二 算法步骤1 分组 变量置换 DFT变换 先将x n 按n的奇偶分为两组 作变量置换 当n 偶数时 令n 2r 当n 奇数时 令n 2r 1 得到 x 2r x1 r x 2r 1 x2 r r 0 N 2 1 则其DFT可化为两部分 前半部分 后半部分 2 代入DFT中 生成两个子序列 x 0 x 2 x 2r 奇数点x 1 x 3 x 2r 1 偶数点 代入DFT变换式 3 求出子序列的DFT 上式得 4 结论1 一个N点的DFT被分解为两个N 2点DFT X1 k X2 k 这两个N 2点的DFT按照 再应用W系数的周期性 求出用X1 k X2 k 表达的后半部的X k N 2 的值 5 求出后半部的表示式 看出 后半部的k值所对应的X1 k X2 k 则完全重复了前半部分的k值所对应的X1 k X2 k 的值 6 结论2 频域中的N个点频率成分为 结论 只要求出 0 N 2 1 区间内的各个整数k值所对应的X1 k X2 k 值 即可以求出 0 N 1 整个区间内全部X k 值 这就是FFT能大量节省计算的关键 7 结论3 由于N 2L 因此N 2仍为偶数 可以依照上面方法进一步把每个N 2点子序列 再按输入n的奇偶分解为两个N 4点的子序列 按这种方法不断划分下去 直到最后剩下的是2点DFT 两点DFT实际上只是加减运算 三 蝶形结 即蝶式计算结构也即为蝶式信号流图上面频域中前 后半部分表示式可以用蝶形信号流图表示 X1 k X2 k 作图要素 1 左边两路为输入 2 右边两路为输出 3 中间以一个小圆表示加 减运算 右上路为相加输出 右下路为相减输出 4 如果在某一支路上信号需要进行相乘运算 则在该支路上标以箭头 将相乘的系数标在箭头旁 5 当支路上没有箭头及系数时 则该支路的传输比为1 例子 求N 23 8点FFT变换 1 先按N 8 N 2 4 做4点的DFT 先将N 8DFT分解成2个4点DFT 可知 时域上 x 0 x 2 x 4 x 6 为偶子序列x 1 x 3 x 5 x 7 为奇子序列频域上 X 0 X 3 由X k 给出X 4 X 7 由X k N 2 给出 a 比较N 8点直接DFT与分解2个4点DFT的FFT运算量 N 8点的直接DFT的计算量为 N2次 64次 复数相乘 N N 1 次 8 8 1 56次 复数相加 共计120次 b 求一个蝶形结需要的运算量 要运算一个蝶形结 需要一次乘法 两次加法 c 分解为两个N 2 4点的DFT的运算量 分解2个N 2点 4点 的DFT 偶数其复数相乘为复数相加为 奇数其复数相乘为复数相加为 d 用2个4点来求N 8点的FFT所需的运算量 再将N 2点 4点 合成N点 8点 DFT时 需要进行N 2个蝶形运算 还需N 2次 4次 乘法及次 8次 加法运算 e 将N 8点分解成2个4点的DFT的信号流图 4点DFT x 0 x 2 x 4 x 6 4点DFT x 1 x 3 x 5 x 7 X 0 X 1 X 2 X 3 X 4 X 5 X 6 X 7 X1 0 X1 1 X1 2 X1 3 X2 0 X2 1 X2 2 X2 3 偶数序列 奇数序列 X 4 X 7 同学们自已写 x1 r x2 r 2 N 2 4点 N 4 2点 FFT a 先将4点分解成2点的DFT 因为4点DFT还是比较麻烦 所以再继续分解 若将N 2 4点 子序列按奇 偶分解成两个N 4点 2点 子序列 即对将x1 r 和x2 r 分解成奇 偶两个N 4点 2点 点的子序列 b 求2点的DFT c 一个2点的DFT蝶形流图 2点DFT 2点DFT x 0 x 4 x 2 x 6 X3 0 X3 1 X4 0 X4 1 X1 0 X1 1 X1 2 X1 3 d 另一个2点的DFT蝶形流图 2点DFT 2点DFT x 1 x 5 x 3 x 7 X5 0 X5 1 X6 0 X6 1 X2 0 X2 1 X2 2 X2 3 同理 3 将N 4 2点 DFT再分解成2个1点的DFT a 求2个一点的DFT 最后剩下两点DFT 它可分解成两个一点DFT 但一点DFT就等于输入信号本身 所以两点DFT可用一个蝶形结表示 取x 0 x 4 为例 b 2个1点的DFT蝶形流图 1点DFT x 0 1点DFT x 4 X3 0 X3 1 进一步简化为蝶形流图 X3 0 X3 1 x 0 x 4 4 一个完整N 8的按时间抽取FFT的运算流图 x 0 x 4 x 2 x 6 x 1 x 5 x 3 x 7 X 0 X 1 X 2 X 3 X 4 X 5 X 6 X 7 m 0 m 1 m 2 四 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级 推进一级 则此系数为后级系数中偶数序号的那一半 4 按时间抽取法 2点DFT 2点DFT 2点DFT 2点DFT 2点DFT 2点DFT 2点DFT 2点DFT 两个2点DFT 两个2点DFT 两个2点DFT 两个2点DFT 两个4点DFT 两个4点DFT 两个N 2点DFT X1 k X2 k X k 由于每一步分解都是基于在每级按输入时间序列的次序是属于偶数还是奇数来分解为两个更短的序列 所以称为 按时间抽取法 x n 五 按时间抽取的FFT算法与直接计算DFT运算量的比较 由前面介绍的按时间抽取的FFT运算流图可见 每级都由N 2个蝶形单元构成 因此每一级运算都需要N 2次复乘和N次复加 每个结加减各一次 这样 N 2M M级运算共需要 复乘次数 复加次数 可以得出如下结论 按时间抽取法所需的复乘数和复加数都是与成正比 而直接计算DFT时所需的复乘数与复加数则都是与N2成正比 复乘数md N2 复加数ad N N 1 N2 例子 看N 8点和N 1024点时直接计算DFT与用基2 按时间抽取法FFT的运算量 看出 当N较大时 按时间抽取法将比直接法快一 二个数量级之多 作业 N 16点的FFT基2 按时间抽取流程图 六 按时间抽取FFT算法的特点 根据DIT基2 FFT算法原理 能得出任何N 2m点的FFT信号流图 并进而得出FFT计算程序流程图 最后总结出按时间抽取法解过程的规律 1 原位运算 in place 2 码位倒读规则 1 原位运算 in place 原位运算的结构 可以节省存储单元 降低设备成本 定义 当数据输入到存储器以后 每一组运算的结果 仍然存放在这同一组存储器中直到最后输出 x 0 x 4 X3 0 X3 1 例子 例 N 8FFT运算 输入 x 0 x 4 x 2 x 6 x 1 x 5 x 3 x 7 A 0 A 1 A 2 A 3 A 4 A 5 A 6 A 7 A 0 A 1 A 2 A 3 A 4 A 5 A 6 A 7 A 0 A 1 A 2 A 3 A 4 A 5 A 6 A 7 A 0 x 0 A 1 x 1 A 2 x 2 A 3 x 3 A 4 x 4 A 5 x 5 A 6 x 6 A 7 x 7 R1 R1 R1 R1 R1 R2 R1 R1 R2 R2 R3 R4 看出 用原位运算结构运算后 A 0 A 7 正好顺序存放X 0 X 7 可以直接顺序输出 2 码位倒读规则 我们从输入序列的序号及整序规律得到码位倒读规则 由N 8蝶形图看出 原位计算时 FFT输出的X k 的次序正好是顺序排列的 即X 0 X 7 但输入x n 都不能按自然顺序存入到存储单元中 而是按x 0 x 4 x 2 x 6 的顺序存入存储单元即为乱序输入 顺序输出 这种顺序看起来相当杂乱 然而它是有规律的 即码位倒读规则 例子 以N 8为例 01234567 000001010011100101110111 自然顺序 二进制码表示 码位倒读 码位倒置顺序 000100010110001101011111 04261537 看出 码位倒读后的顺序刚好是数据送入计算机内的顺序 整序重排子程序 具体执行时 只须将1与4对调送入 3与6对调送入 而0 2 5 7不变 仅需要一个中间存储单元 n01234567 n 04261537 在实际运算时 先按自然顺序将输入序列存入存储单元 再通过变址运算将自然顺序变换成按时间抽取的FFT算法要求的顺序 变址的过程可以用程序安排加以实现 称为 整序 或 重排 采用码位倒读 且注意 1 当n n 时 数据不必调换 2 当n n时 必须将原来存放数据x n 送入暂存器R 再将x n 送入x n R中内容送x n 进行数据对调 3 为了避免再次考虑前面已调换过的数据 保证调换只进行一次 否则又变回原状 n n时 调换 作业 编写N 128点的基2 按时间抽取的FFT算法 要求用C语言编写 并将输入数据放在文件inputdata dat中 输出数据放在outputdata dat文件中 七 按时间抽取的FFT算法的若干变体1 思路 对于任何流图 只要保持各节点所连续的支路及其传输系数不变 则不论节点位置怎么排列所得流图总是等效的 最后所得结果都是x n 的离散付里叶变换的正确结果 只是数据的提取和存放的次序不同而已 2 输入是自然顺序而输出是乱序 将原先中与x 4 水平相邻的所有节点跟x 1 水平相邻的所有节点位置对调 将原先中与x 6 水平相邻的所有节点跟x 3 水平相邻的所有节点位置对调 其余诸节点保持不变 则可得 x 0 x 1 x 2 x 3 x 4 x 5 x 6 x 7 X 0 X 4 X 2 X 6 X 1 X 3 X 5 X 7 它与输入是乱序 输出顺序比较 看出 相同点 运算量一样 不同点 第一是数据存入方式不同 第二是取用系数的顺序不同 2 输入和输出都是自然顺序 对输入自然顺序 输出乱序的第二级 第三级节点作调整 可得输入输出都是顺序 无需数据重排 x 0 x 1 x 2 x 3 x 4 x 5 x 6 x 7 X 0 X 1 X 2 X 3 X 4 X 5 X 6 X 7 1 它失掉了 原位运算 的性质 2 为了变换N点数据 至少需要2N个复数单元 在内存比较紧张时 而对输入数据整序并不困难时一般不用 为了争取速度 才采取这种变体 有一个FIR滤波器处理机 用FFT算法分段过滤信号 每段运算N 1024点 运算一遍需要0 2s 处理机具有两组1024个单元的复数存储器可供交替使用 一组供运算时 另一组可用来存储实时输入的信号序列 用该处理机并配以A D变换器作连续信号的实时过滤 问 1 抽样频率最高可达多少 2 若作两路信号同时过滤时 抽样频率最高是多少 3 在这两种情况下最高可以处理多高频率的信号 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螺旋线的伸展率 W0 1则线外伸 W0 1则线内缩 反时针 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变换 x x x n g n X zk 算法实现 由于输入信号g n 是有限长的 长为N 但序列是无限长的 而计算0 M 1点卷积g k h k 所需要的h n 是取值在n N 1 M 1那一部分的值 因此 可认为h n 是一个有限长序列 长为L N M 1 所以 Chirp z变换为两个有限长序列的线性卷积g k h k 可用圆圈卷积通过FFT来实现 h n 的主值序列可由h n 作周期延拓后取0 n L 1部分值获得 将与g n 作圆周卷积后 其输出的前M个值就是Chirp z变换的M个值 这个圆周卷积的过程可在频域上通过FFT实现 Chirp z变换的计算步骤 4 G k FFT g n L点 5 Y k G k H k L点 6 y n IFFT Y k L点 7 0 k M 1 2 用FFT求的傅里叶变换 H k FFT L点 3 对x n 加权并补零 1 求h n 的主值序列 利用FFT计算Chirp z变换 计计算量估算 乘 乘法数估计 1 2 两步可以事先计算 不必实时计算 3 7 两步两次加权共计N M次复乘 形成Y k 需L次复乘 一个FFT与一个IFFT需Llog2L次乘 所以总乘法数 L N M Llog2L 直接计算乘法数 NMN及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的运算方法 IDFT DFT 以上所讨论的FFT算法可用于IDFT运算 简称为IFFT比较IDFT的定义式 IDFT与DFT的差别 1 把DFT中的每一个系数改为 2 再乘以常数1 N 则以上所讨论的时间抽取或频率抽取的FFT运算均可直接用于IDFT运算 当然 蝶形中的系数应改为 b 第二种方法 完全不需要改动FFT程序 而是直接利用它作IFFT 考虑到故IFFT计算分三步 将X k 取共轭 虚部乘以 1 对直接作FFT 对FFT的结果取共轭并乘以1 N 得x n 2 4FFT应用中的几个问题 2 实数序列的FFT以上讨论的FFT算法都是复数运算 包括序列x n 也认为是复数 但大多数场合 信号是实数序列 任何实数都可看成虚部为零的复数 例如 求某实信号x n 的复谱 可认为是将实信号加上数值为零的虚部变成复信号 x n j0 再用FFT求其离散傅里叶变换 这种作法很不经济