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

离散傅里叶变换 DFT 及其快速算法 DFT的定义 DFT的主要性质 频域采样 快速傅里叶变换 FFT FFT应用 图4 1各种形式的傅里叶变换 几种形式的傅里叶变换 四种傅里叶变换形式的归纳 求有限长序列x n 的DFT的实质是 将有限长序列x n 作周期延拓x n N 求其DFS 取其主值序列 即可得X K 4 1离散傅里叶变换的定义 DFT的引出 例 已知序列x n n 求它的N点DFT 解单位脉冲序列的DFT很容易由DFT的定义式得到 k 0 1 N 1 n 的X k 如图所示 这是一个很特殊的例子 它表明对序列 n 来说 不论对它进行多少点的DFT 所得结果都是一个离散矩形序列 序列 n 及其离散傅里叶变换 例3 5已知x n cos n 6 R12 n 是一个长度N 12的有限长序列 求它的N点DFT 解由DFT的定义式 利用复正弦序列的正交特性式 再考虑到k的取值区间 可得 图3 10有限长序列及其DFT N 16 N 12 Nfft 16 为什么 设x n 为一长度为N的离散序列 其Z变换 DFT和DTFT分别为 则 若令 4 2DFT与序列傅里叶变换 Z变换的关系 表示Z平面单位圆上幅角为的点 对序列的Z变换在单位圆上取值 可得到该序列的傅立叶变换FT 即 若令Z ejw则可得所以有限长序列x n 在Z平面单位圆上的Z变换是该序列x n 的傅立叶变换 也即将Z平面单位圆N等分后的第k点 如图所示 所以x n 的DFTX K 等于它的Z变换X Z 在Z平面单位圆上N个等分点上的采样值 因此 DFT与序列傅里叶变换的关系为 即 X K 也是在w 2 k N处的采样值 如图所示 X k 与X ej 的关系 DFT DFT矩阵表示 DFT矩阵形式为 其中 DFT IDFT矩阵形式为 dftmtx N 函数产生N N的DFT矩阵DN conj dftmtx N N函数产生N N的IDFT矩阵DN 1 DFT 利用MATLAB计算DFT fft x fft x N ifft x ifft x N fft x 计算M点的DFT M是序列x的长度 fft x N 计算N点的DFT M N 将原序列裁为N点计算N点的DFT M N 将原序列补零至N点 然后计算N点DFT 设 x1 n 和x2 n 是两个有限长序列 长度均为N 4 3离散傅立叶变换的基本性质 X1 k DFT x1 n X2 k DFT x2 n 一 线性性质 其各自的离散付里叶变换分别为 式中 a b为任意常数 则 二 圆周移位性质 其过程为 1 序列的圆周移位 x n 的圆周移位定义为y n x n m NRN n 1 将x n 以N为周期进行周期延拓得x n N2 将x n N左移m位 得x n m N3 取其主值序列x n m NRN n 循环移位过程如图所示 循环移位过程 y n x n m NRN n 2 时域圆周移位定理 3 频域循环移位定理 y n 也是一个长度为N的序列 记为 三 圆周卷积定理 1 圆周卷积的定义 2 时域循环卷积定理 3 频域循环卷积定理 当N为偶数时 将上式中的n换成N 2 n可得到 四 DFT的共轭对称性 1 圆周共轭对称序列和共轭反对称序列 设xep n xop n 均为有限长序列 若xep n x ep N n 0 n N 1 则称xep n 为圆周共轭对称序列 若xop n x op N n 0 n N 1 则称xop n 为共轭反对称序列 如图所示 图中 表示对应点为序列取共轭后的值 圆周共轭对称与共轭反对称序列示意图 任何有限长序列x n 都可以表示成其共轭对称分量和共轭反对称分量之和 即 其对应的DFT 则 2 DFT的共轭对称性 如果 即XR k XR N k XI k XI N k 即 X k 的幅频特性是偶对称的 相频特性是奇对称的 若x n 是长度为N的实序列 且X k DFT x n 则 1 X k 圆周共轭对称 X k X N k 0 k N 1 利用DFT的共轭对称性 还可通过计算一个N点DFT 得到两个不同实序列的N点DFT 对于实序列 计算N点DFT 若N 偶 则只要计算前N 2 1点的DFT 其它点由X k X N k 可得到 若N 奇 则只要计算前 N 1 2点的DFT即可 所以X1 k DFT x1 n 1 2 X k X N k X2 k DFT x2 n j1 2 X k X N k 设x1 n 和x2 n 为两个N点的实序列 构成新序列x n x n x1 n jx2 n 对x n 进行DFT X k DFT x n Xep k Xop k Xep k DFT x1 n 1 2 X k X N k Xop k DFT jx2 n 1 2 X k X N k 例 已知一9点实序列的DFT在偶数点的值为X 0 3 1 X 2 2 5 4 6j X 4 1 7 5 2j X 6 9 3 6 3j X 8 5 5 8 0j 确定DFT在奇数点的值 解 X 1 X 9 1 X 8 5 5 8 0j X 3 X 9 3 X 6 9 3 6 3j X 5 X 9 5 X 4 1 7 5 2j X 7 X 9 7 X 2 2 5 4 6j 根据实序列DFT的对称特性X k X N k 可得 可见 当输入信号的频率为q 0时 X K 的N个值中只有X q N 其余皆为零 如果输入信号为若干个不同频率的信号的组合 经离散付里叶变换后 不同的k上 X k 将有一一对应的输出 因此 离散付里叶变换实质上对频率具有选择性 3 4 5选频性 设有复序列x n 0 n N 1 其离散付里叶变换为 其中q为整数 当 0 2 N时 对X k 求IDFT 得xN n IDFT X k 问题 由频域离散采样能否恢复原来的信号 其条件是什么 4 4频率域采样 1 频率域采样定理 设 任意序列x n 的Z变换为 设X z 收敛域包含单位圆 即x n 存在傅里叶变换 对X Z 在单位圆上等间隔采样N点 则得 因为X k 是xN n 以N为周期周期延拓后序列的离散傅里叶级数系数的主值序列 即 将式 3 5 1 代入上式得 求xN n 与x n 之间的关系 说明 xN n 为原序列x n 以N为周期周期延拓后的主值序列 式中 所以 若x n 为M长 则只有当频域采样点数N M时 才有即可由频域采样X K 恢复原序列x n 否则产生时域混叠现象 这就是所谓的频域采样定理 可看出 时域取样时 时域离散 频域周期 同样 频域取样时 频域离散 时域周期 时域抽样定理和频域抽样定理为利用数字化方式 分析和处理信号奠定了理论基础 例 已知序列对x n 的Z变换X z 在单位圆上等间隔采样N点 采样值为 求有限长序列IDFT X k 解 例 设序列x n 8 7 6 5 4 3 2 1 现对x n 的DTFT在一个周期内作N 6点的均匀抽样 得XN K 求XN K 的IDFTxN n 例 已知有限序列x n 1 1 4 3 n 0 1 2 3 序列x n 的DTFT为X ejw 记X ejw 在 w 2pk 3 k 0 1 2 的取样值为X k 求IDFT X k x n x n 3 R3 n 2 1 4 n 0 1 2 解 X ejw 在频域的离散化导致对应的时域序列x n 的周期化 IDFT X k 问题 由X K 能否完整地表达X Z 2 内插公式 设序列x n 长度为M N M 则有 X Z 的内插公式 它说明已知X K 可根据内插公式求任意Z点的X Z 所以X Z 的N个取样值X K 包含了X Z 的全部信息 令 内插函数 称为内插函数 令其分子为零 得 r 0 1 k N 1 即内插函数在单位圆的N等分点上 也即采样点上 有N个零点 而分母为零 则有z WN k 的一个极点 它将和第k个零点相抵消 因而 插值函数 k z 只在本身采样点r k处不为零 在其他 N 1 个采样点r上 r 0 1 N 1 但r k 都是零点 有 N 1 个零点 而它在z 0处还有 N 1 阶极点 进一步化简可得 当z ej 时 内插函数幅度特性与相位特性 N 5 当变量 0时 1 当 i 1 2 N 1 时 0 因而可知 满足以下关系 也就是说 函数在本采样点 而在其他采样点上 函数整个X ej 就是由N个函数分别乘上X k 后求和 所以 在每个采样点上X ej 就精确地等于X k 因为其他点的插值函数在这一点上的值为零 没有影响 各采样点之间的X ej 值由各采样点的加权插值函数在所求 点上的值的叠加得到的 实际应用中 大多数是求解线性卷积 如信号x n 通过系统h n 其输出就是线性卷积y n x n h n 圆周卷积 DFT乘积 线性卷积 FT乘积 而 圆周卷积比起线性卷积 在运算速度上有很大的优越性 它可以采用离散傅立叶变换的乘积实现 而快速付里叶变换 FFT 技术又大大提高了离散傅立叶变换的计算速度 所以利用圆周卷积求线性卷积 4 5用DFT计算线性卷积 设 x n h n 为两个有限长序列 其长度分别为N1和N2 线性卷积的定义为 而循环卷积的定义 要求L max N1 N2 即yc n 的长度为L 可看出y n y n 循环卷积与线性卷积的关系 即yl n 的长度为N1 N2 1 由此可见 循环卷积既可在时域直接计算 也可以在频域计算 由于DFT有快速算法FFT 当N很大时 在频域计算的速度快得多 因而常用DFT FFT 计算循环卷积 若 则由时域循环卷积定理有 Yc k DFT yc n X k H k 循环卷积与线性卷积的关系 所以 循环卷积yc n 等于线性卷积yl n 以L为周期的周期延拓序列的主值序列 因为yl n 的长度为N1 N2 1 因此只有L N1 N2 1时 yl n 以L为周期进行周期延拓才不会产生混叠 yc n yl n 将x n h n 的长度都变为N1 N2 1 则 求x n 与h n 的循环卷积x n h n 即为x n 与h n 的线性卷积x n h n 所以可将求线性卷积的问题转化为求循环卷积的问题 从而转化为求DFT 具体步骤为 1 将x n h n 均延长至N1 N2 1 变为x n h n 2 求x n h n 的DFT X K H K 3 求X K H K 的IDFT y n 则y n x n h n 用DFT求线性卷积 信号x n 通过数字系统h n 后得到输出y n x n h n 有时x n 很长 即长序列 4 6长序列线性卷积的计算 直接用DFT计算的缺点 1 信号要全部输入后才能进行计算 延迟太多 2 内存要求大 3 算法效率不高 需要将x n 分成若干段短序列 设h n 的长度为N x n 为长度为L的长序列 将x n 均匀分成长度为M的若干段 表示为 重叠相加法 可看出 每一项的长度为N M 1 大于序列分段的长度M 所以卷积结果求和时 每段卷积最后N 1个点和下一段的最前面的N 1个点重合 重合相加 可得到最后结果 如图所示 重叠相加法卷积示意图 4 7模拟信号的频谱分析 频谱分析过程 时域抽样 时域加窗 频域抽样 用DFT对连续时间信号进行频谱分析必然是近似的 其近似程度与信号带宽 抽样频率和窗函数宽度等有关 近似公式为 频谱分析存在的问题 混叠现象 栅栏效应 频谱泄漏 混叠现象 时域离散 频域周期 混叠 要求 fs 2fmax 频域抽样 F fs N 1 NTs 1 tp或N fs FF与信号的实际长度有关 指在频域中能辨认的频率 即频域取样中两相邻点间的频率间隔或基波频率 F愈小 频率分辨率愈高 截断效应 x n 可能无限长 将x n 截短 相当于y n x n RN n 则 其中 例如 x n cos 0n 0 4 其频谱为 其频谱图如图所示 所以 截断后序列的频谱与原序列的频谱有差别 1 泄漏现象 截断后 使原来的离散谱线向附近展宽 通常称这种展宽为泄漏 它使频谱模糊 谱分辨率下降 2 谱间干扰 在主谱线两边形成很多旁瓣 引起不同频率间的干扰 影响频谱分辨率 上述两种现象又称为截断效应 截断效应 例 试利用DFT分析一连续信号 已知其最高频率 1000Hz 要求频率分辨率F 2Hz DFT的点数必须为2的整数幂次 确定以下参数 最大的抽样间隔 最少的信号记录时间 最少的DFT点数 解 1 最大的抽样间隔Tmax为 2 最少的信号持续时间Tpmin为 3 由最大的抽样间隔Tmax与最少信号持续时间Tpmin 可得最少DFT点数N为 选择DFT的点数N 1024 以满足其为2的整数幂次 F 例 已知语音信号x t 的最高频率为fm 3 4kHz 用fsam 8kHz对x t 进行抽样 如对抽样信号做N 1600点的DFT 试确定X k 中k 600和k 1200点所分别对应原连续信号的连续频谱点f1和f2 kHz 解 对连续信号x t 按fsam 8kHz进行取样 得到对应的离散序列x n 在利用离散序列x n 的DFTX k 分析连续信号x t 的频谱时 X k 与X jw 存在以下对应关系 当k 600时 由于0 k N 2 1 所以 当k 1200时 由于N 2 k N 所以 N点DFT是在频率区间 0 2 上对信号频谱进行N点等间隔采样 得到的是若干个离散的频谱点X k 且它们限制在基频的整数倍上 这就好像在栅栏的一边通过缝隙看另一边的景象一样 只能在离散点处看到真实的景象 其余部分频谱成分被遮挡 所以称之为栅栏效应 减小栅栏效应方法 尾部补零 使谱线变密 增加频域采样点数 原来漏掉的某些频谱分量就可能被检测出来 栅栏效应 用DFT FFT 对周期信号进行谱分析 假设由模拟信号采样得到周期序列 其周期为N 对进行FT 得到 截取一个周期 即N长 对进行N点DFT 可用表示的频谱结构 如果截取长度M mN m为整数 则 令 令 而 如 x n cos n 6 N 12时的DFT N 24时的DFT N 16时的DFT 4 8快速傅里叶变换 FFT 频谱分析在数字信号处理中用途广泛 如通过语言信号的频谱分析实现语音通讯的频带压缩 声纳信号的频谱分析用以区分水面与水下目标 在各种测量仪器中 频谱分析用得更多 这些都需要DFT运算 虽然频谱分析和DFT运算很重要 但在很长一段时间里 由于DFT运算复杂 并没有得到真正的运用 而频谱分析仍大多采用模拟信号滤波的方法解决 直到1965年首次提出DFT运算的一种快速算法以后 情况才发生了根本 变化 人们开始认识到DFT运算的一些内在规律 从而很快地发展和完善了一套高速有效的运算方法 快速付里变换 FFT 算法 FFT的出现 使DFT的运算大大简化 运算时间缩短一 二个数量级 使DFT的运算在实际中得到广泛应用 DFT是FFT的理论基础 FFT只是DFT的一种快速高效算法 1 直接计算DFT的特点设x n 长度为N 其DFT为 X K 长度也为N 一 DFT运算的特点 可看出正 反变换形式上很相似 所以只讨论正变换 运算时间 设 x n 为复数序列 计算一个X K 值需要N次复数乘及N 1次复数加 则计算一个X K 序列需要N2次复数乘及N N 1 次复数加 所以直接计算DFT的计算量与N的平方成正比 当N较大时 计算量太大 2 提高DFT运算效率的依据 具有周期性和对称性 其周期性表现为 应用了上述性质后 独立的值只有N 2个 为简化DFT的运算提供了有力的依据 其圆周共轭对称性表现为 或写为 此外 二 基2FFT算法 基2时域抽取 Decimationintime FFT算法 基2频率抽取 Decimationinfrequency FFT算法 k 1基2时域抽取FFT算法 DIT FFT 设 x n 为一长度为N的序列 M为正整数按n的奇偶把x n 分解为两个N 2点的子序列 由于 所以 则x n 的DFT为 所以一个N点的DFT可分解为两个N 2点的DFT 由于X1 k 和X2 k 均以N 2为周期 且 而X k 为N点所以 其中X1 k 和X2 k 分别为x1 r 和x2 r 的N 2点DFT 即 用信号流图表示 X k 可表示为 呈蝶形 称为蝶形计算结构 是FFT运算中的一个基本单元 可将上述分解过程用计算流图来表示 N点DFT的一次时域抽取分解图 N 8 通过上述分解后 每个N 2点DFT只需要 N 2 2 N2 4次复数相乘 两个N 2点的DFT需要2 N 2 2 N2 2次复数乘 可见 分解后运算量大约节省了一倍 那么 X1 k 又可表示为 与第一次分解相同 将x1 r 按r的奇偶分解成两个N 4长的子序列x3 l 和x4 l 即 式中 分解图如图所示 将一个N 2点的DFT分解为两个N 4点的DFT 用同样的方法将x2 r 分解成两个N 4长的子序列x5 l 和x6 l 那么 X2 k 可表示为 其中 N点DFT的第二次时域抽取分解图 N 8 一个8点DFT就可分解为四个2点的DFT 如图所示 两点DFT 若N 8 则 N点DIT FFT运算流图 N 8 所以 N 8点的DIT FFT运算流图如图所示 DIT FFT算法与直接计算DFT运算量的比较从图可看出当N 2M时 要经过M级蝶算 每一级蝶算包含N 2个蝶形运算 所以总共需要的蝶形运算为 每个蝶形运算需要一次复数乘和两次复数加法 所以N点的FFT需要 个复数乘 个复数加 例如 N 210 1024时 DFT与FFT复数乘法运算之比为 其中 算法流图规律1 原位计算在FFT信号流图中 每一级里节点是成对出现的 如图所示 可看出第m 1级的两节点xm 1 p xm 1 q 经蝶算后 在m级中的节点序号是不变的 即xm p xm q 且蝶算自成独立单元 即xm p xm q 只与xm 1 p xm 1 q 有关 而与其它节点值无关 同时xm 1 p xm 1 q 也不再参加其它的蝶算 所以可以将计算出来的结果存入原来的单元 而不需要开辟新的地址存放计算结果 这种运算称为原位运算 可以大大的节省存储单元 减少成本 其中 p为对偶节点中上节点的序号q为对偶节点中下节点的序号在对偶节点中 只有下节点Xm q 才乘加权因子 2 旋转因子的变化规律 同一级中参加蝶算的两节点间的间距是相同的 对偶节点的间距可用下式确定 m为该级的序号 r可由下式确定 若 则数据不必交换 如 010若 则须将x n 与x 对调 为了避免重复对调 需要检查是否小于n 若n 说明x n 已与x 调换过 只有当n 时 才将x n 与x 的存放单元对换 如图所示 3 码位整序通常将x n 按顺序输入 通过变址运算 得到它的倒序排列 整序过程为 设N 8 若输入序列x n 的序号n用二进制数n2n1n0表示 则其反序二进制数表示为n0n1n2即 001 100 2 基2频域抽取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 时 当k取奇数 k 2r 1 r 0 1 N 2 1 时 令 则 所以在频域按K将一个N点的DFT 分成偶 奇两部分的N 2点的DFT 图4 2 4DIF FFT一次分解运算流图 N 8 同样还可继续分解 将N 2点的DFT分解成两个N 4点的DFT图4 2 5DIF FFT二次分解运算流图 N 8 图4 2 6DIF FFT运算流图 N 8 图4 2 7DIT FFT的一种变形运算流图 图4 2 8DIT FFT的一种变形运算流图 3 IDFT的高效算法上述FFT算法流图也可以用于离散傅里叶逆变换 InverseDiscreteFourierTransform 简称IDFT 比较DFT和IDFT的运算公式 对上式求共轭 所以由X K 求x n 的IFFT算法如下 1 求X K 的共轭X K 2 求X K 的FFT 得N x n 3 取x n 的共轭并除以N 即得x n 两边再同时取共轭 得 例1 已知序列x n cos 0 48 n cos 0 52 n 1 当0 n 10时 绘制x n 及它的DFT X K 2 0 n 10 后面填90个零时 绘制x n 及它的DFT X K 3 当0 n 100时 绘制x n 及它的DFT X K 4 9Matlab应用 Matlab程序 n 0 1 99 x cos 0 48 pi n cos 0 52 pi n n1 0 1 9 y1 x 1 1 10 0 n 10 后面填90个零时的x n 及 X K y2 x 1 1 10 zeros 1 90 n2 n subplot 323 stem n2 y2 b title x n N 10 补90个零 grid Y2 fft y2 100 magY2 abs Y2 w2 2 pi 100 n2 subplot 324 stem w2 pi magY2 b title DFT N 10 补90个零 grid 0 n 10时 x n 及它的DFT X K subplot 321 stem n1 y1 title x n N 10 grid Y1 fft y1 10 magY1 abs Y1 w1 2 pi 10 n1 subplot 322 stem w1 pi magY1 title DFT N 10 grid 0 n 100时 x n 及它的DFT X K