1、现代数字信号处理及应用,杭州电子科技大学电子信息学院 刘顺兰 电话：13067787010 Email: ,序言,信号：信息的载体 信号处理：从观测信号中获得隐含的信息。 经典信号处理：非参数化信号处理（工具：FFT） 现代信号处理：参数化信号处理 现代信号处理研究的问题 估计(estimation)：参数估计、信道估计、功率谱估计、波达方向估计、特征提取、时频分析、信号检测（多用户检测） 滤波(filtering)：自适应滤波、信号处理的机器学习 辨识(identification)：系统辨识、目标识别、信号分类、反卷积与均衡,课程特点及考核,课程内容离散时间信号和系统、离散时间平稳随机过程
2、、功率谱估计、自适应滤波、参数估计等 课程特点 现代信号处理的主要理论、方法和应用“与前沿接轨” 数学知识（矩阵分析、数理统计、最优化） 考核方式 笔试,教材与参考书,教材：何子述、夏威等.现代数字信号处理及其应用.清华大学出版社.2009参考书：1 S.M. Kay, Modern Spectral Estimation, Prentice-Hall, 19882 张贤达.现代信号处理第二版.清华大学出版社.20023 西蒙赫金.自适应滤波器原理(第4版), 电子工业出版社4何子述、夏威等.现代数字信号处理及其应用习题解答.清华大学出版社.2011,第一章 离散时间信号与系统,本章主要介绍离
3、散时间信号与系统的的基本理论，包括离散时间系统定义及LTI特性，离散时间系统的时域分析和频域分析，z变换，序列的傅里叶变换，离散傅里叶变换（DFT)，快速傅里叶变换（FFT)等。,1.1 离散时间信号序列,离散时间信号的表示：x(n)或x(n)。,图 1-1 离散时间信号x(n)的图形表示,1.1.1 序列的运算,1 序列的移位,图1-2显示了x(n)序列的延时序列w(n)=x(n-2), 即m=2时的情况。,图 1-2 图1-1序列x(n)的延时,2 序列的翻褶 x(-n)如果序列为x(n), 则x(-n)是以n=0的纵轴为对称轴将序列x(n)加以翻褶。x(n)及x(-n)如图1-3(a)、
4、(b)所示。,图 1-3 序列的翻褶 (a) x(n)序列； (b) x(-n)序列,3 序列的和两序列的和是指同序号n的序列值逐项对应相加而构成的一个新序列。,4 序列的乘积 两序列相乘是指同序号n的序列值逐项对应相乘。,5 序列的标乘序列x(n)的标乘是指x(n)的每个序列值乘以常数c。,6 累加,7 差分运算前向差分 x(n)=x(n+1)-x(n)后向差分 x(n)=x(n)-x(n-1)由此得出 x(n)=x(n-1),1.1.2 几种常用序列1 单位脉冲序列(n),单位采样序列如图1-4所示。,(1-1),图 1-4 (n)序列,2 单位阶跃序列u(n),如图 1-5 所示。,(1
5、-2),图 1-5 u(n)序列,(n)和u(n)间的关系为,这就是u(n)的后向差分。 而,令n-m=k，代入此式可得,这里就用到了累加的概念。,(1-3),(1-4),(1-5),3实指数序列,式中，a为实数。当|a|1时，序列是发散的。a为负数时，序列是摆动的，如图1-6所示。,(1-6),图 1-6 指数序列 (a) |a|1; (c) a=-|a|,4 正弦型序列 x(n)=A sin(n0+) (1-7) 式中: A为幅度; 为起始相位; 0为数字域的频率，它反映了序列变化的速率。 0=0.1时, x(n)序列如图1-7所示，该序列值每20个重复一次循环。,图 1-7正弦序列(0=
6、0.1),5 复指数序列序列值为复数的序列称为复指数序列。 复指数序列的每个值具有实部和虚部两部分。复指数序列是最常用的一种复序列：,(1-8a),或,(1-8b),式中，0是复正弦的数字域频率。,对第一种表示，序列的实部、虚部分别为,如果用极坐标表示，则,因此有：,1.1.3 序列的周期性如果对所有n存在一个最小的正整数N，满足,(1-9),则称序列x(n)是周期性序列，周期为N。,现在讨论上述正弦序列的周期性。 由于,则,若N0=2k, 当k为正整数时，则,这时的正弦序列就是周期性序列，其周期满足N=2k/0（N，k必须为整数）。可分几种情况讨论如下。 （1） 当2/0为正整数时，周期为2
7、/0，见图1-8。（2） 当2/0不是整数，而是一个有理数时（有理数可表示成分数），则,式中，k, N为互素的整数，则 为最小正整数， 序列的周期为N。,（3）当2/0是无理数时，则任何k皆不能使N取正整数。 这时，正弦序列不是周期性的。 这和连续信号是不一样的。 同样，指数为纯虚数的复指数序列的周期性与正弦序列的情况相同。 下面，我们来进一步讨论，如果一个正弦型序列是由一个连续信号采样而得到的，那么，采样时间间隔T和连续正弦信号的周期之间应该是什么关系才能使所得到的采样序列仍然是周期序列呢？,设连续正弦信号xa(t)为,这一信号的频率为f0，角频率0=2f0，信号的周期为T0=1/f0=2/
8、0。 如果对连续周期信号xa(t)进行采样，其采样时间间隔为T， 采样后信号以x(n)表示，则有,如果令0为数字域频率，满足,式中, fs是采样频率。可以看出，0是一个相对频率，它是连续正弦信号的频率f0对采样频率fs的相对频率乘以2，或说是连续正弦信号的角频率0对采样频率fs的相对频率。用0代替0T， 可得,这就是我们上面讨论的正弦型序列。,下面我们来看2/0与T及T0的关系，从而讨论上面所述正弦型序列的周期性的条件意味着什么？,这表明，若要2/0为整数，就表示连续正弦信号的周期T0应为采样时间间隔T的整数倍；若要2/0为有理数，就表示T0与T是互为互素的整数，且有,(1-10),式中，k和
9、N皆为正整数，从而有,即N个采样间隔应等于k个连续正弦信号的周期。,1.1.4 用单位采样序列来表示任意序列设x(m)是一个序列值的集合，其中的任意一个值x(n)可以表示成单位采样序列的移位加权和，即,(1-11),1.1.5 序列的能量序列x(n)的能量E定义为序列各采样样本的平方和， 即,(1-12),1.2 离散时间系统时域分析,一个离散时间系统是将输入序列变换成输出序列的一种运算。若以T来表示这种运算，则一个离散时间系统可由图1-8来表示，即,(1-13),离散时间系统中最重要、 最常用的是“线性时不变系统”。,图 1-8 离散时间系统,1.2.1 线性系统满足叠加原理的系统称为线性系
10、统，即若某一输入是由N个信号的加权和组成，则输出就是系统对这几个信号中每一个的响应的同样加权和组成。 如果系统在x(n)和x2(n)单独输入时的输出分别为y1(n)和y2(n) 即:,那么当且仅当式（1-14a）和式（1-14b）成立时，该系统是线性的,（1-14b）,这两个性质合在一起就成为叠加原理，写成,（1-15）,式(1-15)对任意常数a1和a2都成立。,齐次性,可加性,（1-14a）,1.2.2 时不变系统系统的运算关系T在整个运算过程中不随时间（也即不随序列的延迟）而变化，这种系统称为时不变系统（或称移不变系统）。这个性质可用以下关系表达：若输入x(n)的输出为y(n)， 则将输
11、入序列移动任意位后， 其输出序列除了跟着移位外， 数值应该保持不变，即若 Tx(n)=y(n) 则 Tx(n-m)=y(n-m) （m为任意整数） (1-16) 满足以上关系的系统就称为时不变系统。同时满足线性和时不变性的系统，称为线性时不变系统（linear time-invariant system），常简称为LTI系统。,1.2.3 单位脉冲响应与系统的输入输出关系线性时不变系统可用它的单位脉冲响应来表征。 单位脉冲响应是指输入为单位脉冲序列时系统的输出。一般用h(n)表示单位脉冲响应，即 h(n)=T(n) 有了h(n)我们就可以得到此线性时不变系统对任意输入的输出。 下面讨论这个问题
12、： 设系统输入序列为x(n)，输出序列为y(n)。从式（1-11）已经知道，任一序列x(n)可以写成(n)的移位加权和， 即,则线性时不变系统的输出为,因此,（1-17a）,或,（1-17b）,如图1-9所示。上式称为序列x(n)与h(n)的离散卷积，也称为“线性卷积”或“直接卷积”或简称“卷积”，并以“*”表示之。,图 1-9 线性时不变系统,1.2.5 因果系统所谓因果系统，就是系统此时的输出y(n)只取决于此时， 以及此时以前的输入，即x(n), x(n-1), x(n-2), 。如果系统的输出y(n)还取决于x(n+1), x(n+2), ，也即系统的输出还取决于未来的输入，这样在时间
13、上就违背了因果关系，因而是非因果系统， 也即不现实的系统。根据上述定义，可以知道，y(n)=nx(n)的系统是一个因果系统，而y(n)=x(n+2)+ax(n)的系统是非因果系统。,从式（1-17）卷积公式，我们可以看到线性时不变系统是因果系统的充分必要条件是h(n)=0 n0 (1-18) 依照此定义，我们将n0，x(n)=0 的序列称为因果序列，表示这个因果序列可以作为一个因果系统的单位脉冲响应。,1.2.6 稳定系统稳定系统是指有界输入产生有界输出（BIBO）的系统。如果对于输入序列x(n)，存在一个不变的正有限值Bx，对于所有n值满足 |x(n)|Bx （1-19） 则称该输入序列是有
14、界的。稳定性要求对于每个有界输入存在一个不变的正有限值By，对于所有n值，输出序列y(n)满足 |y(n)|By （1-20）,一个线性时不变系统是稳定系统的充分必要条件是单位脉冲响应绝对可和， 即,(1-21),因果稳定系统：,这种稳定因果系统既是可实现的，又是稳定工作的，因而这种系统正是一切数字系统设计的目标。,1.2.7 常系数线性差分方程连续时间线性时不变系统的输入输出关系常用常系数线性微分方程表示，而离散时间线性时不变系统的输入输出关系除了用式(1-17)表示外，常用以下形式的常系数线性差分方程表示， 即,(1-23),所谓常系数是指决定系统特征的a1,a2,aN, b1, b2,
15、, bM都是常数。差分方程的阶数等于未知序列（指y(n)）变量序号的最高值与最低值之差。 例如, 式（1-23）即为N阶差分方程。注： 当系统只有输入信号，初始状态为零时，系统的响应为零状态响应。反之，当系统输入信号为零，初始状态不为零时，系统的响应称为零输入响应。,MATLAB实现,已知输入和差分方程系数，可利用filter函数对差分方程进行数值求解。如果系统的初始条件为零，则函数的调用格式为：y=filter(b, a, x)例1 求解系统单位脉冲响应的程序example15.m如下：%example15.m:调用filter函数解差分方程y(n)=x(n)+0.5y(n-1) b=1;
16、a=1, -0.5;x=1, zeros(1,31); %x(n)为单位脉冲序列,长度为32y=filter(b,a,x); %计算出单位脉冲响应h(n)n=0:31;stem(n,y);xlabel(n);ylabel(y(n);,图1-10 系统的单位脉冲响应,执行结果如下：,1.3.1 Z变换的定义及收敛域 1. Z变换的定义一个离散序列x(n)的Z变换（ 双边Z变换）定义为,1.3 Z 变 换,式中，z是一个复变量，它所在的复平面称为Z平面。我们常用Zx(n)表示对序列x(n)进行Z变换，也即,(1-24),(1-25),2. Z变换的收敛域对任意给定序列x(n)，使其Z变换收敛的所有
17、z值的集合称为X(z)的收敛域。一般收敛域用环状域表示，即Rx-|z|Rx+ 收敛域是分别以Rx-和Rx+为半径的两个圆所围成的环状域（图1-11 中的斜线部分）。Rx-和Rx+称为收敛半径。当然Rx-可以小到零，R x+可以大到无穷大。 ,图1-11 环形收敛域,分子多项式P(z)的根是X(z)的零点，分母多项式Q(z)的根是X(z)的极点。在极点处Z变换不存在，因此收敛域中没有极点， 收敛域总是用极点限定其边界。 ,常用的Z变换是一个有理函数，用两个多项式之比表示：,收敛域的位置Rx-及Rx+的大小和序列的关系,（1）有限长序列:,其Z变换为,有限长序列的收敛域表示如下：,（2）右边序列：
18、右边序列是指x(n)只在nn1时有值，在nn1时x(n)=0。,右边序列Z变换的收敛域为 Rx-|z| 右边序列及其收敛域如图1-12所示。,图1-12 收敛域 （n10， |z|=除外）,因果序列即n1=0的右边序列。 也就是说，在n0时x(n)有值，n0时x(n)=0.，其Z变换的收敛域包括|z|=。,Z变换收敛域包括|z|=是因果序列的特征。,(1-27),收敛域上函数必须是解析的，因此收敛域内不允许有极点存在。 所以，右边序列的Z变换如果有N个有限极点z1,z2,zN存在， 那么收敛域一定在模值为最大的这一个极点所在圆以外，也即,对于因果序列，处也不能有极点。,（3） 左边序列: 左边
19、序列是指在nn2时x(n)有值，而在nn2时x(n)=0.,如果n20，收敛域应包括z=0， 即|z|Rx+。,左边序列Z变换的收敛域为,对于左边序列，如果序列Z变换有N个有限极点z1, z2, , zN存在，那么收敛域一定在模值为最小的这一个极点所在圆以内， 这样X(z)才能在整个圆内解析， 也即 Rx+=min|z1|, |z2|, , |zN| ,(1-28),（4） 双边序列：一个双边序列可以看作一个右边序列和一个左边序列之和。,如果Rx-Rx+，则无公共收敛区域，X(z)无收敛域，也即在Z平面的任何地方都没有有界的(z)值，因此就不存在Z变换的解析式， 这种变换就没有什么意义。,(1
20、-29),表1-1 几种序列的Z变换,1.3.2 Z反变换已知函数X(z)及其收敛域，反过来求序列的变换称为Z反变换，表示为 x(n)=Z-1X(z) Z反变换的一般公式为,(1-30),图1-13 围线积分路径,一般求Z反变换的常用方法有三种： 围线积分法（留数法）、部分分式展开法和幂级数展开法。,表 1-2 Z变换的主要性质,正变换,(1-31),反变换,(1-32),其收敛条件为,(1-33),绝对可加性是傅里叶变换表示存在的一个充分条件。也就是说， 若序列x(n)绝对可和，则它的傅里叶变换一定存在且连续。,1.4 序列的傅氏变换,1.4.1 序列的傅氏变换的定义,由序列z变换的定义和傅
21、氏变换的定义可得：,(1-34),可见，单位圆上的Z变换即为序列的傅里叶变换，也称为数字序列的频谱。,1.4.2 序列的傅氏变换与z变换的关系,表1-3 序列傅里叶变换的主要性质,表1-3 序列傅里叶变换的主要性质,表1-4 傅里叶变换对,1.5 离散时间系统的频域分析（域和Z域）,在1.2节中已经讨论过，在时域中，一个线性时不变系统完全可以由它的单位脉冲响应h(n)来表示。对于一个给定的输入x(n)，其输出y(n)为,对等式两端取Z变换，得,则,(1-35),我们把H(z)定义为线性时不变系统的系统函数，它是单位脉冲响应的Z变换，即,(1-36),在单位圆上(z=ej)的系统函数就是系统的频
22、率响应H(ej)。,(1-37),1.5.1 因果系统单位脉冲响应h(n)为因果序列的系统称为因果系统， 因此由1.4.1节可知因果系统的系统函数H(z)具有包括z=点的收敛域，即,(1-38),1.5.2 稳定系统由1.2节中的讨论已知，一个线性时不变系统稳定的充分必要条件为h(n)必须满足绝对可和条件，即,而Z变换的收敛域由满足 的那些z值确定，因此稳定系统的系统函数H(z)必须在单位圆上收敛，即收敛域包括单位圆|z|=1，H(ej)存在。,1.5.3 因果稳定系统因果稳定系统是最普遍、最重要的一种系统，它的系统函数H(z)必须在从单位圆到的整个Z域内收敛，即,(1-39),也就是说，系统
23、函数的全部极点必须在单位圆内。,1.5.4 系统函数和差分方程的关系1.2节中已说明，一个线性时不变系统也可以用常系数线性差分方程来表示， 其N阶常系数线性差分方程的一般形式为,若系统起始状态为零，这样就可以直接对上式两端取Z变换， 利用Z变换的线性特性和移位特性可得,系统的频率响应为,（1-40）,这样就得到系统函数为,（1-41）,式中，z=ck是H(z)的零点，z=dk是H(z)的极点，它们都由差分方程的系数ak和bk决定。因此，除了比例常数b0/a0以外，系统函数完全由它的全部零点、极点来确定。 ,（1-42）,由此看出系统函数分子、分母多项式的系数分别就是差分方程的系数。式（1-40
24、）是两个z-1的多项式之比，将其分别进行因式分解，可得,1.5.5 有理系统函数的单位脉冲响应（IIR，FIR）,在线性时不变系统中，分成两类不同的系统： 若系统的单位脉冲响应延伸到无穷长，称之为“无限长单位脉冲响应系统”，简写为IIR系统。 若系统的单位脉冲响应是一个有限长序列， 称之为“有限长单位脉冲响应系统”，简称为FIR系统。,IIR系统函数H(z)至少有一个非零极点。因此如果在有限Z平面出现极点，那么这个系统就是IIR系统。 FIR系统函数H(z)在有限Z平面0|z|收敛。也就是说，H(z)在有限Z平面不能有极点，只存在零点。,1.5.6 系统频率响应的意义为了研究离散线性系统对输入
25、频谱的处理作用， 有必要研究线性系统对复指数或正弦序列的稳态响应，即系统的频域表示法。 对于稳定系统，如果输入序列是一个频率为的复正弦序列： x(n)=ejn -n线性时不变系统的单位脉冲响应为h(n)，则其输出为,式中:,因此y(n)=ejnH(ej) (1-43)上式表明，当线性时不变系统输入是频率为的复正弦序列时，输出为同频复正弦序列乘以加权函数H(ej)。显然，H(e j)描述了复正弦序列通过线性时不变系统后，幅度和相位随频率的变化。换句话说，系统对复正弦序列的响应完全由H(e j)决定。故称H(ej)为线性时不变系统的频率响应。线性时不变系统的频率响应是其单位脉冲响应的傅里叶变换。,
26、线性时不变系统的频率响应H(ej)是以2为周期的连续周期函数, 是复函数。它可以写成模和相位的形式,式中，频率响应的模|H(ej)|叫做振幅响应（或幅度响应）， 频率响应的相位arg|H(ej)|叫做系统的相位响应。 系统频率响应H(ej)存在且连续的条件是h(n)绝对可和， 即要求系统是稳定系统。,线性时不变系统在任意输入情况下，输入与输出两者的傅里叶变换间的关系，可通过对卷积公式两端取傅里叶变换，并利用表1-3傅里叶变换性质10得到 Fy(n)=Fx(n) * h(n) 即Y(ej)=X(ej)H(ej) (1-44) H(ej)就是系统的频率响应。 由式（1-44）得知，对于线性时不变系
27、统，其输出序列的傅里叶变换等于输入序列的傅里叶变换与系统频率响应的乘积。,若对Y(ej)取傅里叶反变换，可求得输出序列为,(1-45),若用极坐标形式表示频率响应，则系统的输入和输出的傅里叶变换的振幅和相位间的关系可表示为:,|Y(ej)|=|H(ej)|X(ej)| argY(ej)=argH(ej)+argX(ej),(1-46),(1-47),数字滤波器的理想幅频特性,1.5.7 频率响应的几何确定法观察式（1-42）可以发现，一个N阶的系统函数H(z)完全可以用它在Z平面上的零、极点确定。由于H(z)在单位圆上的Z变换即是系统的频率响应，因此系统的频率响应也完全可以由H(z)的零、极点
28、确定。频率响应的几何确定法实际上就是利用H(z)在Z平面上的零、极点，采用几何方法直观、定性地求出系统的频率响应。式（1-42）已表示出H(z)的因式分解， 即用零、极点表示为,(1-48),假设M=N，这时用z=ej代入，即得系统的频率响应为,(1-49),在Z平面上，ej-ck可以用一根由零点ck指向单位圆上ej点的向量Ck来表示 Ck=ej-ck 同样，ej-dk可以由极点dk指向单位圆上ej的向量Dk来表示 Dk=ej-dk,因此,(1-50),以极坐标表示有： Ck=Ckejk Dk=Dke jk H(ej)=|H(ej)|ej(),就得到,(1-51),(1-52),这样频率响应的
29、幅度函数就等于各零点至ej点向量长度之积除以各极点至ej点向量长度之积，再乘以常数b0/a0。而频率响应的相位函数等于各零点至ej点向量的相角之和减去各极点至ej点向量相角之和。当频率由0到2时，这些向量的终端点沿单位圆反时针方向旋转一圈，从而可以估算出整个系统的频率响应来。例如，图1-14表示了具有两个极点一个零点的系统以及它的频率响应，这个频率响应不难用几何法加以验证。,图 1-14 频率响应的几何表示法,极点对频谱的峰值影响：从（1-51）式看到，当在某个极点附近时，这时向量最短，出现极小值，因而频率响应在这附近可能出现峰值，同时极点越靠近单位圆，的极小值越小，频率响应出现的峰值就将越尖
30、锐。当极点处在单位圆上时，的极小值为零，在所在点的频率响应将出现。极点在单位圆外则不稳定。 零点对频谱的谷点影响：零点的影响则相反，从（1-51）式看到，越接近某零点，频率响应越低，因此在零点附近，频率响应将出现谷点，零点越接近单位圆，谷点越接近零。当零点处在单位圆上时，谷点为零，也即在零点所在频率上出现传输零点。零点可在单位圆外，不受稳定性约束。,零、极点位置对系统的频率响应的影响,1.6 周期序列的离散傅里叶级数(DFS),设 是一个周期为N的周期序列， 即,r为任意整数,周期序列不是绝对可和的，所以不能用Z变换表示，因为在任何z值下，其Z变换都不收敛，也就是,（1-53）,（1-54）,
31、式中，DFS表示离散傅里叶级数正变换，IDFS表示离散傅里叶级数反变换。WN定义为,周期序列可以用离散傅里叶级数来表示，其定义如下：,（1-55）,可以看出，当0kN-1 时， 是对X(z)在Z平面单位圆上的N点等间隔采样，在此区间之外随着k的变化， 的值呈周期变化。 图1-15画出了这些特点。,式（1-53）中的周期序列 可看成是对 的第一个周期x(n)作Z变换，然后将Z变换在Z平面单位圆上采样而得到的。即,与X(z)的关系：,由于单位圆上的Z变换即为序列的傅里叶变换，所以周期序列 也可以解释为 的一个周期x(n)的傅里叶变换的等间隔采样。 即,这相当于以2/N的频率间隔对傅里叶变换进行采样
32、。,（1-56）,与 的关系：,1.7 有限长序列离散傅里叶变换（DFT）,设x(n)为有限长序列，长度为N，即x(n)只在n=0到N-1点上有值，其他n时，x(n)=0。即,1.7.1 DFT的定义,有限长序列的离散傅里叶变换的定义为:,0kN-1,0nN-1,（1-57）,（1-58）,式(1- 57)称为x(n)的N点离散傅里叶变换(DFT) 式(1-58)称为X(k)的N点离散傅里叶反变换(IDFT),1.7.2 DFT与序列傅里叶变换、Z变换的关系若x(n)是一个有限长序列，长度为N，对x(n)进行Z变换,比较Z变换与DFT，我们看到，当z=W-kN时,即,（1-59）,表明 是Z平
33、面单位圆上幅角为 的点，也即将Z平面单位圆N等分后的第k点，所以X(k)也就是对X(z)在Z平面单位圆上N点等间隔采样值，如图1-16所示。此外， 由于序列的傅里叶变换X(ej)即是单位圆上的Z变换，根据式（1-59）， DFT与序列傅里叶变换的关系为,（1-60）,（1-61）,式（1-60）说明X(k)也可以看作序列x(n)的傅里叶变换X(ej)在区间0, 2上的N点等间隔采样，其采样间隔为N=2/N， 这就是DFT的物理意义。显而易见，DFT的变换区间长度N不同， 表示对X(ej)在区间0, 2上的采样间隔和采样点数不同， 所以DFT的变换结果也不同。,图 1-16 DFT与序列傅里叶变
34、换、Z变换的关系,信号时域采样理论实现了信号时域的离散化，使我们能用数字技术在时域对信号进行处理。而离散傅里叶变换理论实现了频域离散化，因而开辟了用数字技术在频域处理信号的新途径，从而推进了信号的频谱分析技术向更深更广的领域发展。,表 1-5 DFT性质表(序列长皆为点),1.8 快速傅里叶变换(FFT),快速傅里叶变换（FFT）并不是一种新的变换， 而是离散傅里叶变换（DFT）的一种快速算法。 离散傅里叶变换（DFT）在数字信号处理中是非常有用的。例如，在信号的频谱分析、 系统的分析、 设计和实现中都会用到DFT的计算。 但是，由于DFT的计算量太大，很难进行实时处理，所以没有得到真正的运用
35、。 FFT出现后使DFT的运算大大简化，运算时间一般可缩短一二个数量级之多，从而使DFT的运算在实际中真正得到了广泛的应用。,设x(n)为N点有限长序列，其DFT为,k=0, 1, , N-1,（1-62）,反变换（IDFT）为,n=0, 1, , N-1,（1-63）,一般来说，x(n)和WN nk都是复数，X(k)也是复数，因此每计算一个X(k)值，需要N次复数乘法和N-1次复数加法。,而X(k)一共有N个点（k从0取到N-1），所以完成整个DFT运算总共需要N2次复数乘法及N(N-1)次复数加法。一次复数乘法需用四次实数乘法和二次实数加法；一次复数加法需二次实数加法。 因而每运算一个X(
36、k)需4N次实数乘法和2N+2(N-1)=2(2N-1)次实数加法。 所以，整个DFT运算总共需要4N2次实数乘法和2N(2N-1)次实数加法。,可以看到，直接计算DFT，乘法次数和加法次数都是和N2成正比的，当N很大时，运算量是很可观的，有时是无法忍受的。,例3-1 根据式（3-1），对一幅NN点的二维图像进行DFT变换，如用每秒可做10万次复数乘法的计算机，当N=1024时，问需要多少时间（不考虑加法运算时间）？解 直接计算DFT所需复乘次数为（N2)21012次，因此用每秒可做10万次复数乘法的计算机，则需要近3000小时。 这对实时性很强的信号处理来说，要么提高计算速度，而这样，对计算
37、速度的要求太高了。另外，只能通过改进对DFT的计算方法，以大大减少运算次数。,改善途径能否减少运算量，从而缩短计算时间呢？仔细观察DFT的运算就可看出， 利用系数Wnk的以下固有特性，就可减少运算量： （1） WNnk的对称性,（2） WNnk的周期性,（3） WNnk的可约性,另外,这样，利用这些特性，使DFT运算中有些项可以合并，并能使DFT分解为更少点数的DFT运算。而前面已经说到，DFT的运算量是与N2成正比的，所以N越小越有利，因而小点数序列的DFT比大点数序列的DFT的运算量要小。 快速傅里叶变换算法正是基于这样的基本思想而发展起来的。 它的算法形式有很多种，但基本上可以分成两大类
38、，即按时间抽取(ecimationinTime，缩写为DIT）法和按频率抽取（Decimationinrequency， 缩写为DIF）法。,FFT算法与直接计算DFT运算量的比较对于长度为N=2M 的DFT计算，采用FFT方法所需的计算量为,由于计算机上乘法运算所需的时间比加法运算所需的时间多得多，故以乘法为例，直接DFT复数乘法次数是N2，FFT复数乘法次数是(N/2) log2N。 直接计算DFT与FFT算法的计算量之比为,当N=2048时，这一比值为372.4，即直接计算DFT的运算量是FFT运算量的372.4倍。当点数N越大时，FFT的优点更为明显。,例3-2 用FFT算法处理一幅N
39、N点的二维图像，如用每秒可做10万次复数乘法的计算机，当N=1024时，问需要多少时间（不考虑加法运算时间）？解 当N=1024点时，FFT算法处理一幅二维图像所需复数乘法约为 次，仅为直接计算DFT所需时间的10万分之一。 即原需要3000小时，现在只需要2 分钟。,1.9 连续时间信号的离散化及其频谱关系,在某些合理条件限制下，一个连续时间信号能用其采样序列来完全给予表示，连续时间信号的处理往往是通过对其采样得到的离散时间序列的处理来完成的。本节将讨论信号的采样过程以及采样后信号的频谱将发生怎样的变换。,采样器可以看成是一个电子开关，它的工作原理可由图1-17（a）来说明。设开关每隔T秒短
40、暂地闭合一次，将连续信号接通， 实现一次采样。如果开关每次闭合的时间为秒，那么采样器的输出将是一串周期为T，宽度为的脉冲。而脉冲的幅度却是重复着在这段时间内信号的幅度。如果以xa(t)代表输入的连续信号，如图1-17（b）所示，以xp(t)表示采样输出信号，它的结构如图1-17（d）所示。显然，这个过程可以把它看作是一个脉冲调幅过程。被调制的脉冲载波是一串周期为T、宽度为的矩形脉冲信号，如图1-9（c）所示，并以p(t)表示，而调制信号就是输入的连续信号。因而有,一般开关闭合时间都是很短的，而且越小，采样输出脉冲的幅度就越准确地反映输入信号在离散时间点上的瞬时值。当T时，采样脉冲就接近于函数性
41、质。,图 1-17 连续时间信号的采样过程,1.9.1 理想采样理想采样就是假设采样开关闭合时间无限短，即0的极限情况。此时，采样脉冲序列p(t)变成冲激函数序列s(t)，如图 1-17(e)所示。这些冲激函数准确地出现在采样瞬间，面积为1。采样后，输出理想采样信号的面积（即积分幅度）则准确地等于输入信号xa(t)在采样瞬间的幅度。理想采样过程如图1-17（f）所示。 冲激函数序列s(t)为,(1-64),以 表示理想采样的输出，以后我们都以下标a表示连续信号（或称模拟信号），如xa(t)；而以它的顶部符号（）表示它的理想采样，如 。这样我们就可将理想采样表示为,(1-65),把式（1-64）
42、代入式（1-65），得,(1-66),(1-67),1.9.2 理想采样信号的频谱我们首先看看通过理想采样后信号频谱发生了什么变化。 由于在连续时间信号与系统中已学过，式（1-65）表示时域相乘， 则频域（傅里叶变换域）为卷积运算。 所以由式（1-65）可知， 若各个信号的傅里叶变换分别表示为:,(1-68),(1-69),(1-70),则应满足,(1-71),或者,(1-72),(1-73),(1-74),其中,一般称fs为频率，单位为赫兹(Hz)，s为角频率，单位为弧度/秒; 习惯上都统称为“频率”。 它们的区别由符号f及来识别。,由此看出，一个连续时间信号经过理想采样后，其频谱将沿着频率
43、轴以采样频率s=2/T 为间隔而重复，这就是说频谱产生了周期性延拓，如图1-18所示。也就是说，理想采样信号的频谱， 是Xa(j)的周期延拓函数，其周期为s，而频谱的幅度则受1/T加权，由于T是常数，所以除了一个常数因子外，每一个延拓的谱分量都和原频谱分量相同。因此只要各延拓分量与原频谱分量不发生频率混叠，则有可能恢复出原信号。也就是说，如果xa(t)是限带信号，其频谱如图1-18（a）所示，且最高频谱分量h不超过s/2，即,(1-75),那么原信号的频谱和各次延拓分量的谱彼此不重叠，如图1-18（c）所示。 这时采用一个截止频率为s/2的理想低通滤波器， 就可得到不失真的原信号频谱。也就是说
44、，可以不失真地还原出原来的连续信号。,图 1-18 时域采样后， 频谱的周期延拓 （a）原始限带信号频谱; （b） 采样函数频谱; （c）已采样信号频谱（s2h）;（d）已采样信号频谱（s2h）,如果信号的最高频谱h超过s/2，则各周期延拓分量产生频谱的交叠，称为混叠现象，如图1-18（d）所示。由于Xa(j)一般是复数，所以混叠也是复数相加。为了简明起见，在图1-18中我们将Xa(j)作为标量来处理。 我们将采样频率之半（s/2）称为折叠频率，即,它如同一面镜子，当信号频谱超过它时，就会被折叠回来，造成频谱的混叠。,(1-76),由此得出结论：要想采样后能够不失真地还原出原信号， 则采样频率
45、必须大于两倍信号谱的最高频率（s2h），这就是奈奎斯特采样定理。 即 fs2fh 频率h一般称为奈奎斯特频率，而频率2h称为奈奎斯特率。 采样频率必须大于奈奎斯特率。 在实际工作中，为了避免频谱混淆现象发生，采样频率总是选得比奈奎斯特频率更大些，例如选到(34)h。同时为了避免高于折叠频率的杂散频谱进入采样器造成频谱混淆，一般在采样器前加入一个保护性的前置低通滤波器，称为防混叠滤波器，其截止频率为s/2，以便滤除掉高于s/2 的频率分量。,(1-77),离散时间信号,1.9.3 离散时间信号（序列）的频谱,(1-78),其傅里叶变换为：,可看成是连续时间信号的采样，即,理想采样信号为：,对其两
46、边进行傅里叶变换，可得：,(1-79),比较式（1-78）、式（1-79）可以发现，离散时间信号傅里叶变换,(1-80),之间有如下关系：,与连续时间信号傅里叶变换,或,(1-81),1.10 利用FFT分析时域连续信号频谱,1.10.1 基本步骤,图 1-19 时域连续信号离散傅里叶分析的处理步骤,xc(t)通过A/D变换器转换(忽略其幅度量化误差)成采样序列x(n)，其频谱用X(ej)表示，它是频率的周期函数，即,(1-82),式中，Xc(j)或 为xc(t)的频谱。,由于进行FFT的需要，必须对序列x(n)进行加窗处理，即v(n)=x(n)w(n)，加窗对频域的影响，用卷积表示如下:,最
47、后是进行FFT运算。 加窗后的DFT是,0kN-1,(1-83),(1-84),有限长序列v(n)=x(n)w(n)的DFT相当于v(n)傅里叶变换的等间隔采样。,(1-85),V(k)便是sc(t)的离散频率函数。因为DFT对应的数字域频率间隔为=2/N，且模拟频率和数字频率间的关系为=， 其中=2f。所以，离散的频率函数第k点对应的模拟频率为:,(1-86),(1-87),由式（1-87）很明显可看出，数字域频率间隔=2/N对应的模拟域谱线间距为,(1-88),谱线间距，又称频谱分辨率（单位：Hz）。所谓频谱分辨率是指可分辨两频率的最小间距。它的意思是，如设某频谱分析的F=5Hz，那么信号
48、中频率相差小于5 Hz的两个频率分量在此频谱图中就分辨不出来。,长度N=16 的时间信号v(n)=(1.1)nR16(n)的图形如图 1-20(a)所示, 其16点的DFT V(k)的示例如图 1 - 20(b)所示。 其中T为采样时间间隔（单位：s）；fs为采样频率（单位：Hz）；tp为截取连续时间信号的样本长度（又称记录长度，单位：s）；F为谱线间距，又称频谱分辨率（单位：Hz）。 注意：V(k)示例图给出的频率间距F及N个频率点之间的频率fs为对应的模拟域频率(单位： Hz)。,图 1 - 20,由图可知：,(1-89),(1-90),在实际应用中, 要根据信号最高频率fh和频谱分辨率F的要求， 来确定T、tp和N的大小。 （1）首先，由采样定理，为保证采样信号不失真，fs2fh(fh为信号频率的最高频率分量，也就是前置低通滤波器阻带的截止频率)， 即应使采样周期T满足,
