离散系统分析和离散傅里叶变换

第四章 离散系统分析和离散傅里叶变换 4-1概述 在上一章中我们已经介绍了连续时间信号（周期的或非周期的）的傅里叶变换。在第一、二章中介绍了离散信号和离散系统的概念，在这一章中主要讨论离散信号的傅里叶变换。 4-2离散信号的傅里叶变换 时域抽样定理告诉我们，连续时间信号可以由它的样本值恢复出来，即 当抽样频率 给定时，抽样函数 就确定了，唯一与信号相关的是信号的样本值 ，换句话说传载 中信息的是样本值 。因此研究连续时间信号 中的信息，就转变为研究样本值 中的信息。当抽样频率 给定时， 也就一定了，样本值 就可以抽象为序列 ，也就是说离散信号的数学抽象是序列。以后我们就用序列 表示离散信号（样本值）。 由于序列的变量是整数变量，与连续信号的变量不同，因此对序列的处理方法与连续时间变量的处理方法也必定不同。先来看看序列的傅里叶变换，连续非周期时间信号 的傅里叶变换为 假定 是非周期的，仿照连续时间信号的傅里叶变换形式可以定义序列的傅里叶变换： （4-1） （4-2） 式中 为数字角频率。（4-1）式和（4-2）式构成了序列的傅里叶变换对，前者称为序列的傅里叶正变换，后者称为序列的傅里叶逆变换。注意到序列傅里叶正变换公式是个和式，这是因为序列 的变量是离散的整数，序列的傅里叶逆变换公式是个积分式，由此也说明序列的傅里叶变换是 的连续函数，也就是说，离散信号的傅里叶变换是频域中连续的函数。此外因 所以任何序列的傅里叶变换都是以 为周期的频域连续函数。 序列的傅里叶变换具有如下性质： 1. 线性特性 若    ，    则        （4-3） 式中a和b均为常数。 2. 时间位移特性 若    则        （4-4） 式中 为任意整数。 3. 频率位移特性 若    则        （4-5） 式中 为任意常数。 4. 对称特性 若 为实数序列，且有 则称 为偶序列（even sequence），通常用下标e表示偶序列，即 。 若 为实数序列，且    则称 为奇序列（odd sequence），通常用下标o表示奇序列，即 。 任何序列都可以表示为偶序列与奇序列之和，即 （4-6） 其中            （4-7） （4-8） 若 为复数序列，且其实部为偶对称，虚部为奇对称，即 则称此序列为共轭对称序列（conjugate symmetric sequence），通常表示为 。 若 为复数序列，且其实部为奇对称，虚部为偶对称，即 则称此序列为共轭反对称序列（conjugate ant symmetric sequence），通常表示为 。 任意复数序列 均可表示为共轭对称序列与共轭反对称序列之和，即 （4-9） 其中            （4-10） （4-11） 实际上，（4-9）式与（4-6）式是等价的，当 为实数序列时，（4-9）式就变成（4-6）式了。 若    则            （4-12） （4-13） （4-12）式说明共轭序列的傅里叶变换等于原序列傅里叶变换的共轭函数的反函数。（4-13）式说明反序列的共轭序列的傅里叶变换等于原序列傅里叶变换的共轭函数。这个性质再一次表明了时域和频域的对称性。 4-3周期序列的傅里叶级数（DFS） 我们知道一个周期为T的连续时间信号 可以展开成傅里叶级数，即 式中    ， 。 对于一个周期为N的离散信号 ，也可以展开成傅里叶级数，注意到连续时间信号展开成傅里叶级数是将信号表示成一系列基波频率 整倍数频率上复指数函数 的加权和。由此我们想到一个周期序列也展开成其基波频率 整倍数频率上复指数 的加权和。比较 和 这两个复值数函数表达式，可以看出有两点不同，一是连续时间信号 的周期T是个模拟量，而周期序列的周期N则为整数值；二是连续时间信号 的自变量t是连续时间变量，而离散时间信号 的自变量n是离散变量（整数值）正是因为存在着这种差别，决定了周期离散信号的傅里叶级数与周期连续信号的傅里叶级数有着本质的区别。在周期连续信号的傅里叶级数表达式中，复指数（谐波分量） 有无穷多个，这表现在傅里叶级数是n无穷和式。然而对于周期离散信号的复指数（谐波分量） 只有N个独立分量，这是因为 同理可以推导出 以上二式说明复指数 既是变量k的周期序列，也是变量n的周期序列，周期均为N。因此周期信号只能分解在独立的N个 分量上，即有 （4-14） 为了与非周期序列加以区别，周期信号序列加注下标 “p”表示周期含义。（4-14）式是仿照周期连续时间信号的傅里叶级数形式得出的周期序列傅里叶级数展开式，现在的问 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 是这样的展开是否可行，即能否找到满足（4-14）式的一组系数 。 用 乘以（4-14）式两边，并对n从0到N-1求和，即 上式右边交换求和次序，有 上式中方括弧中的和式由正交关系求出，即： 式中m为整数，方括弧中的和式只有当 或 时，取非零值N，由于后一个和式变量k的取值范围为 ，所以m必须取零值（即 ），这就是说只有当 时，方括弧中的和式取非零值，于是 （4-15） 以上分析表明，（4-14）式中的系数 可以严格地由（4-15）式求出，也就是说（4-14）式表述的关系是存在的。将（4-14）式和（4-15）式略作修改就是周期序列的傅里叶级数表达式，即 (4-16) 式中 称为旋转因子， 为傅里叶级数的系数，在这里写成序列形式，它由下式给出： （4-17） 注： 因为 所以 注意到按照我们前面推导的结果 因子应该乘以（4-17）式，而在这里将这个因子放在（4-16）式中了，这是信号处理理论中的习惯没有特殊的含义；另外也看到我们除了将傅里叶级数的系数写成序列形式外，还加注了下标“p”，这是因为 周期序列 的傅里叶级数系数 也是以N为周期的周期序列。 任意给定一个周期序列 都可以由（4-17）式求出它的傅里叶级数的系数序列 ，也就是说，时域中的一个周期序列必定与频域中的一个周期序列一一对应，在信号处理理论中通常称（4-17）式为周期序列 的离散傅里叶级数变换（简写为DFS），即 （4-18） 而（4-16）式称为离散傅里叶级数的逆变换（简写为IDFS），即 （4-19） 我们可以把 看成时域序列 的频域表示，反之 也可看成一个频域序列 的时域表示，这就是说， 与 由（4-16）式和（4-17）式构成了时域与频域的映射关系。 现在讨论离散傅里叶级数的性质。设 和 都是周期为N的 如果把周期连续时间信号的傅里叶级数的系数 看成周期序列在频域中的映射 （即 ，则我们可以得出如下关系： 1.连续、非周期时域信号 ←映射→非周期、连续频域信号 ，它由傅里叶变换构成映射关系，即 （4-20） （4-21） 2.离散、非周期时域信号 ←映射→周期、连续频域信号 ，它由序列的傅里叶变换构成映射关系，即 （4-22） （4-23） 3.连续、周期时域信号 ←映射→非周期、离散频域信号 ，它由周期函数的傅里叶级数展开式构成映射关系，即 （4-24） （4-25） 4.离散、周期时域信号 ←映射→周期、离散频域信号 ，它由离散傅里叶级数变换构成映射关系，即 （4-26） （4-27） 以上分析实际上包含了所有可能的信号形式，注意上述映射关系有这样的对称关系：如果信号在时域中是连续的，则它的频域表达式一定是非周期的，反之若信号在频域中是连续的，则它的时域表达式一定是非周期的；如果信号在时域中是离散的，则信号在频域中的表达式一定是周期的，反之如果信号在频域中是离散的，则信号在时域中的表达式是周期的。 4-4离散傅里叶变换（DFT） 如果一个信号的时域表达式是离散的，而且是有限时宽，即 （4-28） 上式表明序列 仅在区间[0,N-1]上取非零值，通常称 为有限长序列或N点序列。事实上，在 工程 路基工程安全技术交底工程项目施工成本控制工程量增项单年度零星工程技术标正投影法基本原理 我们一次观察信号总是在有限时宽范围内进行的，这就是说一次观察信号常常是有限时宽的，对于离散信号就是有限长序列。对于信号的表述无论是在时域，还是在频域一次只能表示有限长度的信号，即我们希望对一个有限时宽的信号，它的频域表示也是个有限长的，在离散情况下，一个有限长的时域序列能否表示为一个有限长的频域序列，这就是离散傅里叶变换要解决的问题。 在介绍离散傅里叶变换之前，先讨论周期序列与有限长序列的关系。一个N点序列 ，若以N为周期做周期展开就构成一个周期为N的周期序列 ，表示一个N点序列 周期性延拓的数学描述为： （4-29） 式中 称为n对N取余数，也就是n被N除可得一个整数商m和一个介于0与N之间的整数余数l，即 （4-30） 式中m为整数（可正可负），l也为整数且 。n对N取余数就等于l，即 （4-31） 例如：给定N=8时， 当n=18时，因为 ，即 ，所以 ； 当n=-18时，因为 ，即 ，所以 ； 当n=4时，因为 ，即 ，所以 ； 当n=--4时，因为 ，即 ，所以 ； 一个N点序列通过周期性延拓可以得到一个周期序列。反之，一个周期序列取其主值区间内的值可以得到一个N点序列，即 （4-32） 式中 为一个矩形序列。 以上分析说明，一个周期序列与一个N点序列有唯一的对应关系，即（4-32）式；反之，一个N点序列也与一个周期序列有唯一对应关系，即（4-29）式。 图4-1离散傅里叶变换 综上所述，可以看到时域与频域之间存在着这样的关系：一个时域N点序列 与一个时域周期序列 存在着对应关系，而这个时域周期序列 通过离散傅里叶级数与频域中的一个周期序列 存在着对应关系，这个频域周期序列又与频域中的一个N点序列 对应，反之亦然。 图4-1清楚表示地表示了上述关系，从图中可以看到时域中的一个N点序列 确实与频域中的一个N点序列 有唯一的对应关系，这个关系由离散傅里叶变换确定。 (4-32) (4-33) （4-32）式和（4-33）式定义N点序列的离散傅里叶变换，从定义式可以看出N点序列 离散傅里叶变换 也是一个N点序列，但要注意从前面引入离散傅里叶变换的过程来看，离散傅里叶变换是隐含周期性的，因为取主值使得周期性表现不出来。 一个序列的傅里叶变换定义为： （4-34） 对于一个N点序列它的非零值区间为[0,N-1]，所以傅里叶变换为 （4-35） 前面我们已经说明 是周期为2π的ω连续函数。若对连续变量ω做离散化处理（即频域抽样），即令 ，这里k为整数，则（4-35）式可以写成： 在上式中离散信号的序列符号加了下标“p”，这是因为做这样的离散化处理得到的是一个周期序列，将上式与（4-32）是比较可以看出在主值区间内[0,N-1]，上述序列傅里叶变换的抽样值与序列的离散傅里叶变换是相等的，即 （4-36） 由此可见，N点序列的离散傅里叶变换实际上就是N点序列傅里叶变换的抽样（抽样间隔为 ）序列的主值序列，这一事实也表明了序列的傅里叶变换与序列的离散傅里叶变换是两个不同的概念，切莫混为一谈。 现在我们将要讨论离散傅里叶变换的主要性质。 设    1. 线性特性 （4-37） 式中a和b均为常数。注意上式中的 和 必须是等长度的N点序列。 2. 逆变换的另一种形式 （4-38） 式中“*”表示取共轭，这个公式意义在于告诉我们求序列的离散傅里叶变换及其逆变换可以用同一个程序来实现。 3. 圆周位移特性 序列的圆周移位（有时又称为循环移位）的定义为： （4-39） 图4-2序列圆周位移 式中m为常数， 为N点序列。上式告诉我们，所谓圆周位移是将一个N点序列作周期延拓，然后移位并取主值序列，由此可见一个N点序列作圆周移位后仍为N点序列。 序列的圆周移位可以用图4-2来说明，当 时，相当于坐标右移（或图形左移）m位；当 时，相当于坐标左移（或图形右移）m位；注意序列的位移 ，相当于将原坐标原点移至 处。 4. 圆周卷积特性 设 和 均为N点序列，则 与 的圆周卷积定义为 （4-40） 上式说明，两个N点序列的圆周卷积仍是一个N点序列。通常序列的圆周卷积用符号表示，即 图4-3序列的圆周卷积 现在举例说明圆周卷积的计算方法。设 和 均为6点序列，如图4-3所示。 5．利用圆周卷积求线性卷积 由离散卷积定理可知，两个序列的圆周卷积可以通过它们的离散傅里叶乘积的反离散傅里叶变换得到，离散傅里叶变换和反离散傅里叶变换又可用其快速算法FFT来计算，所以如果线性卷积能够用圆周卷积来实现，则可利用高效的FFT算法计算线性卷积，从而可以大大提高计算效率。下面就来讨论如何用圆周卷积计算线性卷积的问题。 已知离散线性非移变系统的单位抽样响应为 ，当给定输入 时，系统的输出可以由一下卷积求出， 假定 是个N点序列， 为M点序列，则可以肯定 与 的线性卷积也一定是个有限长序列。这是因为： 而 的非零值区间为 的非零值区间为 将这两个不等式相加，有 ，这就是 可能取非零值的区间，其长度为 ，在这个区间以外，不是 等于零，就是 等于零，因而卷积 。 根据以上分析可知，N点序列 与M点序列 的线性卷积得到一个 点序列 。为了寻找圆周卷积与循环卷积的关系，令 即 补M-1个零点构成 点有限长序列 ，令 即 补N-1个零点构成 点有限长序列 。这样 和 均为 点的有限长序列，它们的圆周卷积也 是个 点序列与 一样长。我们现在逐点比较 与 的关系。 图4-4给定序列 的周期延拓 图4-5线性卷积与圆周卷积 假定 和 的图形如图4-4所示。 因为            线性卷积 圆周卷积 根据线性卷积和圆周卷积公式分别计算 和 ，如图4-5所示。 当 时，圆周卷积 与 在 区间上有非零值交点，并且相交的情况如同线性卷积 与 在 区间上非零值相交的情况。所以在这个区间内圆周卷积与线性卷积相等。 当 时，圆周卷积 与 在 上有非零值交点，并且相交的情况如同线性卷积 与 在同一区间上非零值相交的情况。所以在这个区间内圆周卷积与线性卷积相等。 当 时，圆周卷积 与 在 上有非零值交点，并且相交的情况如同线性卷积 与 在同一区间上非零值相交的情况。所以在这个区间内圆周卷积与线性卷积相等。 总上分析可见， ，即在这种情况下线性卷积与圆周卷积完全相等。一般来说，若 是个N点序列， 是个M点序列，则通过补零的方法将 和 延长为L（ ）点的序列，则 与 的L点圆周卷积就等于 与 的线性卷积。 4-5 快速傅里叶变换 如前所述，快速傅里叶变换（FFT）是DFT的快速算法，其运算次数比按DFT的定义直接计算要大为减少。下面分析FFT的算法原理。 以N=4的DFT为例，按其定义用矩阵表示，有 (4-35) 式中矩阵W的各元素均为复数，故欲求X(k)的每一个值，均需做4次复数乘法和3次复数加法。要计算X(k)的4个值，需做16次复数乘法和12次复数加法。推而广之，计算N点的DFT，需要N2次复数乘与N(N-1)次复数加(如考虑W0 =1，需作的复数乘法的次数将略  有减少）。N越大，计算量增加得越多，当N=2048时，复数乘运算达到419万次，这样大的计算量，不可能对信号实时处理。 然而仔细观察矩阵W发现，它的矩阵元素Wnk具有周期和对称的特性，因此W的许多元素都是相同的，从而为简化DFT计算提供了条件。 （1） Wnk的周期性 （4-36） 当N=4时，则W6=W2，W9=W1 （2） Wnk的对称性 因为 故有 (4-37) 若N=4，则W3= -W-1，W2= -W0 利用Wnk 的周期性与对称性，式(4-35)中的矩阵W可简化为 可见矩阵中有许多元素是相同的。如果把序列x(n)分解为若干短序列，并与W矩阵元素巧妙地结合起来计算DFT，就可能简化运算，这就是FFT的基本思路。 下面推导FFT的算法。 设N=2β(β为整数)，则 将x(n)分解成n为偶数(even)与奇数(odd)的两序列之和，上式可写成 当n为偶数时，令n=2r；n为奇数时，令n=2r+1(r为整数)，则 由于                              ，将此结果代入上式就有 (4-38) 若记 (4-39) 则式(4-38)可改写作 (4-40) 上式表明：求N点的DFT被分解为求两个N/2点的DFT，即H(k)和G(k)。以N=8的DFT[x(n)]为例，利用(4-40)则有 (4-41) 根据式(4-39)，G(k)，H(k)都为4点DFT，它们均是以4为周期的，故有G(k+4)=G(k)，H(k+4)=H(k)。再考虑W8k 的对称性W8k+4 = -W8k (0≤k≤3)，则式(4-41)可写为 (4-42) 式(4-42)可用流图来表示，流图的形式如图4.8右半部所示。根据输入输出数据的运算结构可把流图分为四个基本运算单元，称为蝶形运算单元，如图中由G(1)，H(1)，X(1)，X(5)所构成的结构就是一个蝶形运算单元[示如图4.9(a)]。由图看出，蝶形运算单元有如下规律：左方两节点为输入节点，代表输入数值；右方两节点为输出节点，表示输出数值得叠加。运算是自左向右运行的。线旁标注的加权系数W81 ， -W81 与相应的输入数值作乘法运算，即 图4.8      把8点DFT分解为两个4点DFT的流图 X(1)与X(5)取值都仅与G(1)，H(1)有关，因此把X(1)与X(5)，和H(1)与G(1)称为对偶节点，计算对偶节点的数值[如X(1)和X(5)]，仅需做一次复数乘法[W81 H(1)]、两次复数加法[G(1)+ W81 H(1)、G(1)- W81 H(1)]。蝶形运算单元也可画成图4.9(b)的形式。 图4.9      蝶形运算单元 G(k)是4点DFT，G(k)各点数值的计算可由式(4-39)将x(2r)按r的奇偶继续分组，把4点DFT进一步分解成两个2点的DFT，即 r为偶数时，记r=2l；r为奇数时，记r=2l+1(l均为正整数)，则G(k)可写成          式中，                    ，                        。因设N=8，故A(k)，B(k) 都是2点DFT，它们与G(k) 的关系为 有A（k），B(k)计算G(k)的蝶形流图见图4.10(a)。 (a)                                        (b) 图4.10    4点DFT分解为2点DFT蝶形流图 与G(k)的分解类似，也可将H(k)分解为2点DFT组合，即 式中，                        ，                        也都是2点DFT，由C(k)， D(k)求H(k)的蝶形流图见图4.10(b)。 至此，A(k)，B(k)，C(k)，D(k)都已是2点DFT，它们均是由原始数据x(n)的两点通过一个蝶形运算单元计算而成，如对于A(k)，因 当N=8时，A(0)和A(1)可写成 这已是不可再分解的最简DFT运算。 综上所述，可得出8点DFT的蝶形流图如图4.11所示，它的逐级分解的框图见图4.12。 由图4.11所示结构得出的算法即快速傅里叶变换（FFT）。 图4.11  8点DFT蝶形流图 图4.12    8点DFT逐级分解运算框图 由图4.11可见，蝶形流图的输出序列X(k)是按k由小到大顺序排列的，而输入序列x(n) 则是按所谓的码位倒序排列的。 如果把十进制数换成二进制数，然后把这些二进制的首位至末位的顺序颠倒过来再重新换算成十进制数，这个十进制数的排列即码位倒序排列。N=8的自然顺序与码位倒序的比较适于表4.2。 表4.2    自然顺序与码位倒序  （N=8）   二进制数 二进制数的码位倒序 码位倒序后的十进制数 0 000 000 0 1 001 100 4 2 010 010 2 3 011 110 6 4 100 001 1 5 101 101 5 6 110 011 3 7 111 111 7         码位倒序的排列规律是由FFT算法所决定的。如果要输入按自然顺序排列，则输出就会变成码位倒序排列；如果输出、输入均要按自然序排列，则蝶形流图的形状会发生扭曲，造成不能即位运算，计算机内存增加等新问题。 由图4.11可计算出FFT的运算工作量。每按奇偶分解一次即可得出一级蝶形流图，当N=8时，共有㏒28=3级蝶形流图。每级蝶形流图都有N/2=4个蝶形单元，故三级共有蝶形运算单元 N/2㏒2N=4×3=12（个） 由前分析已知，每一蝶形运算单元需作一次复数乘和两次复数加运算，故8点FFT总运算次数中，复数乘为 N/2㏒2N=12（次） 复数加为 N㏒2N=24（次） 直接按DFT的定义计算，8点DFT所需总运算次数为