离散傅里叶变换ppt课件

第3章离散傅里叶变换 DFT 3 1离散傅里叶变换的定义及物理意义3 2离散傅里叶变换的基本性质3 3频率域采样3 4DFT的应用举例习题与上机题 傅里叶变换和Z变换是数字信号处理中常用的重要数学变换 对于有限长序列 还有一种更为重要的数学变换 即本章要讨论的离散傅里叶变换 DiscreteFourierTransform DFT DFT之所以更为重要 是因为其实质是有限长序列傅里叶变换的有限点离散采样 从而实现了频域离散化 使数字信号处理可以在频域采用数值运算的方法进行 这样就大大增加了数字信号处理的灵活性 更重要的是 DFT有多种快速算法 统称为快速傅里叶变换 FastFourierTransform FFT 从而使信号的实时处理和设备的简化得以实现 因此 时域离散系统的研究与应用在许多方面取代了传统的连续时间系统 所以说 DFT不仅在理论上有重要意义 而且在各种信号的处理中亦起着核心作用 本章主要讨论DFT的定义 物理意义 基本性质以及频域采样和DFT的应用举例等内容 3 1离散傅里叶变换的定义及物理意义 3 1 1DFT的定义 设x n 是一个长度为M的有限长序列 则定义x n 的N点离散傅里叶变换为 3 1 1 X k 的离散傅里叶逆变换 InverseDiscreteFourierTransform IDFT 为 式中 N称为DFT变换区间长度 N M 通常称 3 1 1 式和 3 1 2 式为离散傅里叶变换对 为了叙述简洁 常常用DFT x n N和IDFT X k N分别表示N点离散傅里叶变换和N点离散傅里叶逆变换 下面证明IDFT X k 的唯一性 3 1 2 把 3 1 1 式代入 3 1 2 式 有 由于 所以 在变换区间上满足下式 IDFT X k N x n 0 n N 1 由此可见 3 1 2 式定义的离散傅里叶逆变换是唯一的 例3 1 1 x n R4 n 求x n 的4点和8点DFT 解设变换区间N 4 则 设变换区间N 8 则 由此例可见 x n 的离散傅里叶变换结果与变换区间长度N的取值有关 对DFT与Z变换和傅里叶变换的关系及DFT的物理意义进行讨论后 上述问题就会得到解释 3 1 2DFT与傅里叶变换和Z变换的关系 设序列x n 的长度为M 其Z变换和N N M 点DFT分别为 比较上面二式可得关系式 3 1 3 或 3 1 4 3 1 3 式表明序列x n 的N点DFT是x n 的Z变换在单位圆上的N点等间隔采样 3 1 4 式则说明X k 为x n 的傅里叶变换X ej 在区间 0 2 上的N点等间隔采样 这就是DFT的物理意义 由此可见 DFT的变换区间长度N不同 表示对X ej 在区间 0 2 上的采样间隔和采样点数不同 所以DFT的变换结果不同 上例中 x n R4 n DFT变换区间长度N分别取8 16时 X ej 和X k 的幅频特性曲线图如图3 1 1所示 由此容易得到x n R4 n 的4点DFT为X k DFT x n 4 4 k 这一特殊的结果在下面将得到进一步解释 图3 1 1R4 n 的FT和DFT的幅度特性关系 3 1 3DFT的隐含周期性 前面定义的DFT变换对中 x n 与X k 均为有限长序列 但由于的周期性 使 3 1 1 和 3 1 2 式中的X k 隐含周期性 且周期均为N 对任意整数m 总有 所以 3 1 1 式中 X k 满足 实际上 任何周期为N的周期序列都可以看做长度为N的有限长序列x n 的周期延拓序列 而x n 则是的一个周期 即 为了以后叙述简洁 当N大于等于序列x n 的长度时 将 3 1 5 式用如下形式表示 3 1 7 式中x n N表示x n 以N为周期的周期延拓序列 n N表示模N对n求余 即如果n MN n10 n1 N 1 M为整数 则 n N n1 例如 则有 所得结果符合图3 1 2 a 和 b 所示的周期延拓规律 图3 1 2x n 及其周期延拓序列 现在解释DFT R4 n 4 4 k 根据DFT第二种物理解释可知 DFT R4 n 4表示R4 n 以4为周期的周期延拓序列R4 n 4的频谱特性 因为R4 n 4是一个直流序列 只有直流成分 即零频率成分 3 1 4用MATLAB计算序列的DFTMATLAB提供了用快速傅里叶变换算法FFT 算法见第4章介绍 计算DFT的函数fft 其调用格式如下 Xk fft xn N 调用参数xn为被变换的时域序列向量 N是DFT变换区间长度 当N大于xn的长度时 fft函数自动在xn后面补零 函数返回xn的N点DFT变换结果向量Xk 当N小于xn的长度时 fft函数计算xn的前面N个元素构成的N长序列的N点DFT 忽略xn后面的 元素 Ifft函数计算IDFT 其调用格式与fft函数相同 可参考help文件 例3 1 2 设x n R4 n X ej FT x n 分别计算X ej 在频率区间 0 2 上的16点和32点等间隔采样 并绘制X ej 采样的幅频特性图和相频特性图 解由DFT与傅里叶变换的关系知道 X ej 在频率区间 0 2 上的16点和32点等间隔采样 分别是x n 的16点和32点DFT 调用fft函数求解本例的程序ep312 m如下 例3 1 2程序ep312 m DFT的MATLB计算 xn 1111 输入时域序列向量xn R4 n Xk16 fft xn 16 计算xn的16点DFT Xk32 fft xn 32 计算xn的32点DFT 以下为绘图部分 省略 程序集中有 程序运行结果如图3 1 3所示 图3 1 3程序ep312 m运行结果 3 2离散傅里叶变换的基本性质 3 2 1线性性质如果x1 n 和x2 n 是两个有限长序列 长度分别为N1和N2 且 y n ax1 n bx2 n 式中 a b为常数 取N max N1 N2 则y n 的N点DFT为 Y k DFT y n N aX1 k bX2 k 0 k N 1 3 2 1 其中X1 k 和X2 k 分别为x1 n 和x2 n 的N点DFT 显然 y n 是长度为N的有限长序列 观察图3 2 1可见 循环移位的实质是将x n 左移m位 而移出主值区 0 n N 1 的序列值又依次从右侧进入主值区 循环移位 就是由此得名的 由循环移位的定义可知 对同一序列x n 和相同的位移m 当延拓周期N不同时 y n x n m NRn n 则不同 请读者画出N M 6 m 2时 x n 的循环移位序列y n 波形图 图3 2 1x n 及其循环移位过程 2 时域循环移位定理 设x n 是长度为M M N 的有限长序列 y n 为x n 的循环移位 即 则 3 2 3 其中 证明 令n m n 则有 由于上式中求和项以N为周期 因此对其在任一周期上的求和结果相同 将上式的求和区间改在主值区 则得 3 频域循环移位定理 如果X k DFT x n N0 k N 1 Y k X k l NRN k 则 3 2 4 3 2 4 式的证明方法与时域循环移位定理类似 直接对Y k X k l NRN k 进行IDFT即得证 3 2 3循环卷积定理 时域循环卷积定理是DFT中最重要的定理 具有很强的实用性 已知系统输入和系统的单位脉冲响应 计算系统的输出 以及FIR滤波器用FFT实现等 都是基于该定理的 下面首先介绍循环卷积的概念和计算循环卷积的方法 然后介绍循环卷积定理 1 两个有限长序列的循环卷积 设序列h n 和x n 的长度分别为N和M h n 与x n 的L点循环卷积定义为 3 2 5 式中 L称为循环卷积区间长度 L max N M 上式显然与第1章介绍的线性卷积不同 为了区别线性卷积 用表示循环卷积 用表示L点循环卷积 即yc n h n x n 观察 3 2 5 式 x n m L是以L为周期的周期信号 n和m的变化区间均是 0 L 1 因此直接计算该式比较麻烦 计算机中采用矩阵相乘或快速傅里叶变换 FFT 的方法计算循环卷积 下面介绍用矩阵计算循环卷积的公式 当n 0 1 2 L 1时 由x n 形成的序列为 x 0 x 1 x L 1 令n 0 m 0 1 L 1 由式 3 2 5 中x n m L形成x n 的循环倒相序列为 与序列x n 进行对比 相当于将第一个序列值x 0 不动 将后面的序列反转180 再放在x 0 的后面 这样形成的序列称为x n 的循环倒相序列 令n 1 m 0 1 L 1 由式 3 2 5 中x n m L形成的序列为 观察上式等号右端序列 它相当于x n 的循环倒相序列向右循环移一位 即向右移1位 移出区间 0 L 1 的序列值再从左边移进 再令n 2 m 0 1 L 1 此时得到的序列又是上面的序列向右循环移1位 依次类推 当n和m均从0变化到L 1时 得到一个关于x n m L的矩阵如下 3 2 6 上面矩阵称为x n 的L点 循环卷积矩阵 其特点是 1 第1行是序列 x 0 x 1 x L 1 的循环倒相序列 注意 如果x n 的长度M L 则需要在x n 末尾补L M个零后 再形成第一行的循环倒相序列 2 第1行以后的各行均是前一行向右循环移1位形成的 3 矩阵的各主对角线上的序列值均相等 有了上面介绍的循环卷积矩阵 就可以写出式 3 2 5 的矩阵形式如下 3 2 7 按照上式 可以在计算机上用矩阵相乘的方法计算两个序列的循环卷积 这里关键是先形成循环卷积矩阵 上式中如果h n 的长度N L 则需要在h n 末尾补L N个零 例3 2 1 计算下面给出的两个长度为4的序列h n 与x n 的4点和8点循环卷积 解按照式 3 2 21 写出h n 与x n 的4点循环卷积矩阵形式为 h n 与x n 的8点循环卷积矩阵形式为 h n 和x n 及其4点和8点循环卷积结果分别如图3 2 2 a b c 和 d 所示 请读者计算验证本例的8点循环卷积结果等于h n 与x n 的线性卷积结果 后面将证明 当循环卷积区间长度L大于等于y n h n x n 的长度时 循环卷积结果就等于线性卷积 图3 2 2序列及其循环卷积波形 3 2 8 证明直接对 3 2 8 式两边进行DFT 则有 令n m n 则有 因为上式中是以N为周期的 所以对其在任一个周期上求和的结果不变 因此 由于 因此 即循环卷积亦满足交换律 3 2 4复共轭序列的DFT 设x n 是x n 的复共轭序列 长度为N X k DFT x n N 则 且X N X 0 证明根据DFT的唯一性 只要证明 3 2 11 式右边等于左边即可 3 2 11 又由X k 的隐含周期性 有X N X 0 用同样的方法可以证明 3 2 12 3 2 5DFT的共轭对称性 1 有限长共轭对称序列和共轭反对称序列 为了区别于傅里叶变换中所定义的共轭对称 或共轭反对称 序列 下面用xep n 和xop n 分别表示有限长共轭对称序列和共轭反对称序列 则二者满足如下关系式 3 2 13a 3 2 13b 当N为偶数时 将上式中的n换成N 2 n 可得到 上式更清楚地说明了有限长序列共轭对称序列是关于n N 2点对称 容易证明 如同任何实函数都可以分解成偶对称分量和奇对称分量一样 任何有限长序列x n 都可以表示成其共轭对称分量和共轭反对称分量之和 即 2 DFT的共轭对称性 1 如果将x n 表示为 x n xr n jxi n 3 2 17 其中 那么 由 3 2 11 式和 3 2 16a 式可得 由DFT的线性性质即可得 3 2 20 其中 Xop k DFT xr n 是X k 的共轭对称分量 Xop k DFT jxi n 是X k 的共轭反对称分量 2 如果将x n 表示为 3 2 21 其中 是x n 的共轭对称分量 是x n 的共轭反对称分量 那么 由 3 2 12 式可得 因此 3 2 22 其中 综上所述 可总结出DFT的共轭对称性质 如果序列x n 的DFT为X k 则x n 的实部和虚部 包括j 的DFT分别为X k 的共轭对称分量和共轭反对称分量 而x n 的共轭对称分量和共轭反对称分量的DFT分别为X k 的实部和虚部乘以j 另外 请读者根据上述共轭对称性证明有限长实序列DFT的共轭对称性 见本章习题题7 设x n 是长度为N的实序列 且X k DFT x n N 则X k 满足如下对称性 1 X k 共轭对称 即X k X N k k 0 1 N 1 3 2 23 2 如果x n 是偶对称序列 即x n x N n 则X k 实偶对称 即X k X N k 3 2 24 3 如果是奇对称序列 即x n x N n 则X k 纯虚奇对称 即X k X N k 3 2 25 实际中经常需要对实序列进行DFT 利用上述对称性质 可减少DFT的运算量 提高运算效率 例如 计算实序列的N点DFT时 当N 偶数时 只需计算X k 的前面N 2 1点 而N 奇数时 只需计算X k 的前面 N 1 2点 其他点按照 3 2 23 式即可求得 例如 X N 1 X 1 X N 2 X 2 这样可以减少近一半运算量 例3 2 2 利用DFT的共轭对称性 设计一种高效算法 通过计算一个N点DFT 就可以计算出两个实序列x1 n 和x2 n 的N点DFT 解构造新序列x n x1 n jx2 n 对x n 进行DFT 得到 由 3 2 17 3 2 18 和 3 2 19 式得到 所以 由X k 可以求得两个实序列x1 n 和x2 n 的N点DFT 3 3频率域采样 时域采样定理告诉我们 在一定条件下 可以由时域离散采样信号恢复原来的连续信号 那么能不能也由频域离散采样恢复原来的信号 或原连续频率函数 其条件是什么 内插公式又是什么形式 本节就上述问题进行讨论 设任意序列x n 的Z变换为 且X z 的收敛域包含单位圆 即x n 存在傅里叶变换 在单位圆上对X z 等间隔采样N点 得到 显然 3 3 1 式表示在区间 0 2 上对x n 的傅里叶变换X ej 的N点等间隔采样 将X k 看做长度为N的有限长序列xN n 的DFT 即 下面推导序列xN n 与原序列x n 之间的关系 并导出频域采样定理 3 3 1 由DFT与DFS的关系可知 X k 是xN n 以N为周期的周期延拓序列的离散傅里叶级数系数的主值序列 即 因此 所以 3 3 3 3 3 2 式 3 3 3 说明 X z 在单位圆上的N点等间隔采样X k 的N点IDFT是原序列x n 以N为周期的周期延拓序列的主值序列 综上所述 可以总结出频域采样定理 如果序列x n 的长度为M 则只有当频域采样点数N M时 才有xN n IDFT X k x n 即可由频域采样X k 恢复原序列x n 否则产生时域混叠现象 满足频域采样定理时 频域采样序列X k 的N点IDFT是原序列x n 所以必然可以由X k 恢复X z 和X ej 下面推导用频域采样X k 表示X z 和X ej 的内插公式和内插函数 设序列x n 长度为M 在频域 0 2 上等间隔采样N点 N M 则有 式 3 3 6 称为用X k 表示X z 的内插公式 k z 称为内插函数 将z ej 代入 3 3 4a 式 并进行化简 可得 3 3 7 3 3 8 在数字滤波器的结构与设计中 我们将会看到 频域采样理论及有关公式可提供一种有用的滤波器结构和滤波器设计途径 3 3 7 式有助于分析FIR滤波器频率采样设计法的逼近性能 例3 3 1 长度为26的三角形序列x n 如图3 3 1 a 所示 编写MATLAB程序验证频域采样理论 解解题思想 先计算x n 的32点DFT 得到其频谱函数X ej 在频率区间 0 2 上等间隔32点采样X32 k 再对X32 k 隔点抽取 得到X ej 在频率区间 0 2 上等间隔16点采样X16 k 最后分别对X16 k 和X32 k 求IDFT 得到 绘制x16 n 和x32 n 波形图验证频域采样理论 MATLAB求解程序ep331 m如下 数字信号处理 第三版 第3章例3 3 1程序ep331 频域采样理论验证 M 26 N 32 n 0 M xa 0 M 2 xb ceil M 2 1 1 0 xn xa xb 产生M长三角波序列x n Xk fft xn 512 512点FFT x n X32k fft xn 32 32点FFT x n x32n ifft X32k 32点IFFT X32 k 得到x32 n X16k X32k 1 2 N 隔点抽取X32k得到X16 k x16n ifft X16k N 2 16点IFFT X16 k 得到x16 n 以下绘图部分省略 程序运行结果如图3 3 1所示 图3 3 1 a 和 b 分别为X ej 和x n 的波形 图3 3 1 c 和 d 分别为X ej 的16点采样 X16 k 和x16 n IDFT X16 k 16波形图 图3 3 1 e 和 f 分别为X ej 的32点采样 X32 k 和x32 n IDFT X32 k 32波形图 由于实序列的DFT满足共轭对称性 因此频域图仅画出 0 上的幅频特性波形 本例中x n 的长度M 26 从图中可以看出 当采样点数N 16M时 无时域混叠失真 x32 n IDFT X32 k x n 图3 3 1频域采样定理验证 则由DFT的时域循环卷积定理有 由此可见 循环卷积既可以在时域直接计算 也可以按照图3 4 1所示的计算框图在频域计算 由于DFT有快速算法 当L很大时 在频域计算循环卷积的速度快得多 因而常用DFT FFT 计算循环卷积 图3 4 1用DFT计算循环卷积的原理框图 在实际应用中 为了分析时域离散线性时不变系统或者对序列进行滤波处理等 需要计算两个序列的线性卷积 与计算循环卷积一样 为了提高运算速度 也希望用DFT FFT 计算线性卷积 而DFT只能直接用来计算循环卷积 因此 下面先导出线性卷积和循环卷积之间的关系以及循环卷积与线性卷积相等的条件 最后得出用图3 4 1线性卷积的条件 对照 3 4 1 式可以看出 上式中 即 3 4 3 图3 4 2线性卷积与循环卷积波形图 综上所述 取L N M 1 则可按照如图3 4 1所示的计算框图用DFT FFT 计算线性卷积 其中DFT和IDFT通常用快速算法 FFT 来实现 故常称其为快速卷积 实际上 经常遇到两个序列的长度相差很大的情况 例如M N 若仍选取L N M 1 以L为循环卷积区间 并用上述快速卷积法计算线性卷积 则要求对短序列补很多零点 而且长序列必须全部输入后才能进行快速计算 因此要求存储容量大 运算时间长 并使处理延时很大 不能实现实时处理 况且在某些应用场合 序列长度不定或者认为是无限长 如电话系统中的语音信号和地震检测信号等 显然 在要求实时处理时 直接套用上述方法是不行的 解决这个问题的方法是将长序列分段计算 这种分段处理方法有重叠相加法和重叠保留法两种 下面只介绍重叠相加法 重叠保留法作为本章习题题21 留给读者讨论 设序列h n 长度为N x n 为无限长序列 将x n 等长分段 每段长度取M 则 3 4 4a 于是 h n 与x n 的线性卷积可表示为 3 4 4b 式中 3 4 4b 式说明 计算h n 与x n 的线性卷积时 可先计算分段线性卷积yk n h n xk n 然后把分段卷积结果叠加起来即可 如图3 4 3所示 每一分段卷积yk n 的长度为N M 1 因此相邻分段卷积yk n 与yk 1 n 有N 1个点重叠 必须把重叠部分的yk n 与yk 1 n 相加 才能得到正确的卷积序列y n 显然 可用图3 4 1所示的快速卷积法计算分段卷积yk n 其中L N M 1 由图3 4 3可以看出 当第二个分段卷积y1 n 计算完后 叠加重叠点便可得输出序列y n 的前2M个值 同样道理 分段卷积yi n 计算完后 就可得到y n 第i段的M个序列值 因此 这种方法不要求大的存储容量 且运算量和延时也大大减少 最大延时TDmax 2MTs To Ts是系统采样间隔 To是计算1个分段卷积所需时间 一般要求To MTs 这样 就实现了边输入边计算边输出 如果计算机的运算速度快 可以实现实时处理 图3 4 3用重叠相加法计算线性卷积时域关系示意图 用DFT计算分段卷积yk n 的方法如下 1 i 0 L N M 1 计算并保存H k DFT h n L 2 读入xk n x n RM n kM 构造变换区间 0 L 1 上的序列 实际中就是将xi n 的M个值存放在长度为M的数组中 并计算 3 4 n 0 1 2 L 1 5 计算 6 i i 1 返回 2 应当说明 一般x n 是因果序列 假设初始条件y 1 n 0 MATLAB信号处理工具箱中提供了一个函数fftfilt 该函数用重叠相加法实现线性卷积的计算 调用格式为 y fftfilt h x M 式中 h是系统单位脉冲响应向量 x是输入序列向量 y是系统的输出序列向量 h与x的卷积结果 M是由用户选择的输入序列x的分段长度 缺省M时 默认输入序列x的分段长度M 512 例3 4 1 假设h n R5 n x n cos n 10 cos 2 n 5 u n 用重叠相加法计算y n h n x n 并画出h n x n 和y n 的波形 解h n 的长度为N 5 对x n 进行分段 每段长度为M 10 计算h n 和x n 的线性卷积的MATLAB程序如下 例3 4 1重叠相加法的MATLAB实现程序 ep341 m Lx 41 N 5 M 10 Lx为信号序列x n 长度 hn ones 1 N hn1 hnzeros 1 Lx N 产生h n 其后补零是为了绘图好看 n 0 L 1 xn cos pi n 10 cos 2 pi n 5 产生x n 的Lx个样值 yn fftfilt hn xn M 调用fftfilt用重叠相加法计算卷积 以下为绘图部分 省略 运行程序画出h n x n 和y n 的波形如图3 4 4所示 请读者从理论上证明y n 的稳态波形是单一频率的正弦波 运行绘图程序fig345 m可以得到用重叠相加法求解本例题的xk n yk n 和y n y0 n y1 n y2 n y3 n 如图3 4 5所示 图3 4 4例3 4 1的求解程序运行结果 图3 4 5重叠相加法时域波形 3 4 2用DFT对信号进行谱分析 1 用DFT对连续信号进行谱分析 工程实际中 经常遇到连续信号xa t 其频谱函数Xa j 也是连续函数 为了利用DFT对xa t 进行频谱分析 先对xa t 进行时域采样 得到x n xa nT 再对x n 进行DFT 得到的X k 则是x n 的傅里叶变换X ej 在频率区间 0 2 上的N点等间隔采样 这里x n 和X k 均为有限长序列 然而 由傅里叶变换理论知道 若信号持续时间有限长 则其频谱无限宽 若信号的频谱有限宽 则其持续时间必然为无限长 所以严格地讲 持续时间有限的带限信号是不存在的 因此 按采样定理采样时 上述两种情况下的采样序列x n xa nT 均应为无限长 不满足DFT的变换条件 实际上对频谱很宽的信号 为防止时域采样后产生频谱混叠失真 可用预滤波器滤除幅度较小的高频成分 使连续信号的带宽小于折叠频率 对于持续时间很长的信号 采样点数太多 以致无法存储和计算 只好截取有限点进行DFT 由上述可见 用DFT对连续信号进行频谱分析必然是近似的 其近似程度与信号带宽 采样频率和截取长度有关 实际上从工程角度看 滤除幅度很小的高频成分和截去幅度很小的部分时间信号是允许的 因此 在下面分析中 假设xa t 是经过预滤波和截取处理的有限长带限信号 设连续信号xa t 持续时间为Tp 最高频率为fc 如图3 4 6 a 所示 xa t 的傅里叶变换为 对xa t 以采样间隔T 2fc 1 即fs 1 T 2fc 采样得x n xa nT 设共采样N点 并对Xa jf 作零阶近似 t nT dt T 得 显然 Xa jf 仍是f的连续周期函数 x n 和X jf 如图3 4 5 b 所示 对X jf 在区间 0 Fs 上等间隔采样N点 采样间隔为F 如图3 4 5 c 所示 参数Fs Tp N和F满足如下关系式 由于NT Tp 所以 3 4 5 3 4 6 将f kF和式 3 4 5 代入X jf 中可得Xa jf 的采样 0 k N 1 令 则 3 4 8 3 4 7 同理 由 可以推出 图3 4 6用DFT分析连续信号谱的原理示意图 3 4 7 式说明 可以通过对连续信号采样并进行DFT再乘以T 近似得到模拟信号频谱的周期延拓函数在第一个周期 0 fs 上的N点等间隔采样 如图3 4 6所示 对满足假设的持续时间有限的带限信号 在满足时域采样定理时 包含了模拟信号频谱的全部信息 k 0 1 2 N 2 表示正频率频谱采样 k N 2 1 N 2 2 N 1 表示负频率频谱采样 值得注意 如果xa t 持续时间无限长 上述分析中要进行截断处理 所以会产生所谓的截断效应 从而使谱分析产生误差 本节最后将讨论上述误差问题产生的原因及改进措施 下面举例说明截断效应 理想低通滤波器的单位冲激响应ha t 及其频响函数Ha f 如图3 4 7 a b 所示 图3 4 7 a 中只画出ha t 所截取的一段 图中 Ha kF 如图3 4 7 c 中黑点所示 由图可见 低频部分近似理想低通频响特性 而高频误差较大 且整个频响都有波动 这些误差就是由于对ha t 截断所产生的 所以通常称之为截断效应 为减少这种截断误差 可适当加长Tp 增加采样点数N或用窗函数处理后再进行DFT 有关窗函数的内容将在FIR数字滤波器设计中详细叙述 图3 4 7用DFT计算理想低通滤波器的频响曲线 例3 4 2 对实信号进行谱分析 要求谱分辨率F 10Hz 信号最高频率fc 2 5kHz 试确定最小记录时间Tpmin 最大的采样间隔Tmax 最少的采样点数Nmin 如果fc不变 要求谱分辨率提高1倍 最少的采样点数和最小的记录时间是多少 解 因此Tpmin 0 1s 因为要求Fs 2fc 所以 为使用DFT的快速算法FFT 希望N符合2的整数幂 为此选用N 512点 为使频率分辨率提高1倍 即F 5Hz 要求 用快速算法FFT计算时 选用N 1024点 上面分析了为提高谱分辨率 又保持谱分析范围不变 必须增长记录时间Tp 增加采样点数 应当注意 这种提高谱分辨率的条件是必须满足时域采样定理 甚至采样速率Fs取得更高 2 用DFT对序列进行谱分析 我们知道单位圆上的Z变换就是序列的傅里叶变换 即 X ej 是 的连续周期函数 如果对序列x n 进行N点DFT得到X k 则X k 是在区间 0 2 上对X ej 的N点等间隔采样 频谱分辨率就是采样间隔2 N 因此序列的傅里叶变换可利用DFT 即FFT 来计算 对周期为N的周期序列 由 2 3 10 式知道 其频谱函数为 其中 由于以N为周期 因而X ej 也是以2 为周期的离散谱 每个周期有N条谱线 第k条谱线位于 2 N k处 代表的k次谐波分量 而且 谱线的相对大小与成正比 由此可见 周期序列的频谱结构可用其离散傅里叶级数系数表示 由 DFT 的隐含周期性知道 截取的主值序列 并进行N点DFT 得到 所以可用X k 表示的频谱结构 3 4 16 如果截取长度M等于的整数个周期 即M mN m为正整数 即 令n n iN i 0 1 m 1 n 0 1 N 1 则 3 4 17 因为 所以 3 4 18 由此可见 XM k 也能表示的频谱结构 只是在k im时 表示的i次谐波谱线 其幅度扩大m倍 而其他k值时 XM k 0 当然 X i 与XM im 对应点频率是相等的 所以 只要截取的整数个周期进行DFT 就可得到它的频谱结构 达到谱分析的目的 如果的周期预先不知道 可先截取M点进行DFT 即 再将截取长度扩大1倍 截取 比较XM k 和X2M k 如果二者的主谱差别满足分析误差要求 则以XM k 或X2M k 近似表示的频谱 否则 继续将截取长度加倍 直至前后两次分析所得主谱频率差别满足误差要求 设最后截取长度为iM 则XiM k0 表示 2 iM k0点的谱线强度 在很多实际应用中 并非整个单位圆上的频谱都很有意义 例如 对于窄带信号 往往只希望对信号所在的一段频带进行谱分析 这时便希望采样能密集地在这段频带内进行 而带外部分可完全不予考虑 另外 有时希望采样点不局限于单位圆上 例如 语音信号处理中 常常需要知道系统极点所对应的频率 如果极点位置离单位圆较远 则其单位圆上的频谱就很平滑 如图3 4 8 a 所示 这时很难从中识别出极点对应的频率 图3 4 8单位圆与非单位圆谱分析示意图 例如 要求计算序列在半径为r的圆上的频谱 那么N个等间隔采样点为 的频谱分量为 令 则 上式说明 要计算x n 在半径为r的圆上的N点等间隔频谱分量 可以先对x n 乘以r n 再计算N点DFT FFT 即可得到 若要求x n 分布在该圆的有限角度内的N点等间隔频谱采样 可以在尾部补M N个零 仍按 3 4 10 式用 DFT分析整个圆上的M点等间隔频谱 最后只取所需角度内的N点等间待间隔采样即可 显然 这种方法的计算量大 效率低 下面要介绍的Chirp Z变换可使这种谱分析的运算量大大减少 3 4 19 式中 A0和W0为实数 当k 0时 有 3 4 22 由此可见 3 4 20 式中 A决定谱分析起始点z0的位置 W0的值决定分析路径的盘旋趋势 0则表示两个相邻分析点之间的夹角 如果W01时 向内盘旋 如图3 4 9所示 图3 4 9 Chirp Z变换分析频点分布 将zk代入Z变换公式得到 利用下面的关系式 得到 令 3 4 23 3 4 24 则 3 4 25 3 4 25 式说明 长度为N的序列x n 的M点Chirp Z变换可通过预乘得到y n 再计算y n 与h n 的卷积v k 最后乘以这三个步骤得到 计算方框图如图3 4 10所示 图中 h n 看成一个数字网络的单位脉冲响应 输出v n y n h n n k V k 当W0 1时 是一个频率随时间线性增长的复指数序列 在雷达系统中 这样的信号称作线性调频信号 并用专用词汇Chirp表示 因此对上述变换起名为线性调频Z变换 简称Chirp Z变换 CZT 图3 4 10 Chirp Z变换的计算框图 图3 4 11 Chirp Z变换中hL n 的形成 而v n 的非零区间为 N 1 N M 2 为了用循环卷积代替线性卷积计算出v n 在 0 M 1 区间上的M个序列值 以L为周期进行周期延拓时 0 M 1 区间上不能有混叠 所以循环卷积区间长度L应大于或等于N M 1 一般用快速卷积法 FFT算法 计算 故应选择L N M 1 同时又满足L 2m m为自然数 的最小值 若选择L N M 1 那么y n 尾部应补M 1个零 并将h n 从 N 1 到M 1所截取的一段序列以L为周期进行周期延拓 取主值序列形成hL n 如图3 4 11 c 所示 L 这时可以用快速卷积法计算如上构造的两个序列y n 和hL n 的循环卷积 应当注意 当选择L 2m N M 1时 y n 应补L N个零点 而h n 应从M L到M 1区间上截取 或按上述区间 N 1 M 1截取后在 N 1点前面补L N M 1 个零点 后 以L为周期进行周期延拓 综上所述 可归纳出具体计算步骤如下 1 形成hL n 序列 2 3 4 Y k DFT y n 0 k L 1 5 计算Y k Y k 6 V k IDFT Y m H m 0 k L 1 7 与标准DFT FFT 算法相比较 Chirp Z变换有以下特点 1 输入序列长度N和输出序列长度不需要相等 且二者均可为素数 2 分析频率点zk的起始点z0及相邻两点的夹角 是任意的 即频率分辨率是任意的 因此可从任意频率上开始 对输入数据进行窄带高分辨率的谱分析 3 谱分析路径可以是螺旋形的 4 当 时 zk均匀分布在单位圆上 此时Chirp CD 2 Z变换就是序列的DFT 因此可以说 DFT是Chirp Z变换的特例 总之 Chirp Z变换用作谱分析时 具有灵活 适应性强和运算效率高等优点 4 用DFT进行谱分析的误差问题 DFT 实际中用FFT计算 可用来对连续信号和数字信号进行谱分析 在实际分析过程中 要对连续信号采样和截断 有些非时限数据序列也要截断 由此可能引起分析误差 下面分别对可能产生误差的三种现象进行讨论 1 混叠现象 对连续信号进行谱分析时 首先要对其采样 变成时域离散信号后才能用DFT FFT 进行谱分析 采样速率Fs必须满足采样定理 否则会在 对应模拟频率f Fs 2 附近发生频谱混叠现象 这时用DFT分析的结果必然在f Fs 2附近产生较大误差 因此 理论上必须满足Fs 2fc fc为连续信号的最高频率 对Fs确定的情况 一般在采样前进行预滤波 滤除高于折叠频率Fs 2的频率成分 以免发生频率混叠 现象 2 栅栏效应 我们知道 N点DFT是在频率区间 0 2 上对时域离散信号的频谱进行N点等间隔采样 而采样点之间的频谱函数是看不到的 这就好像从N个栅栏缝隙中观看信号的频谱情况 仅得到N个缝隙中看到的频谱函数值 因此称这种现象为栅栏效应 由于栅栏效应 有可能漏掉 挡住 大的频谱分量 为了把原来被 栅栏 挡住的频谱分量检测出来 对有限长序列 可以在原序列尾部补零 对无限长序列 可以增大截取长度及DFT变换区间长度 从而使频域采样间隔变小 增加频域采样点数和采样点位置 使原来漏掉的某些频谱分量被检测出来 对连续信号的谱分析 只要采样速率Fs足够高 且采样点数满足频率分辨率要求 见 3 4 14 式 就可以认为DFT后所得离散谱的包络近似代表原信号的频谱 3 截断效应 实际中遇到的序列x n 可能是无限长的 用DFT对其进行谱分析时 必须将其截短 形成有限长序列y n x n w n w n 称为窗函数 长度为N w n RN n 称为矩形窗函数 根据傅里叶变换的频域卷积定理 有 其中 对矩形窗数w n RN n 有 幅度谱Wg 曲线如图3 4 12所示 Wg 以2 为周期 只画低频部分 图中 2 N的部分称为主瓣 其余部分称为旁瓣 图3 4 12矩形窗的幅度谱 例如 x n cos 0n 0 4 其频谱为 x n 的频谱X ej 如图3 4 13 a 所示 将x n 截断后 y n x n RN n 的幅频曲线如图3 4 13 b 所示 图3 4 13x n cos 0n 加矩形窗前 后的幅频特性 由上述可见 截断后序列的频谱Y ej 与原序列频谱X ej 必然有差别 这种差别对谱分析的影响主要表现在如下两个方面 1 泄露 由图3 4 13 b 可知 原来序列x n 的频谱是离散谱线 经截断后 使原来的离散谱线向附近展宽 通常称这种展宽为泄露 显然 泄露使频谱变模糊 使谱分辨率降低 从图3 4 13可以看出 频谱泄露程度与窗函数幅度谱的主瓣宽度直接相关 在第7章将证明 在所有的窗函数中 矩形窗的主瓣是最窄的 但其旁瓣的幅度也最大 2 谱间干扰 在主谱线两边形成很多旁瓣 引起不同频率分量间的干扰 简称谱间干扰 特别是强信号谱的旁瓣可能湮没弱信号的主谱线 或者把强信号谱的旁瓣误认为是另一频率的信号的谱线 从而造成假信号 这样就会使谱分析产生较大偏差 由于上述两种影响是由对信号截断引起的 因此称之为截断效应 由图3 4 12可以看出 增加N可使Wg 的主瓣变窄 减小泄露 提高频率分辨率 但旁瓣的相对幅度并不减小 为了减小谱间干扰 应用其它形状的窗函数w n 代替矩形窗 窗函数将在FIR数字滤波其设计中介绍 但在N一定时 旁瓣幅度越小的窗函数 其主瓣就越宽 所以 在DFT变换区间 即截取长度 N一定时 只能以降低谱分析分辨率为代价 换取谱间干扰的减小 通过进一步学习数字信号处理的功率谱估计等现代谱估计内容可知 减小截断效应的最好方法是用近代谱估计的方法 但谱估计只适用于不需要相位信息的谱分析场合 最后要说明的是 栅栏效应与频率分辨率是不同的两个概念 如果截取长度为N的一段数据序列 则可以在其后面补N个零 再进行2N点DFT 使栅栏宽度减半 从而减轻了栅栏效应 但是这种截短后补零的方法不能提高频率分辨率 因为截短已经使频谱变模糊 补零后仅使采样间隔变小 但得到的频谱采样的包络仍是已经变模糊的频谱 所以频率分辨率没有提高 因此 要提高频率分辨率 就必须对原始信号截取的长度加长 对模拟信号 就是增加采样时间Tp的长度 习题与上机题1 计算以下序列的N点DFT 在变换区间0 n N 1内 序列定义为 1 x n 1 2 x n n 3 x n n n0 0 n0 N 4 x n Rm n 0 m N 5 6 7 8 x n sin 0n RN n 9 x n cos 0n RN N 10 x n nRN n 2 已知下列X k 求x n IDFT X k 1 2 其中 m为正整数 0 m N 2 3 已知长度为N 10的两个有限长序列 做图表示x1 n x2 n x1 n 与x2 n 的10点和20点循环卷积 4 证明DFT的对称定理 即假设X k DFT x n 证明 DFT X n Nx N k 5 如果X k DFT x n 证明DFT的初值定理 6 设x n 的长度为N 且X k DFT x n 0 k N 1 令h n x n NRmN n m为自然数 H k DFT h n mN 0 k mN 1 求H k 与X k 的关系式 10 证明离散相关定理 若 则 11 证明离散帕塞瓦尔定理 若X k D