第一部分：
绪 论
MATLAB是一个功能十分强大的系统，是集数值计算、图形管理、程序开发为一体的环境。作为强大的科学计算平台，它几乎能满足所有的计算需求。在美国及其他发达国家的理工科院校里，MATLAB已成为了一门必修的课程，在科研院所、大型公司或企业的工程计算部门，MATLAB也是最为普遍的计算工具之一。
MATLAB具有如下的优势和特点：
1、好的工作平台和编程环境随着MATLAB的商业化以及软件本身的不断升级，MATLAB的用户界面也越来越精致，更加接近Windows的标准界面，人机交互性更强，操作更简单。而且新版本的MATLAB提供了完整的联机查询、帮助系统，极大地方便了用户的使用。简单的编程环境提供了比较完备的调试系统，程序不必经过编译就可以直接运行，而且能够及时地报告出现的错误及进行出错原因分析。
2、简单易用的程序语言新版本的MATLAB语言是基于最为流行的C语言基础上的，因此语法特征与C语言极为相似，而且更加简单，更加符合科技人员对数学表达式的书写格式。使之更利于非计算机专业的科技人员使用。而且这种可移植性好、可扩展性极强，这也是MATLAB之所以能够深入到科学研究及工程计算各个领域的重要原因。
3、强大的科学计算及数据处理能力
MATLAB拥有600多个工程中要用的数学运算函数，可以方便地实现用户所需的各种计算功能函数中所使用的算法都是科研和工程计算中的最新研究成果，而且经过了各种优化及容错处理，因此使用起来鲁棒性和可靠性非常高。在通常情况下，可以用它来代替底层编程语言，如C和C++。在计算要求相同的情况下，使用MATLAB的编程工作量会大大减少。MATLAB函数所能解决的问题包括矩阵运算和线性方程组的求解、微分及偏微分方程组的求解、符号运算、傅立叶变换和数据的统计分析、工程中的优化问题、稀疏矩阵运算、复数的各种运算、三角函数和其他初等数学运算、多维数组操作以及建模动态仿真等。
4、出色的图形处理功能
MATLAB自产生之日起就具有方便的数据可视化功能，新版本的MATLAB对整个图形处理的功能做了很大的改进和完善，使它不仅在一般数据可视化软件都具有的功能（例如二维曲线和三维曲面的绘制和处理等）方面更加完善，而且对于一些其他软件所没有的功能（例如图形的光照处理、色度处理以及四维数据的表现等），MATLAB同样表现了出色的处理能力。同时对一些特殊的可视化要求，例如图形动画等，MATLAB也有相应的功能函数，保证了用户不同层次的要求。另外，新版本的MATLAB还着重在图形用户界面（GUI）的制作上做了很大的改善，对这方面有特殊要求的用户也可以得到满足。
5、应用广泛的模块集或工具箱
MATLAB对许多专门的领域都开发了功能强大的模块集或工具箱。一般来说，它们都是由特定的专家开发的，用户可以直接使用工具箱学习、应用和评估不同的方法而不需要自由编写代码。目前，MATLAB已经把工具延伸到科学研究和工程应用的诸多领域，诸如数据采集、数据库接口、概率统计、样条拟合优化算法、偏微分方程求解、神经网络、小波分析、信号处理、图象处理、系统辩识、控制系统设计、LMI控制、鲁棒控制、模型预测、模糊逻辑、金融分析、地图工具、非线性控制设计、实时快速原型及半物理仿真、嵌入式系统开发、定点仿真、DSP与通讯、电力系统仿真等，都在工具箱（Toolbox）家族中有了自己的一席之地。
6、使用的程序接口和开发平台新版本的MATLAB可以利用MATLAB编译器和C/C++数学库和图形库，将自己的MATLAB程序自动转换为独立于MATLAB运行的C和C++代码。另外，MATLAB网页服务程序还容许在Web应用中使用自己的MATLAB数学和图形程序。
7、模块化的设计和系统级的仿真
Simulink是MATLAB的一个分支产品，主要用来实现对工程问题的模型化和动态仿真。在世界范围内的模型化浪潮的背景下Simulink恰恰体现了模块化设计和系统级仿真的具体思想，使得建模仿真如同搭积木一样简单。Simulink对仿真的实现可以应用于动力系统、信号控制、通信设计、金融财会及生物医学等各个领域的研究中。
一、基本概念
1、变量和其他高级语言一样，MATLAB也是使用变量来保存信息。变量由变量名表示，变量的命名应遵循如下规则：
1）变量名必须以字母开头。
2）变量名可以由字母、数字和下划线混合组成。
3）变量名区分字母大小写。
4）变量名的字符长度不应超过31个。
在MATLAB中还存在着一些系统默认的固定变量，如表2-1所示，即在MATLAB中语句中若出现固定变量名，则系统就将赋予默认值。
变量名
默认值
i
虚数单位
j
虚数单位
pi
圆周率
inf
无穷大
MATLAB的变量分为字符变量和数值变量两种，字符变量必须用但引号括起来。例如，用户可输入：
a=‘happly new year’
则表示将字符串 ‘happly new year’ 赋值给字符变量a.。
若用户输入：b=365
则表示将数值365赋值给数值变量b。
和其他高级语不同的是，MATLAB使用变量时不需要预先对变量类型进行说明，MATLAB会自动根据所输入的数据来决定变量的饿数据类型和分配存储空间。
2、数值在MATLAB内部，每个数据元素都是用双精度来表示和存储的，大约有16位有效数字。其数值有效范围约为。
但在进行数据输入输出时，MATLAB却可以用不同的格式。如果参加运算的每个元素均为整数，则MATLAB将用不加小数点的纯整数格式显示运算结果，否则，按默认的输出格式显示结果。MATLAB的默认格式为short格式，该格式显示运算结果为保留小数点后4位有效数字。用户可以通过format命令改变输出格式为long，以得到更多的有效数字（小数点后14位）。需要注意的是，数据输出格式的改变，并不影响该数据在MATLAB内部的存储精度。设置为short很long输出格式的命令分别为：
format short
format long
MATLAB通常用十进制数来表示一个数，亦可用科学计数法来表示一个数。另外，MATLAB还可以进行复数运算，复数可以由如下语句来产生：
C=a+i*b (或c=a+j*b) 将实部a虚部为b的复数赋值给复变量c
C=a*exp(j*b) (或C=a*exp(j*b)) 将模为a辐角为b的复数赋值给复数c
其中I、j是虚数单位。
以下是MATLAB各种合法的输入数据示例：
365 -18 0.000076 8.2467-10
2.67i 1.02E3i -7.6 -3.4502e+20
3、矩阵
矩阵是MALTAB进行数据处理和运算的基本元素。MATLAB的大部分运算或命令都是在矩阵运算的意义下执行的。我们通常意义上的数量（标量）在MATLAB系统中是作为
的矩阵来处理的，而仅有一行或一列的矩阵在MATLAB中称为向量。
4、数组
在MATLAB中，数组也是一个重要的概念，矩阵在某些情况下可视为二阶的数值型数组。但是在MATLAB中，数组和矩阵运算规则却有较大的区别。例如，两矩阵相乘和两数组相乘遵循的运算规则就是完全不相同的。
5、函数
MATLAB为用户提供了丰富且功能各异的函数，用户可以直接用这些函数来进行数据处理。函数由函数名和参数组成，函数调用的格式为：
函数名（参数）
例如，若在MATLAB的命令窗口输入命令：
a=sin(b)
则表示计算b的正弦值并将赋值给变量a.。
二、多项式处理
1、多项式表示多项式在MATLAB中使用降幂系数的行向量表示。例如：多项式x4-12x3+0x2+25x+116表示为：p=[1 -12 0 25 116]
需要注意的是，表示中需要包含零系数的项。按照这种形式，使用函数roots可以找出多项式等于零的根：r=roots(p)
r =11.74728287082272,2.70282074384101,-1.22505180733187 + 1.46720800872890i、-1.22505180733187,- 1.46720800872890I
按照MATLAB的规定，多项式使用行向量，根使用列向量。给出多项式的根，使用poly函数也可以构造出相应的多项式：
u=poly(r)
pp =1.0000 -12.0000 0 25.0000 116.0000
2、多项式运算在MATLAB中，多项式可以很方便的进行运算。函数conv进行乘法运算，deconv进行除法运算。MATLAB没有提供特别的多项式加减法运算，因为这和向量的加减法很相似。运算举例如下：
a=[1 2 3];b=[1 2];
p=conv(a,b) %乘法
c =1 4 7 6
[q,r]=deconv(c,b) %除法
q = 1 2 3 %商多项式
r = 0 0 0 0 %余数多项式
需要注意的是，多项式除法并不一定能够除尽，很多时候需要有余数多项式。
MATLAB还提供了多项式微分、估计值函数。多项式微分使用polyder(p)函数，估计值使用polyval(p,at)函数。
三、数据分析
1、极小化数值分析中，很多应用需要确定函数的极值，即最大值或者最小值。在数学上，可通过确定函数导数为零的点解析的求出极值。但是在很多时候，很难找到导数为零的点，这样就难以解析的求极值。必须通过树枝上寻找函数的极值。MATLAB提供了fmin和fmins两个函数来求极值，它们分别寻找一维和n维函数的极值。它们的应用比较相似，以fmin为例看如何求极值。
为了求取一维函数的极值，需要给出函数的定义和极值范围。比如求函数f(x)=10*exp(-x)*cos(x)的极小值，可以使用下面的方法：
x='10*exp(-x)*cos(x)';
plot(fx,[2,5])
min=fmin(fx,2,5)
xmin =2.35619746669214
图3
上面的fplot命令画给定函数的曲线，例子中的函数曲线为图3-3。fplot命名在搜索极值的过程中，不断的计算函数值，如果函数的计算量比较大，或者搜索区域内有多个极值，那么搜索的过程可能比较长，在一些情况下，也可能找不到极值。如果fplot找不到极值，即停止运行并且提供解释。
Fmins也可以寻找极大值点，只要重新定义函数为-f(x)即可。
fmins函数也是搜索极值，但它搜索向量的标量函数最小值，使用的单纯性法搜索最小值。
在MATLAB的优化工具(Optimization Toolbox)箱中，有更多的扩展的优化算法。
2、求零点寻找函数值过零或者等于某一个常数值也是十分重要的问题，比如在使用bode图判断控制系统稳定性时，需要看幅频特性过零点和相频特性过1800点。一般使用解析算法求解这一问题十分困难，很多时候还是不可能的。MATLAB提供了该问题的数值算法。函数fzero可以寻找一维函数的过零点。例如：
zero=fzero(fx,5)
xzero =
4.71238898038469
zero=fzero(fx,2)
xzero =
1.57079632679490
这表示f(x)函数具有两个过零点，从图中也可以看出。
Fzero函数不但可以寻找过零点，也可以寻找函数值等于常值点，只要重新定于函数为f(x)-c即可。
3、积分
MATLAB提供了三个函数计算函数在有限区域内的积分：trapz、quad和quad8。函数trapz通过计算梯形面积的和近似函数的积分，函数的分割是人为地。例如：
p=1.5:0.01:5;
r=10*exp(-x).*cos(x);
rea=trapz(x,y)
area =-1.07578023141031
quad使用Simpson递归方法，quad8使用Newton-costes递归方法进行数值积分。为了获得更精确的结果，它们在所需的区间都计算被积函数。quad8比quad更精确。这两个函数的使用和fzero相同。
4、微分与积分相反，数值微分十分困难。积分描述了一个函数的整体或宏观的性质，所以积分对函数的形状在小范围的变化不敏感；而微分则描述了函数在一点处的斜率，是函数的微观性质，它对函数的微小变化十分敏感，函数的很小的变化，容易产生相邻点斜率的巨大变化。由于数值微分的固有困难，应当尽量避免使用数值微分，尤其是树言数据获得的数据微分。如果迫切需要，最好先将试验数据进行最小二乘拟合伙这三次样条拟合，然后对拟合函数进行微分。
MATLAB提供了一个有限插分函数diff，可以做数值微分。在应用中需要注意的是，插分后输出数组比原数组少了一个元素。实际表明，使用有限插分近似将放大噪声，导致极差的结果。
5、FFT变换
FFT即快速傅立叶变换，是数据分析的基本方法，它是向量X的离散傅立叶变换由基2的快速变换算法来计算的。如果X的长度不是精确的2次幂则后面使用0填充，lfft(x)是向量x的离散傅立叶变换的逆变换。
在频率轴商会值FFT曲线，要明确FFT结果与实际频率点的关系。设N个数据点，采样频率为fs，则Nyquist频率或n=N/2+1点与实际频率的关系为：
f=(num-1)*fs/n
例如：绘制向量x的频谱：
y=fft(x);
n=length(y)/2; %FFT变换是对称的，取前半部分
y=y(1:n);
f=fs*(0:n-1)/n;
plot(f,abs(y));
figure;
plot(f,(180/pi)*unwrap(atan2(imag(y),real(y))));
需要注意的是fft结果为复数矩阵，为了得到幅频特性，可使用abs函数，使用atan2得到相角，由于有的系统的相角可能大于1800，而相角函数值域在-1800~1800之间，需要使用unwrap函数展开折叠的相角，从而得到相频特性。
6、数据分析函数
MATLAB提供了很多有用的数据分析函数，这些函数数量极大，需要慢慢发掘。下面给出一些常用的函数：
函数名
含义
max
最大值
min
最小值
mean
均值
std
标准方差
sedian
中值
sort
排序分类
sum
元素的总和
prod
元素的乘积
cumrod
元素的累积
cumsum
元素的累加和
Diff
差分函数
Hist
直方图
Tabel
列表
Corr
互相关矩阵
Cov
协方差矩阵
Find
查找逻辑
这些函数可以灵活应用，比如在生成Bode图的频率序列时，使用logspace直接得到了对数等间距的序列，如果想插入几点自己需要的频率点，可以使用下面的命令：
p=logspace(0,2,100);
r=[1.5 3.4 56 78];
t=sort([x f]);
这样生成的频率点集合xf中不但含有1到100的对数等间隔频率点，还包含了频率f中的四个频率点，这个频率在绘制有尖锐幅频特性的bode图时十分有用。
四、常微分方程数值解控制系统的模型很大一部分是常微分方程形式的，求取它们的解析解比较困难，一般使用数值解。常微分方程数值解一般使用逐步积分的方法实现，Runge-Kutta法是应用最多的一种微分方程数值解的方法。MATLAB提供了两种Runge-Kutta法函数：
[t,x]=ode23(fun?t0,tf,x0,tol,trace)
[t,x]=ioe45(fun?,t0,tf,x0,tol,trace)
这两种方法格式相同。其中xfun为定义的常微分方程函数名，该函数必须以为输出，以t、x为输入。输入变量t0、tf为积分的启始和中止时间，单位是秒。x0为初始的状态向量。tol控制结果的精度，可以缺省。一般来说，ode45比ode23运算速度快一些。
考虑描述经典的Var der Pol微分方程
求解中将这高阶方程等价变换为一阶方程组，重新定义变量，变换为：
令 x1=x x2=dx/dt
则 dx1/dt=x2
dx2/dt=u(1-x12)x2-x1
把它写成函数：
function yp=vdp(t,x)
yp(1)=x(2);
yp(2)=2*(1-x(1)^2)*x(2)-x(1)； %令u=2
在命令行求解这个方程：
[t,x]=ode45(dp?0,20,[1 ;1]_);
plot(t,x(:,1),t,x(:,2)); %画出x和dx/dt的时域波形+图 x和dx/dt的时域波形
第二部分：
一 数字信号处理与MATLAB语言基础
1.1简介
MATLAB是一套用于科学工程计算的可视化高性能语言与软件环境.它集数值分析,矩阵运算,信号处理和图形显示于一体,构成一个界面友好的用户环境,在这个环境中,问题与求解都能方便地以数学的语言(主要是矩阵形式)或图形方式表示出来.
MATLAB的含义是矩阵实验室(Matrix Laboratory),该软件是一个交互式系统,其基本元素是无需定义维数的矩阵萦研制的初衷主要是方便矩阵的存取,但以过几十年的扩充和完善,已成为和各类科学研究与工程应用中的标准工具,其典型应用有:数值分析,算法分析,自动控制,数字信号处理,图像处理以及模型仿真等等.
MATLAB包括了被称为Toolbox(工具箱)的各类应用问题的求解工具,本书着重介绍其中的信号处理工具箱.MATLAB信号处理工具箱包含了各类经典的和现代的数字信号处理技术,是一个非常优秀的算法研究与辅助设计工具,它在语音信号处理,生物医学工程,实时控制和雷达信号处理等多个研究领域得到成功的应用.
本章将简要地介绍MATLAB的基本知识,为后面章节中信号处理的理解和学习打好基础,
本书的所有介绍都以MATLAB6.1版本为基础,如使用其他版本,可能有所不同.
1.2.1 MATLAB快速入门
1.2.1 启动MATLAB有两种方法
(1)在WINDOW95/98/2000/NT4.0下,单击任务条上的"开始"按钮,选择"程序"子菜单,再选择"MATLAB"程序组,最后单击"MATLAB6.1"菜单项,就可以启动MATLAB,由于在安装了MATLAB软件后,会自动在桌面上添加一个图标"MATLAB6.1",故也可以直接双击该图标以MATLAB.
(2)在WINDOWS95/98/2000/NT4.0下,单击任务条上"开始"按钮,选择"运行"菜单项,然后在"运行"对话框内输入"MATLAB"命令,然后点击"确定"按钮,就可以MATLAB.
以上两种操作的结果都会出现一个MATLAB的命令窗口,其命令提示符为"?".
1.2.2 MATLAB命令窗口与基本的矩阵操作
MATLAB命令窗口中,在MATLAB提示符下要输入一个4*4矩阵可以按如下 方式输入命令:
X=[1 2 3 4; 5 4 3 2 ;3 4 5 6;7 6 5 4];
或
X=[1 2 3 4
5 4 3 2
3 4 5 6
7 6 5 4 ];
以上两种输入方式的效果是一样的,命令未尾的分号用于禁止显示该命令的执行结果,如果去掉分号,MATLAB就会显示如下的命令的执行结果:
X=[1 2 3 4;5 4 3 2;3 4 5 6;7 6 5 4]
X=
1 2 3 4
5 4 3 2
3 4 5 6
7 6 5 4
要察看X变量的内容,也可以直接输入法X,按回车,MATLAB的显示与上面的结果似,如果要引用X的某几行(以第二,三行为例),则
X([2,3],:)
ans=
5 4 3 2
3 4 5 6
如果要引用X的某几列(以第二,三,四列为例),则
x(：，2：4)
ans=
2 3 4
4 3 2
4 5 6
6 5 4
要表示一个递增右递减的序列,可以用冒号,其格式为
n：s：m
此式表示产生从n到m且步长为s的一系列值,当s省略时,默认步长为1,在前面的取x的列的例子中我们已经用到了该方法,现在,如果要倒过来取列为4,3,2,且行为3,2,1的值,那么命令如下:
x(3：-1：1,4：-1：2)
ans= 6 5 4
2 3 4
4 3 2
求一个矩阵的转置,命令如下:
x'
ans=
1 5 3 7
2 4 4 6
3 3 5 5
4 2 6 4
要求矩阵X的逆,命令如下(本例中的X的矩阵已为奇异阵,求逆的结果已无意义):
y=inv(x)
MATLAB除了基本的标量运算外,还有各种矩阵运算与数组运算.
矩阵与矩阵,矩阵与标量的加法,减法,和乘法在形式上与标量运算类似,只是当矩阵与矩阵作相应的运算时,要求维数相互匹配,有时需要对矩阵进行必要的转置操作,这属于线性代数的基础知识.
MATLAAB为矩阵的除法提供了两种运算,左除(\)和右除(/),如果a为一非奇异矩阵,则a\b和b\a分别等价于:
a/b=inv(a)*b
b/a=b*inv(a)
矩阵乘方的运算符为^,它分为以下几种情况:
(1)当a为方阵,p为大于1的整数,则a^p为a的p次方,即a自乘p次.
(2)当a为方阵,p为非整数时,则a^p涉及到 a的特征向量v和特征值矩阵d,即有ap=v*d.^p/v.
(3)当a为标题而为矩阵时,a^p为计算标量的矩阵幂
(4)当a和p均为矩阵时,不能计算a^p.
数组运算实际上是对应元素的标题运算,其基本运算在表示上是在相应的标题运算符相加一个点"."(加减法不用加),表示对应数组的元素作相应的标题运算,见下面的几个例子,
x=[1 2 3];
y=[4 5 6];
x+y
ans=
5 7 9
x-y
ans=
-3 -3 -3
x.*y
ans=
-4 10 18
x./y
ans=
0.2500 0.4000 0.5000
x.\y
ans=
4.000 2.5000 2.0000
x.^y
ans=
1 32 729
对于MATLAB语言,需要注意是,它是大小写敏感的语言,即变量Abc和abc是两个不同的变量,函数各也不例外,如inv函数,如果写成 INV函数,则成了未定义的函数.
1.2.3 命令行编辑器
在命令行中输入命令时,可以利用功能键方便地修改当前或以前的命令行,以提高命令输入的效率.]
命令行的功能键如表1-1所示.
按 键
功 能
按 键
功 能
↑,Ctrl-P
重新调出上一行
Home,Ctrl-A
光标移到行首
↓,Ctrl-N
重新调出下一行
End,Ctrl-E
光标移到行尾
←,Ctrl-B
光标左移一个字符
Esc
清除命令行
→,Ctrl-F
光标右移一个字符
Del,Ctrl-D
删除光标处的字符
Ctrl-←
光标左移一个字
Backspace
删除光标左处的字符
Ctrl-→
光标左移一个字
Ctrl-K
删除至行尾
1.2.4 MATLAB帮助系统
要学会使用MATLAB,首先要学会的两个命令是,help和 demo.这两个函数构成了基本的MATLAB的帮助,help命令的使用如下(以查看 demo的使用方法为例):
help demo
DEMO Run demonstrations
Type"demo"at the command linetObrowse available demos.
With the optional action argument demo
('matlab'|'toolbox'|'simulink'|'blockset'|'statefiow'),
DEMO opens the demo screen[0 the specified subtopic]
With the optional categoryArg argument,
DEMO opens IO the specified toolbox or category,e.g.
demo toolbox signal
demo matlab language
只要有一定的英文阅读能力,掌握了help,也就相当于掌握了整个MATLAB.至于 demo,则是用于演示MATLAB在各个应用领域中的基本应用,可以作为MATLAB初学者的范例,例如,如果要学习信号处理工具箱的使用,可以键入下面的两条指令:
help signal
demo toolbox signal
MATLAB中的另一类基于图形界面的帮助系统,可以能过访问命令窗口的HELP菜单中的菜单项.常用的有"HelpWindow","Help Desk(HTML)","Examplesand Demos".其作用与前面介绍的两条命令类似,但界面更友好.
1.2.5 MATLAB的搜索路径
MATLAB利用自身的搜索路径来寻找M文件函数,如果你要执行的文件不在搜索路径中,就无法执行.利用path命令可显示当前的MATLAB搜索路径,如下,
path
MATLABPATH
E:\MATLABR11\toolbox\matlab\general
E:\MATLABR11\toolbox\matlab\ops
E:\MATLABR11\toolbox\matlab\lang
E:\MATLABR11\toolbox\matlab\elmat
E:\MATLABR11\toolbox\matlab\elfun
E:\MATLABR11\toolbox\matlab\specfun
E:\MATLABR11\toolbox\matlab\matfun
E:\MATLABR11\toolbox\matlab\datafun
E:\MATLABR11\toolbox\matlab\polyfun
E:\MATLABR11\toolbox\matlab\funfun
E:\MATLABR11\toolbox\matlab\sparfun
E:\MATLABR11\toolbox\matlab\graph2d
E:\MATLABR11\toolbox\matlab\graph3d
E:\MATLABR11\toolbox\matlab\specgraph
E:\MATLABR11\toolbox\matlab\graph9cs
E:\MATLABR11\toolbox\matlab\uitools
E:\MATLABR11\toolbox\matlab\stffun
E:\MATLABR11\toolbox\matlab\iofun
E:\MATLABR11\toolbox\matlab\timefun
E:\MATLABR11\toolbox\matlab\datatypes
E:\MATLABR11\toolbox\matlab\winfun
E:\MATLABR11\toolbox\matlab\demos
......
E:\MATLABR11\toolbox\signal\siggui
E:\MATLABR11\toolbox\signal\sigdemos
E:\MATLATR11\work
E:\MATLABR11\toolbox\locsl
要将指定的路径加入,可以使用path,addpath命令:要删除某个路径,可以使用rmpath命令:要修改路径,可以使用pathtool命令,为了方便,可以只使用pathtool,因为该命令执行后,会弹出一个图形界面,十分易于使用,其他命令的具体使用可以用help命令查看.
1.2.5 MATLAB图形系统
MATLABR提供了强大图形系统功能,本书并不打算详细地介绍该功能,其实,本书中最常用的函数是plot函数.
其具体的使用方法可以用help命令获得,这里以实例说明基本的使用,实例的运行结果如图1-1所示
x=0:0.02:20
y=cos(x);
z=sin(x);
plot(x,y,x,z,':');
1.3 离散时间信号
信号是传递信息的载体,按照信号的特点不同,可将信号表示成一个或几个独立变量的函数。
通常信号可以分为以下几类:
连续信号:在连续时间范围内定义的信号,信号的幅值可以是连续数值,也可以是离散数值
离散信号:时间为离散变量的信号,即独立变量时间被量化,它只在离散时间上给出函数值,是时间上不连续的"序列".离散时间的间隔是均匀的,以了表示,Z(nT)表示此离散时间信号在nT 点上的值,一般直接用.x(n)表示离散时间信号.
1.3.1 信号表示
MATLAB 中的主要数据类型是二维或多维的实矩阵或复矩阵.数字信号处理过程中所用到的基本数据对象(例如:一维信号或序列,多通道信号,二维信号等等)都可以用矩阵来表示.
MATLAB一般把普通的一维抽样数据信号即抽样序列表示成向量形式,向量表示为 1*n或 n*1的矩阵,其中n为序列中抽样点的个数
最简单的把序列引入MATLAB的方法是在命令行中输入一个元素表,例如:
x=[4 3 7 -9 1]
这样就构造了一个表示成行向量的五元素简单实数序列,当然也可转换成列向量式:
x=x'
结果为
x=
4
3
7
-9
1
列向量形式常用于表示单通道信号,因为它可以很自然地扩展到多通道的情况,对多通道的信号矩阵,矩阵的第一列代表一个通道,第一行代表一个抽样点.
一个包含x,2x和x/pi 的三通道信号可表示如下:
y=[x 2x x/pi]
结果为
y=
4.000 8.000 1.2732
3.000 6.000 0.9549
7.000 14.000 2.2282
-9.000 -18.0000 -2.8648
1.000 2.0000 0.3183
1.3.2 波形发生器
许多不同的工具箱函数都可以产生波形,其中大部分函数都需要一个时间向量作为参数,如果选择1000Hz的抽样频率产生波形,则适宜的时间向量如下:
t= (0:0.001:1);
这样构造了一个有1001个元素的行向量,该向量表示时间从0到1秒,以千分之一秒为步长.
1,基本信号序列
(1)单位抽样序列
这一序列可用MATLAB中的zeros 函数实现:
x=[1 zeros(1,n-1)];
(2)单位阶跃矩阵
这一序列可以利用MATLAB的ones函数实现:
x=ones(1,N);
(3)实指数序列
MATLAB实现:
n=0:N-1;
x=a.^n;
(4)复指数序列
MATLAB实现
n=0:N-1;
x=exp((lu+j*w0)*n);
(5)随机序列
MATLAB提供了两种随机信号:
Rand(1,N)产生[0,1]上均匀分布的随机矢量.
Rand(1,N)产生均值为0,方差为1的高斯随机序列,即白噪声序列.
2.基本周期波形
(1)方波
MATLAB工具函数square可以产生图1-2所示的方波
t=0:0.0*pi:5*pi
y=square(t);
plot(t,y);
axis([0 6*pi -1.5 1.5]);
xlabel('时间');
ylabel('幅值');
(2)正弦波
sin函数可以产生正弦信号(如图1-3所示).例如;
t=0:0.01*pi:2*pi;
y=sin(3*pi*t);
plot(t,y);
xlabel('时间t');
ylabel('幅值y');
(3)锯齿波工具箱函数sawtooth函数可以产生锯齿波(图1-4)
MATLAB实现如下;
Fs=1000;
t=0:1/Fs:1; %抽样长度为1s,抽样频率为1000Hz
y=sawtooth(2*pi*100*t); %信号频率为100Hz
plot(t,y);
axis([0 0.1 -1 1]); %画出了0.1S的波形
3 基本非周期波形
工具箱函数chirp能够产生一种扫射频率信号,其特点是信号的瞬时频率随时间按照一定规律变化,图1-5
t=0:1/1000:2; %抽样频率为1KHz,抽样时间为2s
x=chirp(t,0.1,200); %0时刻为DC信号,1s时频率为200Hz
specgram (x,256,1000,256,250);
4 sinc信号
MATLAB实现(1-6)
t=linspace(-6,6);
y=sinc(t);
plot(t,y);
1.3.3 序列的操作
(1)信号加 x(n)=x1(n)+x2(n)
MATLAB实现:
x=x1+x2;
注意:x1+x2序列应该具有相同的长度,位置对应,才能相加.
(2)信号乘 x(n)=x1(n)x2(n)
MATLAB实现:
x=x1.*x2;%数组乘法
(3)改变比例
y=k*x(n)
MATLAB实现:
y=k*x;
(4)折叠 y(n)=x(-n);
MATLAB实现
y=flip(x);
(5)抽样和
MATLAB实现:
y=sum(x(n1:n2));
(6)抽样积
MATLAB实现:
y=prod(x(n1:n2));
(7)信号能量
MATLAB实现:Ex=sum(abs(x).^2);
(8)信号功率
MATLAB实现:Px=sum(abs(x).^2)/N;
例如,以下程序的功能是产生一个叠加信号,并计算它的抽样和,抽样积,信号能量功率等,练习一下序列的各种操作.
程序如下(1-7)
t=0:0.001:1;
x1=sin(2*pi*50*t)+2*sin(2*pi*120*t);
%x1为两个正弦信号的和,一个为50Hz,幅值1;另一为120Hz,幅值2.
x=x1+0.5*randn(size(t));%采用size(t)统一了两个序列的长度
n=length(x);
y1=sum(x(1:50));%抽样和
y2=sum(x(1:50));%抽样积
Ex=sum(abs(x).^2);%信号能量
Px=sum(abs(x).^2)/n,%信号功率
Subplot(2,3,1);
plot(t(1:50),x1(1:50));
title('X1');
subplot(2,3,4);
plot(t(1:50),x(1:50));
title('x');
subplot(2,3,2);
plot(y1,'x');
title('抽样和');
subplot(2,3,3);
plot(y2,'x');
title('能量');
subplot(2,3,6);
plot(Px,'x');
title('功率');
二 卷积及其在matlab中的实现
2.1 卷积的定义对于线性时不变系统,卷积和可表示为:
上式右边对运算称为x[n]和h[n]对卷积运算.用符号记作:
y[n]=x[n]*h[n]
卷积的性质:
交换率性质:a(n)*b(n)=b(n)*a(n)
分配率性质:a(n)*(b(n)+c(n))=a(n)*b(n)+a(n)*c(n)
结合率性质:a(n)*(b(n)*c(n))=(a(n)*b(n))*c(n)
2.2 卷积在matlab中的实现
在MATLAB中,要直接使用其几个内部的函数conv求出两序列的卷积,其格式为:
y=conv(a,b)
该函数默认序列从n=0开始,如果序列是从一负值开始,就不能直接采用conv函数..但可以编写一个新的卷积函数来求这种情况下的卷积
本章节中采用一个新编写的一个卷积函数convzz来实现卷积计算,该函数采用基于"序列阵表格法"的卷积算法
已求序列x1(n)=[1 1 1 -1 -1]与x2(n)=[1 1 1 1 1 1]的卷积为例,如下表所示,x1(n)从左到右排列,x2(n)从上到下排列,表中填入的值是x1(n)与x2(n) 的元素得到的,沿各条对角虚线将各乘积相加,即得到各卷积s(n)之值.s(n)是的s(0)项是x1(0)与x2(0)相交所在的对角线上各项乘积的相加值,由此往右下推移,得到s(1)s(2)....显然表左上角的第一项就是s(0),其值为1,通过简单的运算,就可以得到卷积和
s(n)=x1(n)*x2(n)
=[1 2 3 2 1 0 -1 -2 -3 -2 -1]
若所求的卷积和s(n)不是从 s(0)开始,那么在表中,从s(0)向上推移,即可依次得到s(-1),s(-2),s(-3)…
2.3 卷积函数convzz的应用举例
[ 例1] 求出序列x1=[1 1 2 3 4 5]与x2=[2 4 1 3 4 5]的卷积和。
本例采用新编写的卷积函数convzz求解序列x1和x2的卷积和，程序如下所示：
x1=[1 1 2 3 4 5 ];
x2=[2 4 1 3 4 5];
y=convzz(x1,x2); %计算x1he x2卷积和
la=length(x1); %获得x1，x2的长度
lb=length(x2);
subplot(3,1,1); %画出的x1曲线图
t=0:1:la-1;
stem(t,x1);
title('the x1 sequence '); %标出曲线的标题
subplot(3,1,2); %画出x2的曲线图
t=0:1:lb-1;
stem(t,x2);
title('the x2 sequence ');
subplot(3,1,3); 画出卷积和的曲线图
t=0:1:(la+lb－2);
stem(t,y);
title('the convolution result ');
[例2] 方波信号和三角波信号的卷积
为计算这两个信号对卷积，本文在程序中定义了一个c(n)序列表示方波信号序列，d(n)序列表示三角波信号序列，具体如下
m=0.1; %m=0.1,T=2,
T=2; %m表示采样间隔，T表示周期
K=1; %定义方波信号序列
for t=0:m:T
c(k)=1;
k=k+1;
end
for t=T+m:m:2﹡T
c(k)=0;
k=k+1;
end
n=1;
for t=0:m:2﹡T
d(n)=10﹡t;
n=n+1;
end
for t=2﹡T+m:m,4﹡T
d(n)=0;
n=n+1;
end
y=convzz(c,d); %对两信号进行卷积
lc=length(c);
ld=length(d);
subplot(3,1,1); %画出结果曲线
t=0:1:lc－1;
stem(t,c);
title(‘the square sequence’); %写出图的标题
subplot(3,1,2);
t=0:1:ld－1;
stem(t,d);
title(‘the triangle sequence’);
subplot(3,1,3);
t=0:1:(lc+ld－2);
stem(t,y);
title(‘the convolution result’);
[例3] 指数信号和方波的卷积
为计算这两个信号对卷积，定义一个c(n)序列表示指数信号，d(n)序列表示方波信号，指数信号的底为a=1.5。matlab程序参考如下：
a=1.5;m=7;n=5; %a=1.5,m=7,n=5,
%a 表示幂函数的底，m 表示指数信号采样点%%个数，n表示方波信号的采样点个数
for I= 1:1:m %定义指数信号序列
c(i)a^(I－1);
end
I=m+2:1:2﹡m;
C(i)=0;
j=1:1:n; %定义方波信号序列
d(j)=1;
j=n+1:1:2*n;
d(j)=0;
y=convzz(c,d); %对两信号进行卷积
lc=length (c);
ld=length (d);
subplot(3,1,1);
t=0:1:l c-1; %画出结果曲线图
stem(t,c);
title(‘the exponent sequence’); %写出图的标题
subplot(3,1,2);
t=0:1:ld-1;
title(‘the exponent sequence’);
subplot(3,1,3);
t=0:1:(lc+ld-2);
stem(t,y);
title(‘the convolution result’);
[例4] 用matlab验证卷积的结合律性质结合律,x(n)*h1(n)*h2(n)=[x(n)*h1(n)]*h2(n)
=x(n)*[h1(n)*h2(n)]=[x(n)*h2(n)]*h1(n)
x=[1 2 3 4 5];
h1=[2 2 1 3 4 1];
h2=[2 4 5 3 4 5 4];
x1=convzz(x,h1); %x(n)*h1(n)*h2(n)
y1=convzz(x1,h2);
x2=convzz(h1,h2);
y2=convzz(x,x2); %x(n)*[h1(n)*h2(n)]
y1=Columns 1 through 12
4 20 56 118 219 350 484 600 698 728 650 516
Columns 13 through 16
405 276 121 20
y2=
Columns 1 through 12
4 20 56 118 219 350 484 600 698 728 650 516
Columns 13 thgough 16
405 276 121 20
由上可见,y1=y2,亦即:x(n)*h1(n)*h2(n)=x(n)*[h1(n)*h2(n)],其他情形读自行证明.
2.4 练习
（1）.求出两序列a=[1 1 1 1 1 1 1 1],b=[2 1 3 1 2 3 1 1 2]的卷积和，并画出图形曲线。
（2） 求出一正弦信号和方波信号的卷积和，并画出图形曲线。
（3） 请用matlab证明卷积的交换律和分配律性质。
（4） 试根据俯录中卷积的算法，编写函数实现序列的卷积从某一负值开始。
三 傅立叶变换及其性质
3.1傅立叶变换定义与性质
傅立叶变换就是建立在以时间为自变量的”信号”与以频率为自变量的”频谱函数”之间的某种关系.
对于N点离散序列,其傅立叶变化定义为:
N
X(k)=sum x(n)*exp(-j*2*pi*(k-1)*(n-1)/N),1<=k<=N.
n=1
对N点离散序列,其傅立叶反变换定义为:
N
X(n)=(1/N)sum x(k)*exp(-j*2*pi*(k-1)*(n-1)/N),1<=n<=N.
k=1
3.2傅立叶变换在matlab中的实现
FFT是快速傅立叶变换的简写,它是离散傅立叶变换(DFT)的一种快速算法,在matlab中可直接利用函数fft进行傅立叶正变换,由ifft函数可进行傅立叶反变换.
本书中采用上述定义编写一函数fftzz,实现傅立叶变换,另一函数ifftzz实现傅立叶反变换.两函数的源程序见附录2.1和2.2。
3.3傅立叶变换在matlab中的应用
[例 1] 用函数fftzz实现对序列的a=[1 1 1 1 1 1 1 1]的傅立叶变换
a=[1 1 1 1 1 1 1 1];
p=16;
y=fftzz(a,p); %为频率取样点个数
t=0:1:p-1; %对序列进行傅立叶变换
subplot(2,1,2)
stem(t,abs(y)) %画出输出结果幅度曲线图
subplot(2,1,2)
stem(t,angle(y)) %画出输出结果相位曲线图
【例2】离散时间数字信号傅立叶变换时域扩展性质的MATLAB实现。
令k是一个正整数，定义
x[n/k],当n为k的整数倍
x(k)[n]=
0,当n不为k的整数倍那么这时的x(k)[n]是在x [n]的连续值之间插入（k-1）个零值得到的。直观上看，可以把x(k)[n]看作是减慢了的x [n]。因为，除非n是k的某一倍数，也即n＝rk，否则x(k)[n]都等于0，所以的傅立叶变换可由下式给出：
由于x(k)[rk]＝x [r]，可求得：
由上可知，当k>1时，该信号在时间上被拉开了，从而在时间上就减慢了，而它的傅立叶变换就受到压缩。下面一个矩型脉冲信号序列为例子来说明这一性质。程序如下（结果参见图3-2）。由图上可以看出，随着k增加，在时域上拉开，而其变换则在频域上压缩。
A1＝[111100000000];
A2＝[1010101000000000];
A3＝[10010010010000000000];
P=200;
T1=lenglth(A1);
T2=lenglth(A2);
T3=lenglth(A3);
Y1=-fftzz(A1,P);
Y2=-fftzz(A2,P);
Y3=-fftzz(A3,P);
T=0:1:P-1;
SUBPLOT(3,2,2);
PLOT(T,Y1),GRID ON;
SUBPLOT(3,2,4);
PLOT(T,Y2),GRID ON;
SUBPLOT(3,2,6);
PLOT(T,Y3),GRID ON;
SUBPLOT(3,2,1);
T1=0;T1-1;
STEM(T1,A1);
SUBPLOT(3,2,3);
T2=0;T2-1;
STEM(T2,A2);
T3=0;T3-1;
SUBPLOT(3,2,5);
STEM(T3,A3);
【例3】用matlab验证傅立叶变换的性质：两序列时域卷积的傅立叶变换等于两序列傅立叶变换的乘积。亦即对x[n]*y[n]的傅立叶变换等于X(jw)Y(jw),X(jw),Y(jw)分别为x[n]、y[n]的傅立叶变换。
先对两个序列进行卷积，然后对卷积结果进行傅立叶变换，保存结果。然后再分别对两个序列求傅立叶变换，将两个序列傅立叶变换对结果相乘。将前后两个运算结果画图进行比较。两序列时域卷积后的傅立叶变换对结果参见图3-3，两序列傅立叶变换和的相乘结果参见图3-4。从两图中可以看出，运算结果是相等的。
C=[1 2 3 4 5 6];
D=[6 5 4 3 2 1];
A=CONVZZ(C,D); ％对两个序列C、D进行卷积运算
P=LENGTH(A);
G=FFTZZ(A,P); ％然后进行傅立叶变换
A=C;
E=FFTZZ(A,P); ％对两个序列先进行傅立叶变换
A=D;
F=FFTZZ(A,P);
H=E,*F; ％将两序列傅立叶变换的结果相乘
T=0:1:P-1; ％将上述两种方法对计算结果画图比较
SUBPLOT(2,1,1)
STEM(T,ABS(G))
TITLE(‘时域卷积后的傅立叶变换的幅度’);
SUBPLOT(2,1,2)
STEM(T,ANGLE(G))
TITLE(‘时域卷积后的傅立叶变换的相位’);
FIGURE(2);
T=0:1:P-1;
SUBPLOT(2,1,1)
STEM(T,ABS(H))
TITLE(‘频域相乘后的幅度’);
SUBPLOT(2,1,2)
STEM(T,ANGLE(H))
TITLE(‘频域相乘后的相位’);
3.4练习
（1）.设定一信号并计算它的fft变换
（2）.模拟信号x(t)=2*sin(4*pi*t)+5*cos(8*pi*t),求N点DFT的幅值谱和相位谱
（3）.比较原序列与经过FFT和IFFT变换后的序列
四,频谱搬移和调制
4.1 频谱搬移和调制的定义
所谓调制，就是利用信号来控制载波的某一参数，使这个参数随该信号而变化。幅度调制是应用较为广泛的一种调制方式，其原理就是使载波的振幅随调制信号的规律而变化，从而达到携带信号进行传输的目的。假设调制信号是普通正弦波：，载波信号为：，那么已调波的数学表达式为：
由上可知，调制后的信号由三个不同频率组成，第一项为未调制的载波；第二项为载波频率与调制频率之和，称为上边频；第三项为载波频率与调制频率之差，称为下边频。显然，调幅过程实际上就是一种频谱搬移过程，经过调制后，调制信号的频谱被搬移到载频附近，成为上边频和下边频。
4.2 频谱搬移和调制在matlab中的实现
编写一个函数lxfft，以实现连续傅立叶变换，其源程序如下：
function y=lxfft(ti,a,p) %ti表示采样间隔，a表示采样后的序列，p为
N=length(a); %频率取样点
for k=0:1:p
y(k+1)=0;
for I=1:1:N
y(k+1)=a(i)*exp(-1*j*(k)*(I-1)*ti)*ti+y(k+1);
end
end
y; %画出变换后的频谱图
I=0:1:p
stem(I,abs(y));
其中，ti表示取样间隔，a表示取样后的序列，p表示频率取样点。先将需实现频谱搬移或调制的信号进行取样，然后经过傅立叶变换，即可看到信号调制前后的频谱关系。
4.3 频谱搬移和调制的应用举例
[例 1] 求正弦信号a=sin（6*t）的频谱。
先将正弦信号进行抽样，然后将抽样的序列进行连续傅立叶变换，即可得出该信号的频谱图。频谱曲线参见图4—2，由图可见，该正弦信号只有单一频率成份。
ti=0.01; %序列的取样间隔
t=0:ti:42; %取样范围
a=sin(6*t); %求出正弦信号的取样序列
p=29; %傅立叶变换频率取样点个数
y=lxfft(ti,a,p); %对其求连续傅立叶变换
[例2] 信号的调制及频谱搬移的实现
已知调制信号为y=sin(t),载波信号为y=5sin(20*t)。
那么已调波的表达式为y=(5+2sin(t))*sin(20*t)，用下面的matlab程序实现，调制波形参见图4-3。
ti=0.01;
t=0:ti:30;
y=(5+2*sin(t))*sin(20*t);
stem(t,y);
频谱搬移的MATLAB实现：
已知调制信号为x1=sin(6*t)，载波信号为x2=sin(20*t)。已调波的表达式为：
y=(5+sin(6*t))*sin(20*t)
先将调制信号和载波信号分别进行连续时间的傅立叶变换，可看到调制前两信号的频谱参见图4-4，调制之后信号频谱曲线参见图4-5。由上图可见信号经过调制后，调制信号的频率分量被搬移到载波频率分量到两边。
ti=0.01;
t=0:ti:100;
p=50;
x1=sin(6*t);
x2=sin(20*t);
subplot(2,1,1);
y1=1xfft(ti,x1,p);
title(’调制信号的频谱’);
subplot(2,1,2);
y2=1xfft(ti,x2,p);
title(‘载波信号的频谱’);
figure(2);
y=(5+sin(6*t))*sin(20*t);
z=1xfft(ti,y,p);
title(‘已调波的频谱’);
4.4练习画出两余弦信号cos(6*t)与cos(50*t)调制前后的频谱图，以及调制后的调幅波。
五,波形分解
5.1波形分解和合成的定义任何信号都是由各种不同频率，幅度和初相的正弦波叠加而成的。对于周期信号，由它的傅氏计数展开式可知，各次谐波的频率是基波频率的整数倍。而非周期信号则包含了从零到无穷大的所有频率成分，每一频率成分的幅度均趋于无限小，但其相对大小则是不同的。利用频率选择滤波器，可以把各个不同频率的谐波分量提取出来，频率选择滤波就是一种专门设计成基本上无失真地通过某些频率，而显著地衰减掉另一些频率的系统。将信号中各个频率分量提取出来的过程就叫做波形分解。而分解出来的各个频率分量的波形相加，就是信号波形的重新合成。
5.2波形分解在matlab中的实现首先，必须先设计滤波器，在matlab中提供了丰富的滤波设计函数，其中cheby2称为切比雪夫II型滤波器设计，格式为：
[b,a]=cheby2(m,Rs,Wn)
其中n为该滤波器的阶数，表示阻带内的纹波（dB），Wn为滤波器的通带。0=<Wn<=1，其中1对应于抽样频率的1/2（Nyquist频率）。当Wn＝【W1 W2】时，cheby2函数产生一个2n阶的数字带通滤波器，通带为W1与W2之间。该函数可以设计模拟或数字的低通、高通、带通、和带阻滤波器。
举例：（1）设抽样频率为1000Hz，要设计一个10阶的低通切比雪夫II型滤波器，截至频率为200Hz，且阻带内的波纹系数为40dB
[b,a]=cheby2(m,Rs,Wn)
（2）设计一个10阶的带通切比雪夫II型滤波器，通带为200-300Hz，抽样频率为1000Hz，阻带纹波系数为40dB
n=10/2
Wn=[200 300]/500;
[b,a]=cheby2(n,40,Wn)
滤波函数filtfilt（b，a，x）利用设计好的滤波器对数据进行滤波。filtfilt（b，a，x）为零位数字滤波格式为：
y＝filtfilt（b，a，x）
通过将输入数据前向和反向处理，以完成零相位数字滤波。它先将数据按顺序滤波，然后将所得结果逆转后反向通过滤波器，这样得到的序列为精确零相位失真，并使滤波器的阶数加倍。
5.3波形分解与合成应用举例对频率为5000Hz的对称周期方波进行分解，并把分解后的波形重新合成，比较分解后合成的信号与源信号的差别。
首先，设计滤波器及滤波函数bandpass（w，x）。参数w表示通带的两个截止频率，x为需要滤波的输入序列。然后利用方波函数square生成一个频率为5000Hz的周期方波，利用设计好的滤波函数bandpass（w，x）对其进行滤波，各次谐波分别用不同颜色来区别，为便于比较，将波形时间轴选为一个方波周期。在波形分解对同时，将分解后各个谐波分量对数据分别保存，然后将这些分量相加，可以比较分解后再合成对波形与原波型对区别。波形分解后各个谐波分量对波型参见图5-1，分解后重新合成对波形参见图5-2
function y=bandpass(w,x)
Wn=w/200000;
[b,a]=cheby2(8/2,40,Wn); %设计一个8阶，阻带内纹波系数为40dB的cheby2II型滤波器
y=filtfilt(b,a,x); %对输入数据进行零相位数字滤波
Fs=400000;%抽样频率
t=0:1/Fs:1;
x=square(2*pi*5000*t); %利用方波产生函数square产生方波
plot(t,x),axis([0.00095 0.00125 –1.5 1.5]),gird on; %画出方波图形曲线
hold on;
for n=1:2:9%利用设计好的滤波器提取各个频率分量
w=[5000*n-3000 5000*n+3000];
y=bandpass(w,x);
switch n
case 1
m=’r’; ％当频率分量不同时，采用不同的线条画图，以示区别
y1=y; ％将滤波后的数据进行保存，以便波形合成时使用
case 3
m=’b’;
y2=y;
case 5
m=’g’;
y3=y;
case 7
m=’y’;
y4=y;
case 9
m=’k’;
y5=y;
end
plot(t,y,m),axis([0.00095 0.00125 –1.5 1.5]);
hold on;
end
Y=y1+y2+y3+y4+y5; %将分解后的谐波分量重新合成
plot(t,y),axis([0.00095 0.00125 –1.5 1.5])； ％画出合成后的波形
5.4练习对周期三角波进行分解，并将分解后的谐波分量重新合成，将分解前的波形和分解后再合成的波形进行比较。
附录表附录1.1
卷积函数convzz
function y=convzz(c,d)
1e=length(c); %获得两个输入序列的长度
1f=length(d);
if 1e<1f %将较长的序列从左到右排列，较短到从
%上到下排列
a=d;
b=c;
else
a=c;
b=d;
end
1a=length(a);
1b=length(b);
if 1a~=prod(size(a)) | 1b~=prod(size(b))
error(A and B must be vectors.’);
end
y=zeros(1,1a+1b-1); %定义一个长度为1a+1b-1到全零序列，用于存
%放卷积结果
for k=1:1b; %第一转折点前点循环运算
i=k; j=1;
while i>0
y(k)=a(i)*b(j)+y(k);
i=i-1; j=j+1;
end
end
for k=(1b+1):1a %第二转折点前的循环运算]
i=k; j=1;
while i>k-1b
y(k)=a(i)*b(j)+y(k);
i=i-1; j=j+1;
end
end
for k=(1a+1):(1a+1b-1) %第二转折点后的循环运算
i=1a; j=k-1a+1;
while i>k-1b
y(k)=a(i)*b(j)+y(k);
i=i-1; j=j+1;
end
end
附录2.1
离散傅立叶变换fftzz及其反变换ifftzz函数
离散傅立叶变换fftzz 函数源程序：
%离散序列的傅立叶变换
function y=fftzz(a,p) %a表示输入的序列，p表示频率取样点个数
N=length(a); %获得a序列的长度
y=zeros(1,p); %定义一个长度为p的全零序列
%外循环
for k=1:1:p
%内循环
for n=1:1:N
y(k)=a(n)*exp(-j*2*pi*(k-1)*(n-1)/p)+y(k); %在某一频率点上对
end %全部序列值运算
if abs(y(k))<10^ (-6) %当y(k)的模小于某一值时，将其置零
y(k)=0; %避免画相位图时出现不确定性
end
end
附录2.2
离散序列傅立叶反变换ifftzz函数源程序
%离散序列傅立叶反变换
function y=ifftzz(b,p)
N=length(b);
Y=zeros(1,p);
for k=1:1:N
for n=1:1:p
y(k)=(1/p)*b(n)*exp(j*2*pi*(k-1)*(n-1)/p)+y(k);
if abs(y(k))<10^(-6)
y(k)=0;
end
end
end
三、实验内容部分实验一 序列的相加和相乘
一、实验内容给出两个序列x (n)和x (n).
x=[0,1,2,3,4,3,2,1,0];n1=[-2:6];
x=[2,2,0,0,0,-2,-2];n=[2:8];
现在求它们的和y及乘积y.
二、实验目的巩固序列的相乘和相加运算。
对MARLAB在同位项相乘和详加中的应用。
三、运行程序参考结果
实验二 序列的移位和周期延拓运算
实验内容已知x(n)=0.8R (n),利用MATLAB生成并图示x(n),x(n-m),x((n))R(n)x((n))表示x(n)以8为周期的延拓和x((n))R(n)，其中N=24,m为一个整常数，0<m<N.
实验目的
1.
2.
三、参考结果图
实验三 信号的采样
实验内容令X(t)=e ,求出并绘制其傅立叶变换X(j)。用三个不同的采样，分别求出并画出离散时间傅立叶变换X(e).三个频率分别为：
(1) f=5kHz; (2) f=2kHz; (3) f=1kHz.
实验目的
1.
2.
三,参考结果图
实验四 基本序列的离散傅立叶变换计算一、实验内容用复正弦序列x (n)=
余弦序列x (n)=
正弦序列x (n)=
分别对N=16和N=8计算序列的N点DFT，并绘制幅频曲线，最后用DFT理解为何两
种N值下的DFT结果差别如此之大。
实验目的
1.
2.
三,参考结果图
实验五 用Z变换求解差分方程一、实验内容令y(n)-2y(n-1)+3y(n-2)=u(n)+ux(n-1)+5u(n-2)-6u(n-3),
初始条件为x(-1)=1,x(-2)=1,y(-1)=-1,y(-2)=1 。
求y(n)。
二、实验目的
1.
2.
三、参考结果图
实验六 快速傅立叶变换
一、实验内容用快速卷积法计算下面两个序列的卷积。并测试直接卷积和快速卷积的时间。
x(n)=
h(n)=
二、实验目的
1.
2.
参考结果图
实验七 用DFT对连续信号作谱分析一、实验内容
用DFT分析的频谱结构。选择不同的截取长度，观察DFT进行频谱分析时存在的截断效应（频谱泄露和谱间干扰）。试用加窗的方法减少谱间干扰。
二、实验目的
1.
2.
三、参考结果图
实验八 用各种窗函数设计FIR数字滤波器实验内容分别用矩形窗和Hamming窗设计线性相位FIR低通滤波器。要求通带截止频率=，单位脉冲响应h(a)的长度N=21。绘制及其幅频响应特性曲线。
实验目的
1.
2.
参考结果图
,
实验九 切比雪夫型低通数字滤波器设计实验内容设计一个切比雪夫型低通数字滤波器设计，指标如下：
通带边界频率：,通带最大衰减：
阻带截止频率： 阻带最小衰减：。
实验目的
1.
2.
参考结果图
实验参考程序一、解：（1）MATLAB程序实现如下：
x1=[0,1,2,3,4,3,2,1,0];ns1=-2; %给定x1及ns1
x2=[2,2,0,0,0,-2,-2];ns2=2; %给定x2及ns2
nf1=ns1+length(x1)-1;nf2=ns2+length(x2)-1;
ny=min(ns1,ns2):max(nf1,nf2); %y(n)的时间变量
xa1=zeros(1,length(ny));xa2=xa1; %延拓序列初始化
xa1(find((ny>=ns1)&(ny<=nf1)==1))=x1; %给xa1赋值x1
xa2(find((ny>=ns2)&(ny<=nf2)==1))=x2; %给xa2赋值x2
ya=xa1+xa2 %序列相加
yp=xa1.*xa2 %序列相乘
subplot(2,2,1),stem(ny,xa1,'.') %绘图
subplot(2,2,2),stem(ny,xa2,'.')
line([ny(1),ny(end)],[0,0]) %画x轴
subplot(2,2,3),stem(ny,ya,'.')
line([ny(1),ny(end)],[0,0])
subplot(2,2,4),stem(ny,yp,'.')
line([ny(1),ny(end)],[0,0])
（2）程序运行结果如图所示：（略）
（3）结论为：
ny=-2 –1 0 1 2 3 4 5 6 7 8
xa1=0 1 2 3 4 3 2 1 0 0 0
xa2=0 0 0 0 2 2 0 0 0 -2 -2
ya= 0 1 2 3 6 5 2 1 0 -2 -2
yp= 0 0 0 0 8 6 0 0 0 0 0
二、解：（1）MATLAB程序实现如下：
clear;close all
N=24;M=8;
m=3;
if (m<1|m>=N-M+1) %检验输入参数n是否合理
fprintf;
end
n=0:N-1;
x1=(0.8).^n; x2=[(n>=0)&(n<M)]; %产生x(n)
xn=x1.*x2;
xm=zeros(1,N); %设定xm的长度
for k=m+1:m+M
xm(k)=xn(k-m);
end
xc=xn(mod(n,8)+1); %产生x(n)的周期延拓
xcm=xn(mod(n-m,8)+1); %产生x(n)移位后的周期延拓
subplot(2,2,1),stem(n,xn,'.')
subplot(2,2,2),stem(n,xm,'.')
subplot(2,2,3),stem(n,xc,'.')
subplot(2,2,4),stem(n,xcm,'.')
三、解：（1）MATLAB程序实现如下：
Dt=0.00005;t=-0.005:Dt:0.005;xa=exp(-1000*abs(t)); %模拟信号
Wmax=2*pi*2000;K=500;k=0:1:K;W=k*Wmax/K;
Xa=xa*exp(-j*t'*W)*Dt;Xa=real(Xa); %连续时间傅立叶变换
W=[-fliplr(W),W(2:501)]; %频率从-Wmax to Wmax
Xa=[fliplr(Xa),Xa(2:501)]; %Xa介于-Wmax和Wmax间
subplot(421);plot(t*1000,xa);
ylabel('xa(t)');
subplot(422);plot(W/(2*pi*1000),Xa*1000);
ylabel('Xa(jw)');
Ts=0.0002;n=-25:1:25;x=exp(-1000*abs(n*Ts)); %离散信号
K=500;k=0:1:K;w=pi*k/K;
X=x*exp(-j*n'*w);X=real(X); %离散时间傅立叶变换
w=[-fliplr(w),w(2:K+1)];
X=[fliplr(X),X(2:K+1)];
subplot(423);stem(n*Ts*1000,x);
ylabel('xl(n)');
subplot(424);plot(w/pi,X);
ylabel('Xl(w)');
Ts=0.0005;n=-10:1:10;x=exp(-1000*abs(n*Ts)); %离散信号
K=500;k=0:1:K;w=pi*k/K;
X=x*exp(-j*n'*w);X=real(X); %离散时间傅立叶变换
w=[-fliplr(w),w(2:K+1)];
X=[fliplr(X),X(2:K+1)];
subplot(425);stem(n*Ts*1000,x);
ylabel('xl(n)');
subplot(426);plot(w/pi,X);
ylabel('Xl(w)');
Ts=0.0010;n=-5:1:5;x=exp(-1000*abs(n*Ts)); %离散信号
K=500;k=0:1:K;w=pi*k/K;
X=x*exp(-j*n'*w);X=real(X); %离散时间傅立叶变换
w=[-fliplr(w),w(2:K+1)];
X=[fliplr(X),X(2:K+1)];
subplot(427);stem(n*Ts*1000,x);
ylabel('xl(n)');
subplot(428);plot(w/pi,X);
ylabel('Xl(w)');
四、解：（1）MATLAB程序实现如下：
clear; close all
N=16; N1=8;
n=0:N-1; k=0:N1-1;
x1n=exp(j*pi*n/8); %产生x1(n)
X1k=fft(x1n,N); %计算N点DFT[x1(n)]
Xk1=fft(x1n,N1); %计算N1点DFT[x1(n)]
x2n=cos(pi*n/8);
X2k=fft(x2n,N); %计算N点DFT[x2(n)]
Xk2=fft(x2n,N1); %计算N1点DFT[x1(n)]
x3n=sin(pi*n/8);
X3k=fft(x3n,N); %计算N点DFT[x3(n)]
Xk3=fft(x3n,N1); %计算N1点DFT[x1(n)]
subplot(2,3,1),
stem(n,abs(X1k),'.'),
subplot(2,3,2),
stem(n,abs(X2k),'.'),
subplot(2,3,3),
stem(n,abs(X3k),'.'),
subplot(2,3,3),
stem(n,abs(X1k),'.'),
subplot(2,3,4),
stem(k,abs(Xk1),'.'),
subplot(2,3,5),
stem(k,abs(Xk2),'.'),
subplot(2,3,6),
stem(k,abs(Xk3),'.')
五、解：（1）MATLAB的程序运行如下：
b=[1 4 5 -6];
a=[1 -2 3];
x0=[1 1];
y0=[-1 1];
xic=filtic(b,a,y0,x0)
bxplus=1;
axplus=[1 -1];
ayplus=conv(a,axplus)
byplus=conv(b,bxplus)+conv(xic,axplus)
[R,P,K]=residuez(byplus,ayplus)
Mp=abs(P)
Ap=angle(P)*180/pi
N=100;
n=0:N-1;
xn=ones(1,N);
yn=filter(b,a,xn,xic);
plot(n,yn)
（2）程序结果如图所示：
（略）
（3）结果为：
xic = 4 2 -6
ayplus = 1 -3 5 -3
byplus = 5 2 -3 0
R = 1.5000 - 4.2426i
1.5000 + 4.2426i
2.0000
P = 1.0000 + 1.4142i
1.0000 - 1.4142i
1.0000
K = 0
Mp = 1.7321
1.7321
1.0000
Ap = 54.7356
-54.7356
0
六、解：（1）MATLAB程序实现如下：
clear;close all
xn=sin(0.4*[1:15]');
hn=0.9.^[1:20]';
M=length(xn);N=length(hn);
nx=1:M;nh=1:N;
%循环卷积等于线性卷积的条件：循环卷积区间长度L>=M+N-1
L=pow2(nextpow2(M+N-1)); %取L为大于等于且最接近(N+M-1)的2的正次幂
tic,%快速卷积计时开始
Xk=fft(hn,L); %L点FFT[x(n)]
Hk=fft(hn,L); %L点FFT[h(n)]
Yk=Xk.*Hk; %频域相乘得Y(k)
yn=ifft(Yk,L); %L点IFFT得到卷积结果y(n)
toc %快速卷积计时结束
subplot(2,2,1),stem(nx,xn,'.');
ylabel('x(n)')
subplot(2,2,2),stem(nh,hn,'.');ylabel('h(n)')
subplot(2,1,2);ny=1:L;
stem(ny,real(yn),'.');ylabel('y(n)')
tic,yn=conv(xn,hn);toc %直接调用函数conv计算卷积与快速卷积比较
实验结果图如下略
（3）结论为：
elapsed_time = 0.0600
elapsed_time = 0
七、解,(1) Matlab程序为：
clear; close all
fs=400; T=1/fs; %采样视频率为400Hz
Tp=0.04;N=Tp*fs; %采样点数N
N1=[N,4*N,8*N]; %设定三种截取长度供调用
st=['|X1(jf)|';'|X4(jf)|';'|X8(jf)|']; %设定三种标注语句供调用
%矩形窗截断
for m=1:3
n=1:N1(m);
%产生采样序列x(n)
xn=cos(200*pi*n*T)+sin(100*pi*n*T)+cos(50*pi*n*T);
Xk=fft(xn,4096); %4096点DFT，用FFT实现
fk=[0:4095]/4096/T;
subplot(3,2,2*m-1)
plot(fk,abs(Xk)/max(abs(Xk)));ylabel(st(m,:))
if m==1 title(' ');end
end
%加hamming窗改善谱间干扰
for m=1:3
n=1:N1(m);
wn=hamming(N1(m)); %调用工具箱函数hamming产生N长hamming窗序列wn
xn=(cos(200*pi*n*T)+sin(100*pi*n*T)+cos(50*pi*n*T)).*wn';
Xk=fft(xn,4096); %4096点DFT，用FFT实现
fk=[0:4095]/4096/T;
subplot(3,2,2*m)
plot(fk,abs(Xk)/max(abs(Xk)));ylabel(st(m,:))
if m==1 title(' ');end
end
八、解,(1) Mtlab的程序为:
clear;close all %理想低通滤波器参数
N=21;wc=pi/4;
n=0:N-1;r=(N-1)/2;
hdn=sin(wc*(n-r))/pi./(n-r); %计算理想低通单位脉冲响应
if rem(N,2)~=0 hdn(r+1)=wc/pi;end %N为奇数时，处理n=r点的0/0型
wn1=boxcar(N); %矩形窗
hn1=hdn.*wn1'; %加窗
%以上两条语句可代以fir函数：hn1=fir1(N-1,wc/pi,boxcar(N));
subplot(2,1,1),
stem(n,hdn,'.'),hw=fft(hdn,512),w=2*[0:511]/512,subplot(2,1,2),plot(w,20*log10(abs(hw))
九、解,(1) Matlab的程序为:
clear;close all
wp=0.2;ws=0.4;Rp=1;Rs=80; %输入指标
[N,wc]=cheb2ord(wp,ws,Rp,Rs) %求滤波器阶次
[B,A]=cheby2(N,Rs,wc) %设计滤波器，得出系数
freqz(B,A) %无左端变量时自动画频率特性图
输出滤波参数如下：
N = 8
Wc = 0.3682
B = Columns 1 through 7
0.0011 0.0008 0.0026 0.0024 0.0032 0.0024 0.0026
Columns 8 through 9
0.0008 0.0011
A = Columns 1 through 7
1.0000 -4.3881 8.9294 -10.8479 8.5205 -4.4046 1.4577
Columns 1 through 7
参考文献
1.MATLAB信号处理详解 陈亚勇 等编著 人民邮电出版社
2.MATLAB 5.X程序设计语言 楼顺天 陈生潭 雷虎民 编著 西安交通大学
出版社
3.信号与系统 ALAN V.OPPENHEIM 刘树棠 译 西安交通大学出版社
4.数字信号处理—基于计算机的方法 Sanjit K.Mitra 清华大学出版社
5.信号与系统解题指南 胡光锐 徐昌庆 等 编著 科学出版社
6.
绪 论
MATLAB是一个功能十分强大的系统，是集数值计算、图形管理、程序开发为一体的环境。作为强大的科学计算平台，它几乎能满足所有的计算需求。在美国及其他发达国家的理工科院校里，MATLAB已成为了一门必修的课程，在科研院所、大型公司或企业的工程计算部门，MATLAB也是最为普遍的计算工具之一。
MATLAB具有如下的优势和特点：
1、好的工作平台和编程环境随着MATLAB的商业化以及软件本身的不断升级，MATLAB的用户界面也越来越精致，更加接近Windows的标准界面，人机交互性更强，操作更简单。而且新版本的MATLAB提供了完整的联机查询、帮助系统，极大地方便了用户的使用。简单的编程环境提供了比较完备的调试系统，程序不必经过编译就可以直接运行，而且能够及时地报告出现的错误及进行出错原因分析。
2、简单易用的程序语言新版本的MATLAB语言是基于最为流行的C语言基础上的，因此语法特征与C语言极为相似，而且更加简单，更加符合科技人员对数学表达式的书写格式。使之更利于非计算机专业的科技人员使用。而且这种可移植性好、可扩展性极强，这也是MATLAB之所以能够深入到科学研究及工程计算各个领域的重要原因。
3、强大的科学计算及数据处理能力
MATLAB拥有600多个工程中要用的数学运算函数，可以方便地实现用户所需的各种计算功能函数中所使用的算法都是科研和工程计算中的最新研究成果，而且经过了各种优化及容错处理，因此使用起来鲁棒性和可靠性非常高。在通常情况下，可以用它来代替底层编程语言，如C和C++。在计算要求相同的情况下，使用MATLAB的编程工作量会大大减少。MATLAB函数所能解决的问题包括矩阵运算和线性方程组的求解、微分及偏微分方程组的求解、符号运算、傅立叶变换和数据的统计分析、工程中的优化问题、稀疏矩阵运算、复数的各种运算、三角函数和其他初等数学运算、多维数组操作以及建模动态仿真等。
4、出色的图形处理功能
MATLAB自产生之日起就具有方便的数据可视化功能，新版本的MATLAB对整个图形处理的功能做了很大的改进和完善，使它不仅在一般数据可视化软件都具有的功能（例如二维曲线和三维曲面的绘制和处理等）方面更加完善，而且对于一些其他软件所没有的功能（例如图形的光照处理、色度处理以及四维数据的表现等），MATLAB同样表现了出色的处理能力。同时对一些特殊的可视化要求，例如图形动画等，MATLAB也有相应的功能函数，保证了用户不同层次的要求。另外，新版本的MATLAB还着重在图形用户界面（GUI）的制作上做了很大的改善，对这方面有特殊要求的用户也可以得到满足。
5、应用广泛的模块集或工具箱
MATLAB对许多专门的领域都开发了功能强大的模块集或工具箱。一般来说，它们都是由特定的专家开发的，用户可以直接使用工具箱学习、应用和评估不同的方法而不需要自由编写代码。目前，MATLAB已经把工具延伸到科学研究和工程应用的诸多领域，诸如数据采集、数据库接口、概率统计、样条拟合优化算法、偏微分方程求解、神经网络、小波分析、信号处理、图象处理、系统辩识、控制系统设计、LMI控制、鲁棒控制、模型预测、模糊逻辑、金融分析、地图工具、非线性控制设计、实时快速原型及半物理仿真、嵌入式系统开发、定点仿真、DSP与通讯、电力系统仿真等，都在工具箱（Toolbox）家族中有了自己的一席之地。
6、使用的程序接口和开发平台新版本的MATLAB可以利用MATLAB编译器和C/C++数学库和图形库，将自己的MATLAB程序自动转换为独立于MATLAB运行的C和C++代码。另外，MATLAB网页服务程序还容许在Web应用中使用自己的MATLAB数学和图形程序。
7、模块化的设计和系统级的仿真
Simulink是MATLAB的一个分支产品，主要用来实现对工程问题的模型化和动态仿真。在世界范围内的模型化浪潮的背景下Simulink恰恰体现了模块化设计和系统级仿真的具体思想，使得建模仿真如同搭积木一样简单。Simulink对仿真的实现可以应用于动力系统、信号控制、通信设计、金融财会及生物医学等各个领域的研究中。
一、基本概念
1、变量和其他高级语言一样，MATLAB也是使用变量来保存信息。变量由变量名表示，变量的命名应遵循如下规则：
1）变量名必须以字母开头。
2）变量名可以由字母、数字和下划线混合组成。
3）变量名区分字母大小写。
4）变量名的字符长度不应超过31个。
在MATLAB中还存在着一些系统默认的固定变量，如表2-1所示，即在MATLAB中语句中若出现固定变量名，则系统就将赋予默认值。
变量名
默认值
i
虚数单位
j
虚数单位
pi
圆周率
inf
无穷大
MATLAB的变量分为字符变量和数值变量两种，字符变量必须用但引号括起来。例如，用户可输入：
a=‘happly new year’
则表示将字符串 ‘happly new year’ 赋值给字符变量a.。
若用户输入：b=365
则表示将数值365赋值给数值变量b。
和其他高级语不同的是，MATLAB使用变量时不需要预先对变量类型进行说明，MATLAB会自动根据所输入的数据来决定变量的饿数据类型和分配存储空间。
2、数值在MATLAB内部，每个数据元素都是用双精度来表示和存储的，大约有16位有效数字。其数值有效范围约为。
但在进行数据输入输出时，MATLAB却可以用不同的格式。如果参加运算的每个元素均为整数，则MATLAB将用不加小数点的纯整数格式显示运算结果，否则，按默认的输出格式显示结果。MATLAB的默认格式为short格式，该格式显示运算结果为保留小数点后4位有效数字。用户可以通过format命令改变输出格式为long，以得到更多的有效数字（小数点后14位）。需要注意的是，数据输出格式的改变，并不影响该数据在MATLAB内部的存储精度。设置为short很long输出格式的命令分别为：
format short
format long
MATLAB通常用十进制数来表示一个数，亦可用科学计数法来表示一个数。另外，MATLAB还可以进行复数运算，复数可以由如下语句来产生：
C=a+i*b (或c=a+j*b) 将实部a虚部为b的复数赋值给复变量c
C=a*exp(j*b) (或C=a*exp(j*b)) 将模为a辐角为b的复数赋值给复数c
其中I、j是虚数单位。
以下是MATLAB各种合法的输入数据示例：
365 -18 0.000076 8.2467-10
2.67i 1.02E3i -7.6 -3.4502e+20
3、矩阵
矩阵是MALTAB进行数据处理和运算的基本元素。MATLAB的大部分运算或命令都是在矩阵运算的意义下执行的。我们通常意义上的数量（标量）在MATLAB系统中是作为
的矩阵来处理的，而仅有一行或一列的矩阵在MATLAB中称为向量。
4、数组
在MATLAB中，数组也是一个重要的概念，矩阵在某些情况下可视为二阶的数值型数组。但是在MATLAB中，数组和矩阵运算规则却有较大的区别。例如，两矩阵相乘和两数组相乘遵循的运算规则就是完全不相同的。
5、函数
MATLAB为用户提供了丰富且功能各异的函数，用户可以直接用这些函数来进行数据处理。函数由函数名和参数组成，函数调用的格式为：
函数名（参数）
例如，若在MATLAB的命令窗口输入命令：
a=sin(b)
则表示计算b的正弦值并将赋值给变量a.。
二、多项式处理
1、多项式表示多项式在MATLAB中使用降幂系数的行向量表示。例如：多项式x4-12x3+0x2+25x+116表示为：p=[1 -12 0 25 116]
需要注意的是，表示中需要包含零系数的项。按照这种形式，使用函数roots可以找出多项式等于零的根：r=roots(p)
r =11.74728287082272,2.70282074384101,-1.22505180733187 + 1.46720800872890i、-1.22505180733187,- 1.46720800872890I
按照MATLAB的规定，多项式使用行向量，根使用列向量。给出多项式的根，使用poly函数也可以构造出相应的多项式：
u=poly(r)
pp =1.0000 -12.0000 0 25.0000 116.0000
2、多项式运算在MATLAB中，多项式可以很方便的进行运算。函数conv进行乘法运算，deconv进行除法运算。MATLAB没有提供特别的多项式加减法运算，因为这和向量的加减法很相似。运算举例如下：
a=[1 2 3];b=[1 2];
p=conv(a,b) %乘法
c =1 4 7 6
[q,r]=deconv(c,b) %除法
q = 1 2 3 %商多项式
r = 0 0 0 0 %余数多项式
需要注意的是，多项式除法并不一定能够除尽，很多时候需要有余数多项式。
MATLAB还提供了多项式微分、估计值函数。多项式微分使用polyder(p)函数，估计值使用polyval(p,at)函数。
三、数据分析
1、极小化数值分析中，很多应用需要确定函数的极值，即最大值或者最小值。在数学上，可通过确定函数导数为零的点解析的求出极值。但是在很多时候，很难找到导数为零的点，这样就难以解析的求极值。必须通过树枝上寻找函数的极值。MATLAB提供了fmin和fmins两个函数来求极值，它们分别寻找一维和n维函数的极值。它们的应用比较相似，以fmin为例看如何求极值。
为了求取一维函数的极值，需要给出函数的定义和极值范围。比如求函数f(x)=10*exp(-x)*cos(x)的极小值，可以使用下面的方法：
x='10*exp(-x)*cos(x)';
plot(fx,[2,5])
min=fmin(fx,2,5)
xmin =2.35619746669214
图3
上面的fplot命令画给定函数的曲线，例子中的函数曲线为图3-3。fplot命名在搜索极值的过程中，不断的计算函数值，如果函数的计算量比较大，或者搜索区域内有多个极值，那么搜索的过程可能比较长，在一些情况下，也可能找不到极值。如果fplot找不到极值，即停止运行并且提供解释。
Fmins也可以寻找极大值点，只要重新定义函数为-f(x)即可。
fmins函数也是搜索极值，但它搜索向量的标量函数最小值，使用的单纯性法搜索最小值。
在MATLAB的优化工具(Optimization Toolbox)箱中，有更多的扩展的优化算法。
2、求零点寻找函数值过零或者等于某一个常数值也是十分重要的问题，比如在使用bode图判断控制系统稳定性时，需要看幅频特性过零点和相频特性过1800点。一般使用解析算法求解这一问题十分困难，很多时候还是不可能的。MATLAB提供了该问题的数值算法。函数fzero可以寻找一维函数的过零点。例如：
zero=fzero(fx,5)
xzero =
4.71238898038469
zero=fzero(fx,2)
xzero =
1.57079632679490
这表示f(x)函数具有两个过零点，从图中也可以看出。
Fzero函数不但可以寻找过零点，也可以寻找函数值等于常值点，只要重新定于函数为f(x)-c即可。
3、积分
MATLAB提供了三个函数计算函数在有限区域内的积分：trapz、quad和quad8。函数trapz通过计算梯形面积的和近似函数的积分，函数的分割是人为地。例如：
p=1.5:0.01:5;
r=10*exp(-x).*cos(x);
rea=trapz(x,y)
area =-1.07578023141031
quad使用Simpson递归方法，quad8使用Newton-costes递归方法进行数值积分。为了获得更精确的结果，它们在所需的区间都计算被积函数。quad8比quad更精确。这两个函数的使用和fzero相同。
4、微分与积分相反，数值微分十分困难。积分描述了一个函数的整体或宏观的性质，所以积分对函数的形状在小范围的变化不敏感；而微分则描述了函数在一点处的斜率，是函数的微观性质，它对函数的微小变化十分敏感，函数的很小的变化，容易产生相邻点斜率的巨大变化。由于数值微分的固有困难，应当尽量避免使用数值微分，尤其是树言数据获得的数据微分。如果迫切需要，最好先将试验数据进行最小二乘拟合伙这三次样条拟合，然后对拟合函数进行微分。
MATLAB提供了一个有限插分函数diff，可以做数值微分。在应用中需要注意的是，插分后输出数组比原数组少了一个元素。实际表明，使用有限插分近似将放大噪声，导致极差的结果。
5、FFT变换
FFT即快速傅立叶变换，是数据分析的基本方法，它是向量X的离散傅立叶变换由基2的快速变换算法来计算的。如果X的长度不是精确的2次幂则后面使用0填充，lfft(x)是向量x的离散傅立叶变换的逆变换。
在频率轴商会值FFT曲线，要明确FFT结果与实际频率点的关系。设N个数据点，采样频率为fs，则Nyquist频率或n=N/2+1点与实际频率的关系为：
f=(num-1)*fs/n
例如：绘制向量x的频谱：
y=fft(x);
n=length(y)/2; %FFT变换是对称的，取前半部分
y=y(1:n);
f=fs*(0:n-1)/n;
plot(f,abs(y));
figure;
plot(f,(180/pi)*unwrap(atan2(imag(y),real(y))));
需要注意的是fft结果为复数矩阵，为了得到幅频特性，可使用abs函数，使用atan2得到相角，由于有的系统的相角可能大于1800，而相角函数值域在-1800~1800之间，需要使用unwrap函数展开折叠的相角，从而得到相频特性。
6、数据分析函数
MATLAB提供了很多有用的数据分析函数，这些函数数量极大，需要慢慢发掘。下面给出一些常用的函数：
函数名
含义
max
最大值
min
最小值
mean
均值
std
标准方差
sedian
中值
sort
排序分类
sum
元素的总和
prod
元素的乘积
cumrod
元素的累积
cumsum
元素的累加和
Diff
差分函数
Hist
直方图
Tabel
列表
Corr
互相关矩阵
Cov
协方差矩阵
Find
查找逻辑
这些函数可以灵活应用，比如在生成Bode图的频率序列时，使用logspace直接得到了对数等间距的序列，如果想插入几点自己需要的频率点，可以使用下面的命令：
p=logspace(0,2,100);
r=[1.5 3.4 56 78];
t=sort([x f]);
这样生成的频率点集合xf中不但含有1到100的对数等间隔频率点，还包含了频率f中的四个频率点，这个频率在绘制有尖锐幅频特性的bode图时十分有用。
四、常微分方程数值解控制系统的模型很大一部分是常微分方程形式的，求取它们的解析解比较困难，一般使用数值解。常微分方程数值解一般使用逐步积分的方法实现，Runge-Kutta法是应用最多的一种微分方程数值解的方法。MATLAB提供了两种Runge-Kutta法函数：
[t,x]=ode23(fun?t0,tf,x0,tol,trace)
[t,x]=ioe45(fun?,t0,tf,x0,tol,trace)
这两种方法格式相同。其中xfun为定义的常微分方程函数名，该函数必须以为输出，以t、x为输入。输入变量t0、tf为积分的启始和中止时间，单位是秒。x0为初始的状态向量。tol控制结果的精度，可以缺省。一般来说，ode45比ode23运算速度快一些。
考虑描述经典的Var der Pol微分方程
求解中将这高阶方程等价变换为一阶方程组，重新定义变量，变换为：
令 x1=x x2=dx/dt
则 dx1/dt=x2
dx2/dt=u(1-x12)x2-x1
把它写成函数：
function yp=vdp(t,x)
yp(1)=x(2);
yp(2)=2*(1-x(1)^2)*x(2)-x(1)； %令u=2
在命令行求解这个方程：
[t,x]=ode45(dp?0,20,[1 ;1]_);
plot(t,x(:,1),t,x(:,2)); %画出x和dx/dt的时域波形+图 x和dx/dt的时域波形
第二部分：
一 数字信号处理与MATLAB语言基础
1.1简介
MATLAB是一套用于科学工程计算的可视化高性能语言与软件环境.它集数值分析,矩阵运算,信号处理和图形显示于一体,构成一个界面友好的用户环境,在这个环境中,问题与求解都能方便地以数学的语言(主要是矩阵形式)或图形方式表示出来.
MATLAB的含义是矩阵实验室(Matrix Laboratory),该软件是一个交互式系统,其基本元素是无需定义维数的矩阵萦研制的初衷主要是方便矩阵的存取,但以过几十年的扩充和完善,已成为和各类科学研究与工程应用中的标准工具,其典型应用有:数值分析,算法分析,自动控制,数字信号处理,图像处理以及模型仿真等等.
MATLAB包括了被称为Toolbox(工具箱)的各类应用问题的求解工具,本书着重介绍其中的信号处理工具箱.MATLAB信号处理工具箱包含了各类经典的和现代的数字信号处理技术,是一个非常优秀的算法研究与辅助设计工具,它在语音信号处理,生物医学工程,实时控制和雷达信号处理等多个研究领域得到成功的应用.
本章将简要地介绍MATLAB的基本知识,为后面章节中信号处理的理解和学习打好基础,
本书的所有介绍都以MATLAB6.1版本为基础,如使用其他版本,可能有所不同.
1.2.1 MATLAB快速入门
1.2.1 启动MATLAB有两种方法
(1)在WINDOW95/98/2000/NT4.0下,单击任务条上的"开始"按钮,选择"程序"子菜单,再选择"MATLAB"程序组,最后单击"MATLAB6.1"菜单项,就可以启动MATLAB,由于在安装了MATLAB软件后,会自动在桌面上添加一个图标"MATLAB6.1",故也可以直接双击该图标以MATLAB.
(2)在WINDOWS95/98/2000/NT4.0下,单击任务条上"开始"按钮,选择"运行"菜单项,然后在"运行"对话框内输入"MATLAB"命令,然后点击"确定"按钮,就可以MATLAB.
以上两种操作的结果都会出现一个MATLAB的命令窗口,其命令提示符为"?".
1.2.2 MATLAB命令窗口与基本的矩阵操作
MATLAB命令窗口中,在MATLAB提示符下要输入一个4*4矩阵可以按如下 方式输入命令:
X=[1 2 3 4; 5 4 3 2 ;3 4 5 6;7 6 5 4];
或
X=[1 2 3 4
5 4 3 2
3 4 5 6
7 6 5 4 ];
以上两种输入方式的效果是一样的,命令未尾的分号用于禁止显示该命令的执行结果,如果去掉分号,MATLAB就会显示如下的命令的执行结果:
X=[1 2 3 4;5 4 3 2;3 4 5 6;7 6 5 4]
X=
1 2 3 4
5 4 3 2
3 4 5 6
7 6 5 4
要察看X变量的内容,也可以直接输入法X,按回车,MATLAB的显示与上面的结果似,如果要引用X的某几行(以第二,三行为例),则
X([2,3],:)
ans=
5 4 3 2
3 4 5 6
如果要引用X的某几列(以第二,三,四列为例),则
x(：，2：4)
ans=
2 3 4
4 3 2
4 5 6
6 5 4
要表示一个递增右递减的序列,可以用冒号,其格式为
n：s：m
此式表示产生从n到m且步长为s的一系列值,当s省略时,默认步长为1,在前面的取x的列的例子中我们已经用到了该方法,现在,如果要倒过来取列为4,3,2,且行为3,2,1的值,那么命令如下:
x(3：-1：1,4：-1：2)
ans= 6 5 4
2 3 4
4 3 2
求一个矩阵的转置,命令如下:
x'
ans=
1 5 3 7
2 4 4 6
3 3 5 5
4 2 6 4
要求矩阵X的逆,命令如下(本例中的X的矩阵已为奇异阵,求逆的结果已无意义):
y=inv(x)
MATLAB除了基本的标量运算外,还有各种矩阵运算与数组运算.
矩阵与矩阵,矩阵与标量的加法,减法,和乘法在形式上与标量运算类似,只是当矩阵与矩阵作相应的运算时,要求维数相互匹配,有时需要对矩阵进行必要的转置操作,这属于线性代数的基础知识.
MATLAAB为矩阵的除法提供了两种运算,左除(\)和右除(/),如果a为一非奇异矩阵,则a\b和b\a分别等价于:
a/b=inv(a)*b
b/a=b*inv(a)
矩阵乘方的运算符为^,它分为以下几种情况:
(1)当a为方阵,p为大于1的整数,则a^p为a的p次方,即a自乘p次.
(2)当a为方阵,p为非整数时,则a^p涉及到 a的特征向量v和特征值矩阵d,即有ap=v*d.^p/v.
(3)当a为标题而为矩阵时,a^p为计算标量的矩阵幂
(4)当a和p均为矩阵时,不能计算a^p.
数组运算实际上是对应元素的标题运算,其基本运算在表示上是在相应的标题运算符相加一个点"."(加减法不用加),表示对应数组的元素作相应的标题运算,见下面的几个例子,
x=[1 2 3];
y=[4 5 6];
x+y
ans=
5 7 9
x-y
ans=
-3 -3 -3
x.*y
ans=
-4 10 18
x./y
ans=
0.2500 0.4000 0.5000
x.\y
ans=
4.000 2.5000 2.0000
x.^y
ans=
1 32 729
对于MATLAB语言,需要注意是,它是大小写敏感的语言,即变量Abc和abc是两个不同的变量,函数各也不例外,如inv函数,如果写成 INV函数,则成了未定义的函数.
1.2.3 命令行编辑器
在命令行中输入命令时,可以利用功能键方便地修改当前或以前的命令行,以提高命令输入的效率.]
命令行的功能键如表1-1所示.
按 键
功 能
按 键
功 能
↑,Ctrl-P
重新调出上一行
Home,Ctrl-A
光标移到行首
↓,Ctrl-N
重新调出下一行
End,Ctrl-E
光标移到行尾
←,Ctrl-B
光标左移一个字符
Esc
清除命令行
→,Ctrl-F
光标右移一个字符
Del,Ctrl-D
删除光标处的字符
Ctrl-←
光标左移一个字
Backspace
删除光标左处的字符
Ctrl-→
光标左移一个字
Ctrl-K
删除至行尾
1.2.4 MATLAB帮助系统
要学会使用MATLAB,首先要学会的两个命令是,help和 demo.这两个函数构成了基本的MATLAB的帮助,help命令的使用如下(以查看 demo的使用方法为例):
help demo
DEMO Run demonstrations
Type"demo"at the command linetObrowse available demos.
With the optional action argument demo
('matlab'|'toolbox'|'simulink'|'blockset'|'statefiow'),
DEMO opens the demo screen[0 the specified subtopic]
With the optional categoryArg argument,
DEMO opens IO the specified toolbox or category,e.g.
demo toolbox signal
demo matlab language
只要有一定的英文阅读能力,掌握了help,也就相当于掌握了整个MATLAB.至于 demo,则是用于演示MATLAB在各个应用领域中的基本应用,可以作为MATLAB初学者的范例,例如,如果要学习信号处理工具箱的使用,可以键入下面的两条指令:
help signal
demo toolbox signal
MATLAB中的另一类基于图形界面的帮助系统,可以能过访问命令窗口的HELP菜单中的菜单项.常用的有"HelpWindow","Help Desk(HTML)","Examplesand Demos".其作用与前面介绍的两条命令类似,但界面更友好.
1.2.5 MATLAB的搜索路径
MATLAB利用自身的搜索路径来寻找M文件函数,如果你要执行的文件不在搜索路径中,就无法执行.利用path命令可显示当前的MATLAB搜索路径,如下,
path
MATLABPATH
E:\MATLABR11\toolbox\matlab\general
E:\MATLABR11\toolbox\matlab\ops
E:\MATLABR11\toolbox\matlab\lang
E:\MATLABR11\toolbox\matlab\elmat
E:\MATLABR11\toolbox\matlab\elfun
E:\MATLABR11\toolbox\matlab\specfun
E:\MATLABR11\toolbox\matlab\matfun
E:\MATLABR11\toolbox\matlab\datafun
E:\MATLABR11\toolbox\matlab\polyfun
E:\MATLABR11\toolbox\matlab\funfun
E:\MATLABR11\toolbox\matlab\sparfun
E:\MATLABR11\toolbox\matlab\graph2d
E:\MATLABR11\toolbox\matlab\graph3d
E:\MATLABR11\toolbox\matlab\specgraph
E:\MATLABR11\toolbox\matlab\graph9cs
E:\MATLABR11\toolbox\matlab\uitools
E:\MATLABR11\toolbox\matlab\stffun
E:\MATLABR11\toolbox\matlab\iofun
E:\MATLABR11\toolbox\matlab\timefun
E:\MATLABR11\toolbox\matlab\datatypes
E:\MATLABR11\toolbox\matlab\winfun
E:\MATLABR11\toolbox\matlab\demos
......
E:\MATLABR11\toolbox\signal\siggui
E:\MATLABR11\toolbox\signal\sigdemos
E:\MATLATR11\work
E:\MATLABR11\toolbox\locsl
要将指定的路径加入,可以使用path,addpath命令:要删除某个路径,可以使用rmpath命令:要修改路径,可以使用pathtool命令,为了方便,可以只使用pathtool,因为该命令执行后,会弹出一个图形界面,十分易于使用,其他命令的具体使用可以用help命令查看.
1.2.5 MATLAB图形系统
MATLABR提供了强大图形系统功能,本书并不打算详细地介绍该功能,其实,本书中最常用的函数是plot函数.
其具体的使用方法可以用help命令获得,这里以实例说明基本的使用,实例的运行结果如图1-1所示
x=0:0.02:20
y=cos(x);
z=sin(x);
plot(x,y,x,z,':');
1.3 离散时间信号
信号是传递信息的载体,按照信号的特点不同,可将信号表示成一个或几个独立变量的函数。
通常信号可以分为以下几类:
连续信号:在连续时间范围内定义的信号,信号的幅值可以是连续数值,也可以是离散数值
离散信号:时间为离散变量的信号,即独立变量时间被量化,它只在离散时间上给出函数值,是时间上不连续的"序列".离散时间的间隔是均匀的,以了表示,Z(nT)表示此离散时间信号在nT 点上的值,一般直接用.x(n)表示离散时间信号.
1.3.1 信号表示
MATLAB 中的主要数据类型是二维或多维的实矩阵或复矩阵.数字信号处理过程中所用到的基本数据对象(例如:一维信号或序列,多通道信号,二维信号等等)都可以用矩阵来表示.
MATLAB一般把普通的一维抽样数据信号即抽样序列表示成向量形式,向量表示为 1*n或 n*1的矩阵,其中n为序列中抽样点的个数
最简单的把序列引入MATLAB的方法是在命令行中输入一个元素表,例如:
x=[4 3 7 -9 1]
这样就构造了一个表示成行向量的五元素简单实数序列,当然也可转换成列向量式:
x=x'
结果为
x=
4
3
7
-9
1
列向量形式常用于表示单通道信号,因为它可以很自然地扩展到多通道的情况,对多通道的信号矩阵,矩阵的第一列代表一个通道,第一行代表一个抽样点.
一个包含x,2x和x/pi 的三通道信号可表示如下:
y=[x 2x x/pi]
结果为
y=
4.000 8.000 1.2732
3.000 6.000 0.9549
7.000 14.000 2.2282
-9.000 -18.0000 -2.8648
1.000 2.0000 0.3183
1.3.2 波形发生器
许多不同的工具箱函数都可以产生波形,其中大部分函数都需要一个时间向量作为参数,如果选择1000Hz的抽样频率产生波形,则适宜的时间向量如下:
t= (0:0.001:1);
这样构造了一个有1001个元素的行向量,该向量表示时间从0到1秒,以千分之一秒为步长.
1,基本信号序列
(1)单位抽样序列
这一序列可用MATLAB中的zeros 函数实现:
x=[1 zeros(1,n-1)];
(2)单位阶跃矩阵
这一序列可以利用MATLAB的ones函数实现:
x=ones(1,N);
(3)实指数序列
MATLAB实现:
n=0:N-1;
x=a.^n;
(4)复指数序列
MATLAB实现
n=0:N-1;
x=exp((lu+j*w0)*n);
(5)随机序列
MATLAB提供了两种随机信号:
Rand(1,N)产生[0,1]上均匀分布的随机矢量.
Rand(1,N)产生均值为0,方差为1的高斯随机序列,即白噪声序列.
2.基本周期波形
(1)方波
MATLAB工具函数square可以产生图1-2所示的方波
t=0:0.0*pi:5*pi
y=square(t);
plot(t,y);
axis([0 6*pi -1.5 1.5]);
xlabel('时间');
ylabel('幅值');
(2)正弦波
sin函数可以产生正弦信号(如图1-3所示).例如;
t=0:0.01*pi:2*pi;
y=sin(3*pi*t);
plot(t,y);
xlabel('时间t');
ylabel('幅值y');
(3)锯齿波工具箱函数sawtooth函数可以产生锯齿波(图1-4)
MATLAB实现如下;
Fs=1000;
t=0:1/Fs:1; %抽样长度为1s,抽样频率为1000Hz
y=sawtooth(2*pi*100*t); %信号频率为100Hz
plot(t,y);
axis([0 0.1 -1 1]); %画出了0.1S的波形
3 基本非周期波形
工具箱函数chirp能够产生一种扫射频率信号,其特点是信号的瞬时频率随时间按照一定规律变化,图1-5
t=0:1/1000:2; %抽样频率为1KHz,抽样时间为2s
x=chirp(t,0.1,200); %0时刻为DC信号,1s时频率为200Hz
specgram (x,256,1000,256,250);
4 sinc信号
MATLAB实现(1-6)
t=linspace(-6,6);
y=sinc(t);
plot(t,y);
1.3.3 序列的操作
(1)信号加 x(n)=x1(n)+x2(n)
MATLAB实现:
x=x1+x2;
注意:x1+x2序列应该具有相同的长度,位置对应,才能相加.
(2)信号乘 x(n)=x1(n)x2(n)
MATLAB实现:
x=x1.*x2;%数组乘法
(3)改变比例
y=k*x(n)
MATLAB实现:
y=k*x;
(4)折叠 y(n)=x(-n);
MATLAB实现
y=flip(x);
(5)抽样和
MATLAB实现:
y=sum(x(n1:n2));
(6)抽样积
MATLAB实现:
y=prod(x(n1:n2));
(7)信号能量
MATLAB实现:Ex=sum(abs(x).^2);
(8)信号功率
MATLAB实现:Px=sum(abs(x).^2)/N;
例如,以下程序的功能是产生一个叠加信号,并计算它的抽样和,抽样积,信号能量功率等,练习一下序列的各种操作.
程序如下(1-7)
t=0:0.001:1;
x1=sin(2*pi*50*t)+2*sin(2*pi*120*t);
%x1为两个正弦信号的和,一个为50Hz,幅值1;另一为120Hz,幅值2.
x=x1+0.5*randn(size(t));%采用size(t)统一了两个序列的长度
n=length(x);
y1=sum(x(1:50));%抽样和
y2=sum(x(1:50));%抽样积
Ex=sum(abs(x).^2);%信号能量
Px=sum(abs(x).^2)/n,%信号功率
Subplot(2,3,1);
plot(t(1:50),x1(1:50));
title('X1');
subplot(2,3,4);
plot(t(1:50),x(1:50));
title('x');
subplot(2,3,2);
plot(y1,'x');
title('抽样和');
subplot(2,3,3);
plot(y2,'x');
title('能量');
subplot(2,3,6);
plot(Px,'x');
title('功率');
二 卷积及其在matlab中的实现
2.1 卷积的定义对于线性时不变系统,卷积和可表示为:
上式右边对运算称为x[n]和h[n]对卷积运算.用符号记作:
y[n]=x[n]*h[n]
卷积的性质:
交换率性质:a(n)*b(n)=b(n)*a(n)
分配率性质:a(n)*(b(n)+c(n))=a(n)*b(n)+a(n)*c(n)
结合率性质:a(n)*(b(n)*c(n))=(a(n)*b(n))*c(n)
2.2 卷积在matlab中的实现
在MATLAB中,要直接使用其几个内部的函数conv求出两序列的卷积,其格式为:
y=conv(a,b)
该函数默认序列从n=0开始,如果序列是从一负值开始,就不能直接采用conv函数..但可以编写一个新的卷积函数来求这种情况下的卷积
本章节中采用一个新编写的一个卷积函数convzz来实现卷积计算,该函数采用基于"序列阵表格法"的卷积算法
已求序列x1(n)=[1 1 1 -1 -1]与x2(n)=[1 1 1 1 1 1]的卷积为例,如下表所示,x1(n)从左到右排列,x2(n)从上到下排列,表中填入的值是x1(n)与x2(n) 的元素得到的,沿各条对角虚线将各乘积相加,即得到各卷积s(n)之值.s(n)是的s(0)项是x1(0)与x2(0)相交所在的对角线上各项乘积的相加值,由此往右下推移,得到s(1)s(2)....显然表左上角的第一项就是s(0),其值为1,通过简单的运算,就可以得到卷积和
s(n)=x1(n)*x2(n)
=[1 2 3 2 1 0 -1 -2 -3 -2 -1]
若所求的卷积和s(n)不是从 s(0)开始,那么在表中,从s(0)向上推移,即可依次得到s(-1),s(-2),s(-3)…
2.3 卷积函数convzz的应用举例
[ 例1] 求出序列x1=[1 1 2 3 4 5]与x2=[2 4 1 3 4 5]的卷积和。
本例采用新编写的卷积函数convzz求解序列x1和x2的卷积和，程序如下所示：
x1=[1 1 2 3 4 5 ];
x2=[2 4 1 3 4 5];
y=convzz(x1,x2); %计算x1he x2卷积和
la=length(x1); %获得x1，x2的长度
lb=length(x2);
subplot(3,1,1); %画出的x1曲线图
t=0:1:la-1;
stem(t,x1);
title('the x1 sequence '); %标出曲线的标题
subplot(3,1,2); %画出x2的曲线图
t=0:1:lb-1;
stem(t,x2);
title('the x2 sequence ');
subplot(3,1,3); 画出卷积和的曲线图
t=0:1:(la+lb－2);
stem(t,y);
title('the convolution result ');
[例2] 方波信号和三角波信号的卷积
为计算这两个信号对卷积，本文在程序中定义了一个c(n)序列表示方波信号序列，d(n)序列表示三角波信号序列，具体如下
m=0.1; %m=0.1,T=2,
T=2; %m表示采样间隔，T表示周期
K=1; %定义方波信号序列
for t=0:m:T
c(k)=1;
k=k+1;
end
for t=T+m:m:2﹡T
c(k)=0;
k=k+1;
end
n=1;
for t=0:m:2﹡T
d(n)=10﹡t;
n=n+1;
end
for t=2﹡T+m:m,4﹡T
d(n)=0;
n=n+1;
end
y=convzz(c,d); %对两信号进行卷积
lc=length(c);
ld=length(d);
subplot(3,1,1); %画出结果曲线
t=0:1:lc－1;
stem(t,c);
title(‘the square sequence’); %写出图的标题
subplot(3,1,2);
t=0:1:ld－1;
stem(t,d);
title(‘the triangle sequence’);
subplot(3,1,3);
t=0:1:(lc+ld－2);
stem(t,y);
title(‘the convolution result’);
[例3] 指数信号和方波的卷积
为计算这两个信号对卷积，定义一个c(n)序列表示指数信号，d(n)序列表示方波信号，指数信号的底为a=1.5。matlab程序参考如下：
a=1.5;m=7;n=5; %a=1.5,m=7,n=5,
%a 表示幂函数的底，m 表示指数信号采样点%%个数，n表示方波信号的采样点个数
for I= 1:1:m %定义指数信号序列
c(i)a^(I－1);
end
I=m+2:1:2﹡m;
C(i)=0;
j=1:1:n; %定义方波信号序列
d(j)=1;
j=n+1:1:2*n;
d(j)=0;
y=convzz(c,d); %对两信号进行卷积
lc=length (c);
ld=length (d);
subplot(3,1,1);
t=0:1:l c-1; %画出结果曲线图
stem(t,c);
title(‘the exponent sequence’); %写出图的标题
subplot(3,1,2);
t=0:1:ld-1;
title(‘the exponent sequence’);
subplot(3,1,3);
t=0:1:(lc+ld-2);
stem(t,y);
title(‘the convolution result’);
[例4] 用matlab验证卷积的结合律性质结合律,x(n)*h1(n)*h2(n)=[x(n)*h1(n)]*h2(n)
=x(n)*[h1(n)*h2(n)]=[x(n)*h2(n)]*h1(n)
x=[1 2 3 4 5];
h1=[2 2 1 3 4 1];
h2=[2 4 5 3 4 5 4];
x1=convzz(x,h1); %x(n)*h1(n)*h2(n)
y1=convzz(x1,h2);
x2=convzz(h1,h2);
y2=convzz(x,x2); %x(n)*[h1(n)*h2(n)]
y1=Columns 1 through 12
4 20 56 118 219 350 484 600 698 728 650 516
Columns 13 through 16
405 276 121 20
y2=
Columns 1 through 12
4 20 56 118 219 350 484 600 698 728 650 516
Columns 13 thgough 16
405 276 121 20
由上可见,y1=y2,亦即:x(n)*h1(n)*h2(n)=x(n)*[h1(n)*h2(n)],其他情形读自行证明.
2.4 练习
（1）.求出两序列a=[1 1 1 1 1 1 1 1],b=[2 1 3 1 2 3 1 1 2]的卷积和，并画出图形曲线。
（2） 求出一正弦信号和方波信号的卷积和，并画出图形曲线。
（3） 请用matlab证明卷积的交换律和分配律性质。
（4） 试根据俯录中卷积的算法，编写函数实现序列的卷积从某一负值开始。
三 傅立叶变换及其性质
3.1傅立叶变换定义与性质
傅立叶变换就是建立在以时间为自变量的”信号”与以频率为自变量的”频谱函数”之间的某种关系.
对于N点离散序列,其傅立叶变化定义为:
N
X(k)=sum x(n)*exp(-j*2*pi*(k-1)*(n-1)/N),1<=k<=N.
n=1
对N点离散序列,其傅立叶反变换定义为:
N
X(n)=(1/N)sum x(k)*exp(-j*2*pi*(k-1)*(n-1)/N),1<=n<=N.
k=1
3.2傅立叶变换在matlab中的实现
FFT是快速傅立叶变换的简写,它是离散傅立叶变换(DFT)的一种快速算法,在matlab中可直接利用函数fft进行傅立叶正变换,由ifft函数可进行傅立叶反变换.
本书中采用上述定义编写一函数fftzz,实现傅立叶变换,另一函数ifftzz实现傅立叶反变换.两函数的源程序见附录2.1和2.2。
3.3傅立叶变换在matlab中的应用
[例 1] 用函数fftzz实现对序列的a=[1 1 1 1 1 1 1 1]的傅立叶变换
a=[1 1 1 1 1 1 1 1];
p=16;
y=fftzz(a,p); %为频率取样点个数
t=0:1:p-1; %对序列进行傅立叶变换
subplot(2,1,2)
stem(t,abs(y)) %画出输出结果幅度曲线图
subplot(2,1,2)
stem(t,angle(y)) %画出输出结果相位曲线图
【例2】离散时间数字信号傅立叶变换时域扩展性质的MATLAB实现。
令k是一个正整数，定义
x[n/k],当n为k的整数倍
x(k)[n]=
0,当n不为k的整数倍那么这时的x(k)[n]是在x [n]的连续值之间插入（k-1）个零值得到的。直观上看，可以把x(k)[n]看作是减慢了的x [n]。因为，除非n是k的某一倍数，也即n＝rk，否则x(k)[n]都等于0，所以的傅立叶变换可由下式给出：
由于x(k)[rk]＝x [r]，可求得：
由上可知，当k>1时，该信号在时间上被拉开了，从而在时间上就减慢了，而它的傅立叶变换就受到压缩。下面一个矩型脉冲信号序列为例子来说明这一性质。程序如下（结果参见图3-2）。由图上可以看出，随着k增加，在时域上拉开，而其变换则在频域上压缩。
A1＝[111100000000];
A2＝[1010101000000000];
A3＝[10010010010000000000];
P=200;
T1=lenglth(A1);
T2=lenglth(A2);
T3=lenglth(A3);
Y1=-fftzz(A1,P);
Y2=-fftzz(A2,P);
Y3=-fftzz(A3,P);
T=0:1:P-1;
SUBPLOT(3,2,2);
PLOT(T,Y1),GRID ON;
SUBPLOT(3,2,4);
PLOT(T,Y2),GRID ON;
SUBPLOT(3,2,6);
PLOT(T,Y3),GRID ON;
SUBPLOT(3,2,1);
T1=0;T1-1;
STEM(T1,A1);
SUBPLOT(3,2,3);
T2=0;T2-1;
STEM(T2,A2);
T3=0;T3-1;
SUBPLOT(3,2,5);
STEM(T3,A3);
【例3】用matlab验证傅立叶变换的性质：两序列时域卷积的傅立叶变换等于两序列傅立叶变换的乘积。亦即对x[n]*y[n]的傅立叶变换等于X(jw)Y(jw),X(jw),Y(jw)分别为x[n]、y[n]的傅立叶变换。
先对两个序列进行卷积，然后对卷积结果进行傅立叶变换，保存结果。然后再分别对两个序列求傅立叶变换，将两个序列傅立叶变换对结果相乘。将前后两个运算结果画图进行比较。两序列时域卷积后的傅立叶变换对结果参见图3-3，两序列傅立叶变换和的相乘结果参见图3-4。从两图中可以看出，运算结果是相等的。
C=[1 2 3 4 5 6];
D=[6 5 4 3 2 1];
A=CONVZZ(C,D); ％对两个序列C、D进行卷积运算
P=LENGTH(A);
G=FFTZZ(A,P); ％然后进行傅立叶变换
A=C;
E=FFTZZ(A,P); ％对两个序列先进行傅立叶变换
A=D;
F=FFTZZ(A,P);
H=E,*F; ％将两序列傅立叶变换的结果相乘
T=0:1:P-1; ％将上述两种方法对计算结果画图比较
SUBPLOT(2,1,1)
STEM(T,ABS(G))
TITLE(‘时域卷积后的傅立叶变换的幅度’);
SUBPLOT(2,1,2)
STEM(T,ANGLE(G))
TITLE(‘时域卷积后的傅立叶变换的相位’);
FIGURE(2);
T=0:1:P-1;
SUBPLOT(2,1,1)
STEM(T,ABS(H))
TITLE(‘频域相乘后的幅度’);
SUBPLOT(2,1,2)
STEM(T,ANGLE(H))
TITLE(‘频域相乘后的相位’);
3.4练习
（1）.设定一信号并计算它的fft变换
（2）.模拟信号x(t)=2*sin(4*pi*t)+5*cos(8*pi*t),求N点DFT的幅值谱和相位谱
（3）.比较原序列与经过FFT和IFFT变换后的序列
四,频谱搬移和调制
4.1 频谱搬移和调制的定义
所谓调制，就是利用信号来控制载波的某一参数，使这个参数随该信号而变化。幅度调制是应用较为广泛的一种调制方式，其原理就是使载波的振幅随调制信号的规律而变化，从而达到携带信号进行传输的目的。假设调制信号是普通正弦波：，载波信号为：，那么已调波的数学表达式为：
由上可知，调制后的信号由三个不同频率组成，第一项为未调制的载波；第二项为载波频率与调制频率之和，称为上边频；第三项为载波频率与调制频率之差，称为下边频。显然，调幅过程实际上就是一种频谱搬移过程，经过调制后，调制信号的频谱被搬移到载频附近，成为上边频和下边频。
4.2 频谱搬移和调制在matlab中的实现
编写一个函数lxfft，以实现连续傅立叶变换，其源程序如下：
function y=lxfft(ti,a,p) %ti表示采样间隔，a表示采样后的序列，p为
N=length(a); %频率取样点
for k=0:1:p
y(k+1)=0;
for I=1:1:N
y(k+1)=a(i)*exp(-1*j*(k)*(I-1)*ti)*ti+y(k+1);
end
end
y; %画出变换后的频谱图
I=0:1:p
stem(I,abs(y));
其中，ti表示取样间隔，a表示取样后的序列，p表示频率取样点。先将需实现频谱搬移或调制的信号进行取样，然后经过傅立叶变换，即可看到信号调制前后的频谱关系。
4.3 频谱搬移和调制的应用举例
[例 1] 求正弦信号a=sin（6*t）的频谱。
先将正弦信号进行抽样，然后将抽样的序列进行连续傅立叶变换，即可得出该信号的频谱图。频谱曲线参见图4—2，由图可见，该正弦信号只有单一频率成份。
ti=0.01; %序列的取样间隔
t=0:ti:42; %取样范围
a=sin(6*t); %求出正弦信号的取样序列
p=29; %傅立叶变换频率取样点个数
y=lxfft(ti,a,p); %对其求连续傅立叶变换
[例2] 信号的调制及频谱搬移的实现
已知调制信号为y=sin(t),载波信号为y=5sin(20*t)。
那么已调波的表达式为y=(5+2sin(t))*sin(20*t)，用下面的matlab程序实现，调制波形参见图4-3。
ti=0.01;
t=0:ti:30;
y=(5+2*sin(t))*sin(20*t);
stem(t,y);
频谱搬移的MATLAB实现：
已知调制信号为x1=sin(6*t)，载波信号为x2=sin(20*t)。已调波的表达式为：
y=(5+sin(6*t))*sin(20*t)
先将调制信号和载波信号分别进行连续时间的傅立叶变换，可看到调制前两信号的频谱参见图4-4，调制之后信号频谱曲线参见图4-5。由上图可见信号经过调制后，调制信号的频率分量被搬移到载波频率分量到两边。
ti=0.01;
t=0:ti:100;
p=50;
x1=sin(6*t);
x2=sin(20*t);
subplot(2,1,1);
y1=1xfft(ti,x1,p);
title(’调制信号的频谱’);
subplot(2,1,2);
y2=1xfft(ti,x2,p);
title(‘载波信号的频谱’);
figure(2);
y=(5+sin(6*t))*sin(20*t);
z=1xfft(ti,y,p);
title(‘已调波的频谱’);
4.4练习画出两余弦信号cos(6*t)与cos(50*t)调制前后的频谱图，以及调制后的调幅波。
五,波形分解
5.1波形分解和合成的定义任何信号都是由各种不同频率，幅度和初相的正弦波叠加而成的。对于周期信号，由它的傅氏计数展开式可知，各次谐波的频率是基波频率的整数倍。而非周期信号则包含了从零到无穷大的所有频率成分，每一频率成分的幅度均趋于无限小，但其相对大小则是不同的。利用频率选择滤波器，可以把各个不同频率的谐波分量提取出来，频率选择滤波就是一种专门设计成基本上无失真地通过某些频率，而显著地衰减掉另一些频率的系统。将信号中各个频率分量提取出来的过程就叫做波形分解。而分解出来的各个频率分量的波形相加，就是信号波形的重新合成。
5.2波形分解在matlab中的实现首先，必须先设计滤波器，在matlab中提供了丰富的滤波设计函数，其中cheby2称为切比雪夫II型滤波器设计，格式为：
[b,a]=cheby2(m,Rs,Wn)
其中n为该滤波器的阶数，表示阻带内的纹波（dB），Wn为滤波器的通带。0=<Wn<=1，其中1对应于抽样频率的1/2（Nyquist频率）。当Wn＝【W1 W2】时，cheby2函数产生一个2n阶的数字带通滤波器，通带为W1与W2之间。该函数可以设计模拟或数字的低通、高通、带通、和带阻滤波器。
举例：（1）设抽样频率为1000Hz，要设计一个10阶的低通切比雪夫II型滤波器，截至频率为200Hz，且阻带内的波纹系数为40dB
[b,a]=cheby2(m,Rs,Wn)
（2）设计一个10阶的带通切比雪夫II型滤波器，通带为200-300Hz，抽样频率为1000Hz，阻带纹波系数为40dB
n=10/2
Wn=[200 300]/500;
[b,a]=cheby2(n,40,Wn)
滤波函数filtfilt（b，a，x）利用设计好的滤波器对数据进行滤波。filtfilt（b，a，x）为零位数字滤波格式为：
y＝filtfilt（b，a，x）
通过将输入数据前向和反向处理，以完成零相位数字滤波。它先将数据按顺序滤波，然后将所得结果逆转后反向通过滤波器，这样得到的序列为精确零相位失真，并使滤波器的阶数加倍。
5.3波形分解与合成应用举例对频率为5000Hz的对称周期方波进行分解，并把分解后的波形重新合成，比较分解后合成的信号与源信号的差别。
首先，设计滤波器及滤波函数bandpass（w，x）。参数w表示通带的两个截止频率，x为需要滤波的输入序列。然后利用方波函数square生成一个频率为5000Hz的周期方波，利用设计好的滤波函数bandpass（w，x）对其进行滤波，各次谐波分别用不同颜色来区别，为便于比较，将波形时间轴选为一个方波周期。在波形分解对同时，将分解后各个谐波分量对数据分别保存，然后将这些分量相加，可以比较分解后再合成对波形与原波型对区别。波形分解后各个谐波分量对波型参见图5-1，分解后重新合成对波形参见图5-2
function y=bandpass(w,x)
Wn=w/200000;
[b,a]=cheby2(8/2,40,Wn); %设计一个8阶，阻带内纹波系数为40dB的cheby2II型滤波器
y=filtfilt(b,a,x); %对输入数据进行零相位数字滤波
Fs=400000;%抽样频率
t=0:1/Fs:1;
x=square(2*pi*5000*t); %利用方波产生函数square产生方波
plot(t,x),axis([0.00095 0.00125 –1.5 1.5]),gird on; %画出方波图形曲线
hold on;
for n=1:2:9%利用设计好的滤波器提取各个频率分量
w=[5000*n-3000 5000*n+3000];
y=bandpass(w,x);
switch n
case 1
m=’r’; ％当频率分量不同时，采用不同的线条画图，以示区别
y1=y; ％将滤波后的数据进行保存，以便波形合成时使用
case 3
m=’b’;
y2=y;
case 5
m=’g’;
y3=y;
case 7
m=’y’;
y4=y;
case 9
m=’k’;
y5=y;
end
plot(t,y,m),axis([0.00095 0.00125 –1.5 1.5]);
hold on;
end
Y=y1+y2+y3+y4+y5; %将分解后的谐波分量重新合成
plot(t,y),axis([0.00095 0.00125 –1.5 1.5])； ％画出合成后的波形
5.4练习对周期三角波进行分解，并将分解后的谐波分量重新合成，将分解前的波形和分解后再合成的波形进行比较。
附录表附录1.1
卷积函数convzz
function y=convzz(c,d)
1e=length(c); %获得两个输入序列的长度
1f=length(d);
if 1e<1f %将较长的序列从左到右排列，较短到从
%上到下排列
a=d;
b=c;
else
a=c;
b=d;
end
1a=length(a);
1b=length(b);
if 1a~=prod(size(a)) | 1b~=prod(size(b))
error(A and B must be vectors.’);
end
y=zeros(1,1a+1b-1); %定义一个长度为1a+1b-1到全零序列，用于存
%放卷积结果
for k=1:1b; %第一转折点前点循环运算
i=k; j=1;
while i>0
y(k)=a(i)*b(j)+y(k);
i=i-1; j=j+1;
end
end
for k=(1b+1):1a %第二转折点前的循环运算]
i=k; j=1;
while i>k-1b
y(k)=a(i)*b(j)+y(k);
i=i-1; j=j+1;
end
end
for k=(1a+1):(1a+1b-1) %第二转折点后的循环运算
i=1a; j=k-1a+1;
while i>k-1b
y(k)=a(i)*b(j)+y(k);
i=i-1; j=j+1;
end
end
附录2.1
离散傅立叶变换fftzz及其反变换ifftzz函数
离散傅立叶变换fftzz 函数源程序：
%离散序列的傅立叶变换
function y=fftzz(a,p) %a表示输入的序列，p表示频率取样点个数
N=length(a); %获得a序列的长度
y=zeros(1,p); %定义一个长度为p的全零序列
%外循环
for k=1:1:p
%内循环
for n=1:1:N
y(k)=a(n)*exp(-j*2*pi*(k-1)*(n-1)/p)+y(k); %在某一频率点上对
end %全部序列值运算
if abs(y(k))<10^ (-6) %当y(k)的模小于某一值时，将其置零
y(k)=0; %避免画相位图时出现不确定性
end
end
附录2.2
离散序列傅立叶反变换ifftzz函数源程序
%离散序列傅立叶反变换
function y=ifftzz(b,p)
N=length(b);
Y=zeros(1,p);
for k=1:1:N
for n=1:1:p
y(k)=(1/p)*b(n)*exp(j*2*pi*(k-1)*(n-1)/p)+y(k);
if abs(y(k))<10^(-6)
y(k)=0;
end
end
end
三、实验内容部分实验一 序列的相加和相乘
一、实验内容给出两个序列x (n)和x (n).
x=[0,1,2,3,4,3,2,1,0];n1=[-2:6];
x=[2,2,0,0,0,-2,-2];n=[2:8];
现在求它们的和y及乘积y.
二、实验目的巩固序列的相乘和相加运算。
对MARLAB在同位项相乘和详加中的应用。
三、运行程序参考结果
实验二 序列的移位和周期延拓运算
实验内容已知x(n)=0.8R (n),利用MATLAB生成并图示x(n),x(n-m),x((n))R(n)x((n))表示x(n)以8为周期的延拓和x((n))R(n)，其中N=24,m为一个整常数，0<m<N.
实验目的
1.
2.
三、参考结果图
实验三 信号的采样
实验内容令X(t)=e ,求出并绘制其傅立叶变换X(j)。用三个不同的采样，分别求出并画出离散时间傅立叶变换X(e).三个频率分别为：
(1) f=5kHz; (2) f=2kHz; (3) f=1kHz.
实验目的
1.
2.
三,参考结果图
实验四 基本序列的离散傅立叶变换计算一、实验内容用复正弦序列x (n)=
余弦序列x (n)=
正弦序列x (n)=
分别对N=16和N=8计算序列的N点DFT，并绘制幅频曲线，最后用DFT理解为何两
种N值下的DFT结果差别如此之大。
实验目的
1.
2.
三,参考结果图
实验五 用Z变换求解差分方程一、实验内容令y(n)-2y(n-1)+3y(n-2)=u(n)+ux(n-1)+5u(n-2)-6u(n-3),
初始条件为x(-1)=1,x(-2)=1,y(-1)=-1,y(-2)=1 。
求y(n)。
二、实验目的
1.
2.
三、参考结果图
实验六 快速傅立叶变换
一、实验内容用快速卷积法计算下面两个序列的卷积。并测试直接卷积和快速卷积的时间。
x(n)=
h(n)=
二、实验目的
1.
2.
参考结果图
实验七 用DFT对连续信号作谱分析一、实验内容
用DFT分析的频谱结构。选择不同的截取长度，观察DFT进行频谱分析时存在的截断效应（频谱泄露和谱间干扰）。试用加窗的方法减少谱间干扰。
二、实验目的
1.
2.
三、参考结果图
实验八 用各种窗函数设计FIR数字滤波器实验内容分别用矩形窗和Hamming窗设计线性相位FIR低通滤波器。要求通带截止频率=，单位脉冲响应h(a)的长度N=21。绘制及其幅频响应特性曲线。
实验目的
1.
2.
参考结果图
,
实验九 切比雪夫型低通数字滤波器设计实验内容设计一个切比雪夫型低通数字滤波器设计，指标如下：
通带边界频率：,通带最大衰减：
阻带截止频率： 阻带最小衰减：。
实验目的
1.
2.
参考结果图
实验参考程序一、解：（1）MATLAB程序实现如下：
x1=[0,1,2,3,4,3,2,1,0];ns1=-2; %给定x1及ns1
x2=[2,2,0,0,0,-2,-2];ns2=2; %给定x2及ns2
nf1=ns1+length(x1)-1;nf2=ns2+length(x2)-1;
ny=min(ns1,ns2):max(nf1,nf2); %y(n)的时间变量
xa1=zeros(1,length(ny));xa2=xa1; %延拓序列初始化
xa1(find((ny>=ns1)&(ny<=nf1)==1))=x1; %给xa1赋值x1
xa2(find((ny>=ns2)&(ny<=nf2)==1))=x2; %给xa2赋值x2
ya=xa1+xa2 %序列相加
yp=xa1.*xa2 %序列相乘
subplot(2,2,1),stem(ny,xa1,'.') %绘图
subplot(2,2,2),stem(ny,xa2,'.')
line([ny(1),ny(end)],[0,0]) %画x轴
subplot(2,2,3),stem(ny,ya,'.')
line([ny(1),ny(end)],[0,0])
subplot(2,2,4),stem(ny,yp,'.')
line([ny(1),ny(end)],[0,0])
（2）程序运行结果如图所示：（略）
（3）结论为：
ny=-2 –1 0 1 2 3 4 5 6 7 8
xa1=0 1 2 3 4 3 2 1 0 0 0
xa2=0 0 0 0 2 2 0 0 0 -2 -2
ya= 0 1 2 3 6 5 2 1 0 -2 -2
yp= 0 0 0 0 8 6 0 0 0 0 0
二、解：（1）MATLAB程序实现如下：
clear;close all
N=24;M=8;
m=3;
if (m<1|m>=N-M+1) %检验输入参数n是否合理
fprintf;
end
n=0:N-1;
x1=(0.8).^n; x2=[(n>=0)&(n<M)]; %产生x(n)
xn=x1.*x2;
xm=zeros(1,N); %设定xm的长度
for k=m+1:m+M
xm(k)=xn(k-m);
end
xc=xn(mod(n,8)+1); %产生x(n)的周期延拓
xcm=xn(mod(n-m,8)+1); %产生x(n)移位后的周期延拓
subplot(2,2,1),stem(n,xn,'.')
subplot(2,2,2),stem(n,xm,'.')
subplot(2,2,3),stem(n,xc,'.')
subplot(2,2,4),stem(n,xcm,'.')
三、解：（1）MATLAB程序实现如下：
Dt=0.00005;t=-0.005:Dt:0.005;xa=exp(-1000*abs(t)); %模拟信号
Wmax=2*pi*2000;K=500;k=0:1:K;W=k*Wmax/K;
Xa=xa*exp(-j*t'*W)*Dt;Xa=real(Xa); %连续时间傅立叶变换
W=[-fliplr(W),W(2:501)]; %频率从-Wmax to Wmax
Xa=[fliplr(Xa),Xa(2:501)]; %Xa介于-Wmax和Wmax间
subplot(421);plot(t*1000,xa);
ylabel('xa(t)');
subplot(422);plot(W/(2*pi*1000),Xa*1000);
ylabel('Xa(jw)');
Ts=0.0002;n=-25:1:25;x=exp(-1000*abs(n*Ts)); %离散信号
K=500;k=0:1:K;w=pi*k/K;
X=x*exp(-j*n'*w);X=real(X); %离散时间傅立叶变换
w=[-fliplr(w),w(2:K+1)];
X=[fliplr(X),X(2:K+1)];
subplot(423);stem(n*Ts*1000,x);
ylabel('xl(n)');
subplot(424);plot(w/pi,X);
ylabel('Xl(w)');
Ts=0.0005;n=-10:1:10;x=exp(-1000*abs(n*Ts)); %离散信号
K=500;k=0:1:K;w=pi*k/K;
X=x*exp(-j*n'*w);X=real(X); %离散时间傅立叶变换
w=[-fliplr(w),w(2:K+1)];
X=[fliplr(X),X(2:K+1)];
subplot(425);stem(n*Ts*1000,x);
ylabel('xl(n)');
subplot(426);plot(w/pi,X);
ylabel('Xl(w)');
Ts=0.0010;n=-5:1:5;x=exp(-1000*abs(n*Ts)); %离散信号
K=500;k=0:1:K;w=pi*k/K;
X=x*exp(-j*n'*w);X=real(X); %离散时间傅立叶变换
w=[-fliplr(w),w(2:K+1)];
X=[fliplr(X),X(2:K+1)];
subplot(427);stem(n*Ts*1000,x);
ylabel('xl(n)');
subplot(428);plot(w/pi,X);
ylabel('Xl(w)');
四、解：（1）MATLAB程序实现如下：
clear; close all
N=16; N1=8;
n=0:N-1; k=0:N1-1;
x1n=exp(j*pi*n/8); %产生x1(n)
X1k=fft(x1n,N); %计算N点DFT[x1(n)]
Xk1=fft(x1n,N1); %计算N1点DFT[x1(n)]
x2n=cos(pi*n/8);
X2k=fft(x2n,N); %计算N点DFT[x2(n)]
Xk2=fft(x2n,N1); %计算N1点DFT[x1(n)]
x3n=sin(pi*n/8);
X3k=fft(x3n,N); %计算N点DFT[x3(n)]
Xk3=fft(x3n,N1); %计算N1点DFT[x1(n)]
subplot(2,3,1),
stem(n,abs(X1k),'.'),
subplot(2,3,2),
stem(n,abs(X2k),'.'),
subplot(2,3,3),
stem(n,abs(X3k),'.'),
subplot(2,3,3),
stem(n,abs(X1k),'.'),
subplot(2,3,4),
stem(k,abs(Xk1),'.'),
subplot(2,3,5),
stem(k,abs(Xk2),'.'),
subplot(2,3,6),
stem(k,abs(Xk3),'.')
五、解：（1）MATLAB的程序运行如下：
b=[1 4 5 -6];
a=[1 -2 3];
x0=[1 1];
y0=[-1 1];
xic=filtic(b,a,y0,x0)
bxplus=1;
axplus=[1 -1];
ayplus=conv(a,axplus)
byplus=conv(b,bxplus)+conv(xic,axplus)
[R,P,K]=residuez(byplus,ayplus)
Mp=abs(P)
Ap=angle(P)*180/pi
N=100;
n=0:N-1;
xn=ones(1,N);
yn=filter(b,a,xn,xic);
plot(n,yn)
（2）程序结果如图所示：
（略）
（3）结果为：
xic = 4 2 -6
ayplus = 1 -3 5 -3
byplus = 5 2 -3 0
R = 1.5000 - 4.2426i
1.5000 + 4.2426i
2.0000
P = 1.0000 + 1.4142i
1.0000 - 1.4142i
1.0000
K = 0
Mp = 1.7321
1.7321
1.0000
Ap = 54.7356
-54.7356
0
六、解：（1）MATLAB程序实现如下：
clear;close all
xn=sin(0.4*[1:15]');
hn=0.9.^[1:20]';
M=length(xn);N=length(hn);
nx=1:M;nh=1:N;
%循环卷积等于线性卷积的条件：循环卷积区间长度L>=M+N-1
L=pow2(nextpow2(M+N-1)); %取L为大于等于且最接近(N+M-1)的2的正次幂
tic,%快速卷积计时开始
Xk=fft(hn,L); %L点FFT[x(n)]
Hk=fft(hn,L); %L点FFT[h(n)]
Yk=Xk.*Hk; %频域相乘得Y(k)
yn=ifft(Yk,L); %L点IFFT得到卷积结果y(n)
toc %快速卷积计时结束
subplot(2,2,1),stem(nx,xn,'.');
ylabel('x(n)')
subplot(2,2,2),stem(nh,hn,'.');ylabel('h(n)')
subplot(2,1,2);ny=1:L;
stem(ny,real(yn),'.');ylabel('y(n)')
tic,yn=conv(xn,hn);toc %直接调用函数conv计算卷积与快速卷积比较
实验结果图如下略
（3）结论为：
elapsed_time = 0.0600
elapsed_time = 0
七、解,(1) Matlab程序为：
clear; close all
fs=400; T=1/fs; %采样视频率为400Hz
Tp=0.04;N=Tp*fs; %采样点数N
N1=[N,4*N,8*N]; %设定三种截取长度供调用
st=['|X1(jf)|';'|X4(jf)|';'|X8(jf)|']; %设定三种标注语句供调用
%矩形窗截断
for m=1:3
n=1:N1(m);
%产生采样序列x(n)
xn=cos(200*pi*n*T)+sin(100*pi*n*T)+cos(50*pi*n*T);
Xk=fft(xn,4096); %4096点DFT，用FFT实现
fk=[0:4095]/4096/T;
subplot(3,2,2*m-1)
plot(fk,abs(Xk)/max(abs(Xk)));ylabel(st(m,:))
if m==1 title(' ');end
end
%加hamming窗改善谱间干扰
for m=1:3
n=1:N1(m);
wn=hamming(N1(m)); %调用工具箱函数hamming产生N长hamming窗序列wn
xn=(cos(200*pi*n*T)+sin(100*pi*n*T)+cos(50*pi*n*T)).*wn';
Xk=fft(xn,4096); %4096点DFT，用FFT实现
fk=[0:4095]/4096/T;
subplot(3,2,2*m)
plot(fk,abs(Xk)/max(abs(Xk)));ylabel(st(m,:))
if m==1 title(' ');end
end
八、解,(1) Mtlab的程序为:
clear;close all %理想低通滤波器参数
N=21;wc=pi/4;
n=0:N-1;r=(N-1)/2;
hdn=sin(wc*(n-r))/pi./(n-r); %计算理想低通单位脉冲响应
if rem(N,2)~=0 hdn(r+1)=wc/pi;end %N为奇数时，处理n=r点的0/0型
wn1=boxcar(N); %矩形窗
hn1=hdn.*wn1'; %加窗
%以上两条语句可代以fir函数：hn1=fir1(N-1,wc/pi,boxcar(N));
subplot(2,1,1),
stem(n,hdn,'.'),hw=fft(hdn,512),w=2*[0:511]/512,subplot(2,1,2),plot(w,20*log10(abs(hw))
九、解,(1) Matlab的程序为:
clear;close all
wp=0.2;ws=0.4;Rp=1;Rs=80; %输入指标
[N,wc]=cheb2ord(wp,ws,Rp,Rs) %求滤波器阶次
[B,A]=cheby2(N,Rs,wc) %设计滤波器，得出系数
freqz(B,A) %无左端变量时自动画频率特性图
输出滤波参数如下：
N = 8
Wc = 0.3682
B = Columns 1 through 7
0.0011 0.0008 0.0026 0.0024 0.0032 0.0024 0.0026
Columns 8 through 9
0.0008 0.0011
A = Columns 1 through 7
1.0000 -4.3881 8.9294 -10.8479 8.5205 -4.4046 1.4577
Columns 1 through 7
参考文献
1.MATLAB信号处理详解 陈亚勇 等编著 人民邮电出版社
2.MATLAB 5.X程序设计语言 楼顺天 陈生潭 雷虎民 编著 西安交通大学
出版社
3.信号与系统 ALAN V.OPPENHEIM 刘树棠 译 西安交通大学出版社
4.数字信号处理—基于计算机的方法 Sanjit K.Mitra 清华大学出版社
5.信号与系统解题指南 胡光锐 徐昌庆 等 编著 科学出版社
6.
