实验部分matlab在数字信号处理中的应用.docx
《实验部分matlab在数字信号处理中的应用.docx》由会员分享，可在线阅读，更多相关《实验部分matlab在数字信号处理中的应用.docx（51页珍藏版）》请在冰豆网上搜索。
实验部分matlab在数字信号处理中的应用
Matlab在数字信号处理中的应用（基础）
一、数据类型：
1、整数：
Matlab支持8位，16位，32位和64位的有符号和无符号整数数据类型。
如：
x=int8（50）;%指定x的数据类型为int8.
x=50
2、浮点数：
matlab的默认数据类型是双精度类型（double），为了节省存蓄空间，matlab也支持单精度数据类型的数组。
Realmin（‘single’）
Ans=
1.1755e-038
Realmax（‘double’）
Ans=
2.2251e-308
3、复数：
matlab中虚数单位由i或者j表示。
Z=6+7j
另一种创建复数的方法可以通过complex（）函数，complex（）函数的调用格式：
C=complex（a,b），返回结果c为复数，实部是a，虚部是b。
二、数组的创建
1、一维数组的创建：
创建一维行向量，只需要把所有数组元素用空格或者逗号分隔，并用方括号吧所有数组元素括起来即可。
如用分号，即为列向量。
创建等差的一维数组：
格式Var=start-val:
step:
stop-val。
如果步长是1，可以省略。
2、二维数组的创建；在创建二维数组时，用逗号或者空格区分同一行的不同元素，用分号或者软回车区分不同的行。
三、函数流程控制
1、顺序结构。
2、判断语句（if---elseif---else----end）.
3、循环语句（for----end）
四、作图
1、二维图：
plot（x,y,linespec），linespec参数，用于对图像外观属性的控制，包括线条的形状，颜色和点的形状，颜色。
stem（x,y）;绘制脉冲杆图图形。
Stairs（x,y）;绘制阶梯图图形。
2、图像子窗口：
subplot（m,n,p），将图像分为m╳n个子区域，在第p个区域中绘制图像。
3、坐标轴：
axis（xmin,xmax,ymin,ymax）.指定当前图像中x轴和y轴的范围。
4、图形注释：
1）标题：
title（‘图形名字’）。
2）坐标轴名：
xlabel（‘x轴的名称’），ylabel（‘y轴的名称’）。
特殊符号的输入:
\alpha的输入，则自动转变成，а
实验一、几种典型离散时间序列
Matlab中处理的数组，将下标放在变量后面的小括号内，且约定从1开始递增。
例如：
x=[5,4,3,2,1,0]，表示x
（1）=5,x
（2）=4,x（3）=3,x（4）=2,x（5）=1,x（6）=0。
要表示一个下标不由1开始的数组x（n），一般应采用两个矢量，如：
n=[-3:
5];
x=[1,-1,3,2,0,-2,-1,2,1];这表示一个含有9个点的矢量，n为一组时间矢量，对应x有：
x（-3）=1,x（-2）=-1.。
。
。
。
。
。
连续信号作图使用plot（）函数，绘制线性图。
离散信号作图使用stem（）函数，绘制脉冲杆图。
一些常用的函数：
abs（）:
求绝对值（幅值）。
调用格式：
y=abs（x）。
length（）:
取某一变量的长度（采样点数）。
调用格式：
N=length（n），取n的点数，赋值给N。
real（）：
取一个复数的实部，调用格式：
x=real（h）；取复数h的实部，赋值给变量x。
imag（）:
取复数的虚部，调用格式：
x=imag（h）；取复数h的实部，赋值给变量y
x=sawtooth（t）;类似于sin（t），产生周期为2pi，幅值从-1到+1的锯齿波。
x=sowtooth（t,width）;产生三角波，其中width（0x=square（t）;产生类似于sin（t）,周期为2pi，幅值我1的方波，x=square（t,duty）,产生指定周期的矩形波，其中duty用于指定脉冲宽度与整个周期的比例。
rand（n,m）;产生一组具有n行m列的随机信号。
1、单位冲激序列：
1）利用零序列：
x=zeros（1,N），生成一个1╳N维的零向量。
2）利用逻辑关系表达式产生单位冲激序列：
x=[（n-n0）==0]；只在n=n0的地方产生1.
例：
MATLAB程序如下：
%采样逻辑关系求脉冲序列。
n1=-5;n2=5;n0=0;
n=n1:
n2;
x=[n==n0];
%作图部分
stem（n,x,’filled’）;
axis（[n1,n2,0,1.1*max（x）]）;
title（‘单位脉冲序列’）；
xlabel（‘时间（n）’）;ylabel（‘幅度x（n）’）;
%采样零序列求脉冲序列。
n1=-5;n2=5;k=0;
n=n1:
n2;
nt=length（n）;%求采样点n的个数（长度）。
nk=abs（k-n1）+1;
x=zeros（1,nt）;
x（nk）=1;
%作图同上。
2、单位阶跃序列：
1）利用1序列：
x=ones（1,N），产生一个1╳N维的全1向量。
2）利用逻辑关系表达式产生单位阶跃序列：
x=[（n-n0）>=0]。
Matlab程序：
n=0:
49;
x=ones（1,50）;
closeall;
stem（n,x）;
title（‘单位阶跃信号序列’）；
3、单位矩形序列：
1）x=ones（1,N），
2）利用逻辑关系表达式产生：
x=[（（n-n0>=0）&（n-nf<=0））]。
matlab程序：
N=10;
n=0:
49;
x=sign（sign（N-1-n）+1）;
closeall;%关闭所有打开的图形窗口
stem（n,x）;
注：
sign（x）,符号函数，当x大于0时值为1，当x等于0时值为0，当x小于0时值为-1.
4、正弦序列：
x=a*sin（omega*n+thwlta）；
x=a*sin（2*pi*f0/Fs*n+thelta）；
例：
频率为1.振幅为1的正弦信号，在窗口中显示2个周期的信号波形，并对该信号的一个周期进行32点采样。
获得离散信号。
做出连续信号和离散信号的图形。
MATLAB程序如下：
f=1;Um-1;nt=2;%频率，振幅，周期的个数。
N=32;T=1/f;%采样点数，周期
dt=T/N;%采样时间间隔
n=0:
nt*N-1;
tn=n*dt;
x=Um*sin（2*f*pi*tn）;
%作图部分
subplot（2,1,1）,plot（tn,x）;
axis（[0,nt*T,1.1*min（x）,1.1*max（x）]）;
ylabel（‘连续正弦信号x（t）’）;
subplot（2,1,2）,stem（tn,x）;
axis（[0,nt*T,1.1*min（x）,1.1*max（x）]）;
ylabel（‘离散正弦序列x（n）’）;
5、实指数序列
x（n）=a^n;
例：
编写产生a=1/2和a=2的实指数连续信号和离散信号序列的程序。
MATLAB程序如下：
n1=-10;n2=10;a1=1/2;a2=2;
na1=-10:
0;na2=0:
10;
x1=a1.^na1;
x2=a2.^na2;
%作图部分
subplot（2,2,1）,plot（na1,x1）;%作连续图形
title（‘实指数原信号（a<1）’）;
subplot（2,2,2）,stem（na1,x1,’filled’）;%作离散图形。
title（‘实指数序列（a<1）’）;
subplot（2,2,3）,plot（na2,x2）;%作连续图形
title（‘实指数原信号（a>1）’）;
subplot（2,2,4）,stem（na2,x2,’filled’）;%作离散图形。
title（‘实指数序列（a>1）’）;
6、复指数序列：
x=exp（（sigma+jomega）*n）;
7、矩形波序列：
y=rectpuls（t,width）.该函数产生一个幅度为1宽度为width，且以t=0为对称轴的矩形脉冲信号，width的默认值为1.
y=square（t,DUTY）,产生一个周期为2*pi,幅值为+1（-1）的周期性方波信号。
其中DUTY表示信号的占空比。
默认值为0.5.
例：
矩形脉冲信号的波形图：
2（0<=t<=1）
f（t）=
0（t<1,t>1）
MATLAB程序如下：
t=-0.5:
0.01:
3;
t0=0.5;width=1;
ft=2*rectpuls（t-t0,width）;
plot（t,ft）;
gridon;
axis（[-0.5,3,0.2,2.2]）;
title（‘矩形脉冲信号’）；
例：
产生一个频率为10HZ，占空比为30%的周期方波信号。
MATLAB程序如下：
t=0:
0.001:
3;
y=square（2*pi*10*t,30）;
plot（t,y）;
aixs（0,0.3,-1.2,1.2）;
title（‘周期方波信号’）；
实验二、序列的基本运算
1、序列的加法和乘法
x=x1+x2;x=x1.*x2
例：
已知x1（n）=u（n+2）（-4x2（n）=u（n-4）（-5求：
x（n）=x1（n）+x2（n）
MATLAB程序如下：
n1=-4:
6;n01=-2;
x1=[（n1-n01）>=0];
n2=-5:
8;n02=4;
x2=[[（n2-n02）>=0];
%用0值来扩展它们的序列号，变成相同的起点和终点，原来的值不变。
n=min（[n1,n2]）:
max（[n1,n2]）;
N=length（n）;
y1=zeros（1,N）;y2=zeros（1,N）;
y1（find（（n>=min（n1））&（n<=max（n1））））=x1;
y2（find（（n>=min（n2））&（n<=max（n2））））=x2;
x=y1+y2;
subplot（3,1,1）,stem（n,y1）;
axis（[min（n）,max（n）,1.1*min（x）,1.1*max（x）]）;
subplot（3,1,2）,stem（n,y2）;
axis（[min（n）,max（n）,1.1*min（x）,1.1*max（x）]）;
subplot（3,1,3）,stem（n,x）;
axis（[min（n）,max（n）,1.1*min（x）,1.1*max（x）]）;
注：
序列的乘法与上程序相同。
2、序列的翻转：
翻转运算用fliplr（）函数实现，设序列x（n），样值向量x和位置向量nx表示，则翻转之后的序列y（n）的样值向量y和位置向量ny表示，则
y=fliplr（x）;ny=-fliplr（nx）;
3、序列的移位
设序列x（n），样值向量x和位置向量nx表示，移位n0之后的序列y（n）的样值向量y和位置向量ny表示，则
y=x;
ny=nx+n0;
例：
已知一正弦信号：
x（n）=2sin（2pi*n/10）求其移位信号x（n-2）在-2MATLAB程序如下：
n=-2:
10;n0=2；
x=2*sin（2*pi*n/10）;%建立原信号x（n）
x1=2*sin（2*pi*（n-n0）/10）;%建立x（n-2）信号
subplot（2,1,1）,stem（n,x,’filled’）;
ylabel（‘x（n）’）;
subplot（2,1,2）,stem（n,x1,’filled’）;
4、离散序列卷积：
MATLAB提供一个conv函数
功能：
进行两个序列间的卷积运算。
调用格式：
y=conv（x,h）,用于求取两个有限长度序列x和h的卷积，y的长度等于x和h的长度之和减1。
注意：
conv函数默认两个信号的时间序列从n=0开始的。
例：
已知两个信号序列；
（0（02、如果两个信号不是从n=0开始的，则采用[y,ny]=conv_new（x,nx,h,nh）函数。
其中x是输入序列，nx是它的序列号，h是另一个序列，nh是它的序列号，y是卷积和，ny是它的序列号。
%conv_new.m实现任意位置序列卷积运算，返回值是卷积和的值y和时间向量ny。
function[y,ny]=conv_new（x,nx,h,nh）
n1=nx
（1）+nh
（1）;
n2=nx（length（x））+nh（length（h））;
ny=[n1:
n2];
y=conv（x,h）;
例：
x=[3,11,7,0,-1,4,2];nx=[-3:
3];
h=[2,3,0,-5,2,1,];nh=[-1:
4];
[y,ny]=conv_new（x,nx,h,nh）;
subplot（3,1,1）,stem（nx,x）;
axis（[min（nx）,max（nx）,1.1*min（x）,1.1*max（x）]）;
subplot（3,1,2）,stem（nh,h）;
axis（[min（nh）,max（nh）,1.1*min（h）,1.1*max（h）]）;
subplot（3,1,3）,stem（ny,y）;
axis（[min（ny）,max（ny）,1.1*min（y）,1.1*max（y）]）;
实验三：
离散系统的冲激响应和阶跃响应
1、impz（）：
功能：
求解数字系统的冲激响应。
调用格式：
[h,t]=impz（b,a）;求解数字系统的冲激响应h，取样点数为缺省值。
[h,t]=impz（b,a,n）;求解数字系统的冲激响应h，取样点数为n值。
impz（b,a）;在当前窗口用stem（t,h）函数出图。
2、dstep（）:
功能：
求解数字系统的阶跃响应。
调用格式：
[h,t]=dstep（b,a）;求解数字系统的阶跃响应h，取样点数为缺省值。
[h,t]=dstep（b,a,n）;求解数字系统的阶跃响应h，取样点数为n值。
dstep（b,a）；在当前窗口用stairs（t,h）函数出图。
3、filter子函数
功能：
对数字系统的输入信号进行滤波处理。
因为一个离散系统可以看作是一个滤波器，系统的输出就是输入经过滤波器滤波的结果。
调用格式：
y=filtet（b,a,x），对于由矢量b,a决定的数字系统（b和a分别表示系统函数H（z）对应的分子项和分母项系数构成的数组，而且分母系数要归一化处理。
）当输入信号为x时，对x中的数据进行滤波，结果存于y中，长度取max（na,nb）.
[y,zf]=filter（b,a,x）;除得到结果矢量y外，还得到x的最终状态矢量zf。
y=filter（b,a,x,zi）；可在zi中指定x的初始状态。
4、filtic子函数
功能：
为filter子函数选择初始条件。
调用格式：
zi=filtic（b,a,y,x）；求给定输入x和y时的初始状态。
zi=filtic（b,a,y）；求x=0，给定输入y时的初始状态。
其中，x和y分别是表示过去的输入和输出。
例：
已知一个因果系统的差分方程为6y（n）+2y（n-2）=x（n）+3x（n-1）+3x（n-2）+x（n-3）满足初始条件y（-1）=0,x（-1）=0.求系统的单位冲激响应和阶跃响应。
将上述方程对y（n）项系数进行归一化，得到其系统函数分子和分母系数
a0=1,a1=0,a2=1/3,a3=0
b0=1/6,b1=1/2,b2=1/2,b3=1/6
用impz（）函数的MATLAB程序（取N=32点作图）
a=[1,0,1/3,0];
b=[1/6,1/2,1/2,1/6];
N=32;
n=0:
N-1;
hn=impz（b,a,n）;
gn=dstep（b,a,n）;
subplot（1,2,1）,stem（n,hn,’filled’）;
title（‘系统的单位阶跃响应’）；
ylabel（‘h（n）’）;xlabel（‘n’）;
axis（[0,N-1,-1.1*min（hn）,1.1*max（hn）]）;
subplot（1,2,2）,stem（n,gn,’k’）;
title（‘系统的单位阶跃响应’）；
ylabel（‘g（n）’）;xlabel（‘n’）;
axis（[0,N-1,-1.1*min（gn）,1.1*max（gn）]）;
用filter（）函数的MATLAB程序（取N=32点作图）
a=[1,0,1/3,0];
b=[1/6,1/2,1/2,1/6];
xi=filtic（b,a,0,0）;
N=32;
n=0:
N-1;
x1=[n==0];%单位冲激信号
hn=filter（b,a,x1,xi）;
subplot（1,2,1）,stem（n,hn,’filled’）;
title（‘系统的单位阶跃响应’）；
ylabel（‘h（n）’）;xlabel（‘n’）;
axis（[0,N-1,-1.1*min（hn）,1.1*max（hn）]）;
x2=[n>=0];%单位阶跃信号
gn=filter（b,a,x2,xi）;
subplot（1,2,2）,stem（n,gn,’k’）;
title（‘系统的单位阶跃响应’）；
ylabel（‘g（n）’）;xlabel（‘n’）;
axis（[0,N-1,-1.1*min（gn）,1.1*max（gn）]）;
注：
1、hold（）控制当前图形是否刷新的双向切换开关。
holdon使当前轴及图形保持而不被刷新，准备接受此后将绘制的新曲线。
holdoff使当前轴及图形不再具备刷新的性质。
3、pause（）暂停执行文件，等待用户按任意键继续，pause（n）在继续执行之前，暂停n秒。
实验四：
离散LSI系统的时域响应
对于离散LSI系统的响应，MATLAB为我们提供了多种求解方法：
（1）用conv子函数进行卷积积分，求任意输入的系统零状态响应。
（2）用dlsim子函数求任意输入的系统零状态响应。
（3）用filter和filtic子函数求任意输入的系统完全响应。
1、dlsim子函数
功能：
求解离散系统的响应。
调用格式：
y=dlsim（b,a,x）,求输入信号为x时系统的响应。
其中，b和a分别表示系统函数H（z）中，由对应的分子项和分母项系数构成的数组，而且分母系数要归一化处理。
例：
（书本P96页15题）已知一个用以下差分方程表示的线性移不变因果系统为：
当激励
时，求系统的响应。
2、用filter和filtic子函数求LSI系统对任意输入的响应
Filter子函数
功能：
对数字系统的输入信号进行滤波处理。
调用格式：
y=filtet（b,a,x），对于由b和a决定的数字系统（b和a分别表示系统函数H（z）中，由对应的分子项和分母项系数构成的数组，而且分母系数要归一化处理。
）当输入信号为x时，对x中的数据进行滤波，结果存于y中，长度取max（na,nb）.
[y,zf]=filter（b,a,x）;除得到结果矢量y外，还得到x的最终状态矢量zf。
y=filter（b,a,x,zi）；可在zi中指定x的初始状态。
Filtic子函数
功能：
为filter子函数选择初始条件。
调用格式：
zi=filtic（b,a,y,x）；求给定输入x和y时的初始状态。
zi=filtic（b,a,y）；求x=0，给定输入y时的初始状态。
其中，x和y分别是表示过去的输入和输出。
例：
已知一个因果系统的差分方程为：
6y（n）-2y（n-4）=x（n）-3x（n-2）+3x（n-4）-x（n-6）
满足初始条件y（-1）=0,x（-1）=0，求系统的单位冲激响应和单位阶跃响应，时间轴上N取32点作图。
例：
已知一个IIR数字低通滤波器的系统函数为：
=
输入一个矩形信号序列x=square（n/5）（-2）
作业：
书本p96页16题。
作业：
1、已知离散线性时不变的系统函数，请分别用impz和dstep子函数，filter和filtic子函数两种方法求解系统的冲激响应和阶跃响应。
（1）
（2）
2、一个LSI系统的系统函数表示式为：
满足初始条件y（-1）=5,y（-2）=5，试用filtet和filtic子函数求此系统的输入序列x（n）为下列信号时的零输入，零状态以及完全响应。
（1）x（n）=
（2）x（n）=
实验五：
z变换及其应用
1、ztrans子函数
功能：
返回无限长序列函数x（n）的z变换。
调用格式：
z=ztrans（x）；求无限长序列函数x（n）的z变换X（z），返回z变换的表达式。
2、iztrans子函数
功能：
求函数X（z）的z反变换x（n）.
调用格式：
x=iztrans（X（z））；求函数X（z）的z反变换x（n），返回z反变换的表达式。
3、syms子函数
功能：
定义多个符号对象。
调用格式：
symsabwo;把字符abwo定义为基本的符号对象。
4、residuez子函数
功能：
有理多项式的部分分式展开。
调用格式：
[rpc]=residuez（b,a）；把b（z）/a（z）展开成部分分式的形式。
[b,a]=residuez（rpc）；根据部分分式的rpc数组，返回有理多项式。
其中：
b,a为按降幂排列的多项式的分子和分母的系数组；r为余数数组；p为极点数组；c为无穷限多项式系数数组。
有理多项式如下：
X（z）=
注意：
利用ztrans（）子函数时，它只给出了z变换的表达式，而没有给出收敛域。
另外，由于这一功能还不尽完美，因而有的序列的z变换还求不出来，z的反变换也存在同样的问题。
例：
用部分分式法求解函数
的z反变换，写出h（n）的表达式，并用图形与impz求得的结果相比较。
MATLAB程序：
%求z的反变换
b=[0,1,0];a=[1,-12,36];
[r,p,c]=residuez（b,a）
%由此可知，这个多项式含有重极点，多项式分解后表示为：
H（z）=-0.1667/（1-6z-1）+0.1667/（1-6z-1）2
=-0.1667/（1-6z-1）+0.1667z/6*6z-1/（1-6z-1）2
根据时域位移性质，可写出z反变换公式
h（n）=-0.1667（6）nu（n）+0.1667/6*（n+1）6n+1u（n+1）.
%作图
N=8;n=0:
N-1;
h=r
（1）*p
（1）.^n.*[n>=0]+r
（2）.*（n+1）.*p
（2）.^n.*[n-1>=0];
subplot（1,2,2）,stem（n,h）;
title（‘用部分分式法求反变换h（n）’）;
%用冲击响应法求h（n）
h1=impz（b,a,N）;
subplot（1,2,2）,stem（n,h1）;
title（‘用impz（）求反变换h（n）’）;
升级会员 