MATLAB实验二修改.docx

MATLAB实验二修改.docx

《MATLAB实验二修改.docx》由会员分享，可在线阅读，更多相关《MATLAB实验二修改.docx（42页珍藏版）》请在冰豆网上搜索。

MATLAB实验二修改

实验二信号的表示及其基本运算

一、实验目的

1、掌握连续信号及其MATLAB实现方法；

2、掌握离散信号及其MATLAB实现方法

3、掌握离散信号的基本运算方法，以及MATLAB实现

4熟悉应用MATLAB实现求解系统响应的方法

4、了解离散傅里叶变换的MATLAB实现

5、了解IIR数字滤波器设计

6、了解FIR数字滤波器设计1

二、实验设备

计算机，Matlab软件

三、实验内容

（一）、连续信号及其MATLAB实现

1、单位冲激信号

例1.1：

单位冲击信号的MATLAB实现程序如下：

t1=-4;

t2=4;

t0=0;

dt=0.01;

t=t1:

dt:

t2;

n=length（t）;

x=zeros（1,n）;

x（1,（-t0-t1）/dt+1）=1/dt;

stairs（t,x）;

axis（[t1,t2,0,1.2/dt]）;

2、任意函数

例1.2：

用MATLAB画出如下表达式的脉冲序列

3单位阶跃函数

例1.3：

用MATLAB实现单位阶跃函数

clearall;

t=-0.5:

0.001:

1;

t0=0;

u=stepfun（t,t0）;

plot（t,u）

axis（[-0.51-0.21.2]）

4斜坡函数

例1.4：

用MATLAB实现g（t）=3（t-1）

clearall;

t=0:

0.01:

3;

B=3;

t0=1;

u=stepfun（t,t0）;

n=length（t）;

fori=1:

n

u（i）=B*u（i）*（t（i）-t0）;

end

plot（t,u）

axis（[-0.23.1-0.26.2]）

5抽样信号抽样信号Sa（t）=sin（t）/t在MATLAB中用sinc函数表示。

定义为

t=-3*pi:

pi/100:

3*pi;

ft=sinc（t/pi）;

plot（t,ft）;

gridon;

axis（[-10,10,-0.5,1.2]）;%定义画图范围，横轴，纵轴

title（'抽样信号'）%定义图的标题名字

6指数函数

例1.5：

用MATLAB实现

7正弦函数

例1.6：

用MATLAB实现正弦函数f（t）=3cos（10πt+1）

8虚指数信号

例虚指数信号调用格式是f=exp（（j*w）*t）

t=0:

0.01:

15;

w=pi/4;

X=exp（j*w*t）;

Xr=real（X）;%取实部

Xi=imag（X）;%取虚部

Xa=abs（X）;%取模

Xn=angle（X）;%取相位

subplot（2,2,1）,plot（t,Xr）,axis（[0,15,-（max（Xa）+0.5）,max（Xa）+0.5]）,

title（'实部'）;

subplot（2,2,3）,plot（t,Xi）,axis（[0,15,-（max（Xa）+0.5）,max（Xa）+0.5]）,

title（'虚部'）;

subplot（2,2,2）,plot（t,Xa）,axis（[0,15,0,max（Xa）+1]）,title（'模'）;

subplot（2,2,4）,plot（t,Xn）,axis（[0,15,-（max（Xn）+1）,max（Xn）+1]）,title（'相角'）;

%subplot（m,n,i）命令是建立m行n列画图窗口，并指定画图位置i

9复指数信号

例复指数信号调用格式是f=exp（（a+j*b）*t）

t=0:

0.01:

3;

a=-1;b=10;

f=exp（（a+j*b）*t）;

subplot（2,2,1）,plot（t,real（f））,title（'实部'）

subplot（2,2,3）,plot（t,imag（f））,title（'虚部'）

subplot（2,2,2）,plot（t,abs（f））,title（'模'）

subplot（2,2,4）,plot（t,angle（f））,title（'相角'）

（二）、离散信号及其MATLAB实现

1、单位冲激序列

例2.1：

用MATLAB产生64点的单位冲激序列

clearall;

N=64;

x=zeros（1,N）;

x

（1）=1;

xn=0:

N-1;

stem（xn,x）

axis（[-16501.1]）

2、任意序列

例2.2：

用MATLAB画出如下表达式的脉冲序列

3、单位阶跃序列

例2.3：

用MATLAB实现单位阶跃函数

4、斜坡序列

例2.4：

用MATLAB实现g（n）=3（n-4）点数为32的斜坡序列

clearall;

N=32;

k=4

B=3;

t0=1;

x=[zeros（1,k）ones（1,N-k）];

fori=1:

N

x（i）=B*x（i）*（i-k）;

end

xn=0:

N-1;

stem（xn,x）

axis（[-132090]）

5、正弦序列

例2.5：

用MATLAB实现幅度A=3，频率f=100，初始相位Φ=1.2，点数为32的正弦信号

6、实指数序列

例2.6：

用MATLAB实现

，点数为32的实指数序列

clearall;

N=32;

A=3;

a=0.7;

xn=0:

N-1;

x=A*a.^xn;

stem（xn,x）

7、复指数序列

例2.7：

用MATLAB实现幅度A=3，a=0.7，角频率ω=314，点数为32的实指数序列

clearall;

N=32;

A=3;

a=0.7;

w=314;

xn=0:

N-1;

x=A*exp（（a+j*w）*xn）;

stem（xn,x）

8、随机序列

利用MATLAB产生两种随机信号：

rand（1,N）在区间上产生N点均匀分布的随机序列

randn（1,N）产生均值为0，方差为1的高斯随机序列，即白噪声序列

例2.8：

用MATLAB产生点数为32的均匀分布的随机序列与高斯随机序列

clearall;

N=32;

x_rand=rand（1,N）;

x_randn=randn（1,N）;

xn=0:

N-1;

figure

（1）;

stem（xn,x_rand）

figure

（2）;

stem（xn,x_randn）

（三）、离散信号的基本运算

1、信号的延迟

给定离散信号x（n），若信号y（n）定义为：

y（n）=x（n-k），那么y（n）是信号x（n）在时间轴上右移k个抽样周期得到的新序列。

例3.1：

正弦序列y（n）=sin（100n）右移3个抽样周期后所得的序列，MATLAB程序如下：

clearall;

N=32;

w=100;

k=3;

x1=zeros（1,k）;

xn=0:

N-1;

x2=sin（100*xn）;

figure

（1）

stem（xn,x2）

x=[x1x2];

axis（[-1N-1.11.1]）

N=N+k;

xn=0:

N-1;

figure

（2）

stem（xn,x）

axis（[-1N-1.11.1]）

利用for循环语句实现周期延迟.

2、信号相加

若信号

，值得注意的是当序列

和

的长度不相等或者位置不对应时，首先应该使两者的位置对齐，然后通过zeros函数左右补零使其长度相等后再相加

例3.2：

用MATLAB实现两序列相加

clearall;

n1=0:

3

x1=[20.50.91];

figure

（1）

stem（n1,x1）

axis（[-1802.1]）

n2=0:

7

x2=[00.10.20.30.40.50.60.7];

figure

（2）

stem（n2,x2）

axis（[-1800.8]）

n=0:

7;

x1=[x1zeros（1,8-length（n1））];

x2=[zeros（1,8-length（n2））,x2];

x=x1+x2;

figure（3）

stem（n,x）

axis（[-1802.1]）

已知f1（t）=sinwt,f2（t）=sin8wt,w=2pi,求f1（t）+f2（t）和f1（t）f2（t）的波形图

3、信号相乘

信号序列

和

相乘所得信号

的表达式为：

这是样本与样本之间的点乘运算，在MATLAB中可采用“.*”来实现，但是在信号序列相乘之前，应对其做与相加运算一样的操作。

例3.3：

用MATLAB实现上例中两序列相乘

clearall;

n1=0:

3

x1=[20.50.91];

figure

（1）

stem（n1,x1）

axis（[-1802.1]）

n2=0:

7

x2=[00.10.20.30.40.50.60.7];

figure

（2）

stem（n2,x2）

axis（[-1800.8]）

n=0:

7;

x1=[x1zeros（1,8-length（n1））];

x2=[zeros（1,8-length（n2））,x2];

x=x1.*x2;

figure（3）

stem（n,x）

axis（[-1800.35]）

4、信号翻转

信号翻转的表达式为：

y（n）=x（-n），在MATLAB中可以用fliplr函数实现此操作

例3.4：

用MATLAB实现“信号相加”中的

序列翻转

clearall;

n=0:

3

x1=[20.50.91];

x=fliplr（x1）;

stem（n,x）

axis（[-1402.1]）

5、信号和

对于N点信号

，其和的定义为：

例3.5：

用MATLAB实现“信号相加”中的

序列和

clearall;

n=0:

3

x1=[20.50.91];

x=sum（x1）

6、信号积

对于N点信号

，其积的定义为：

例3.5：

用MATLAB实现“信号相加”中的

序列积

clearall;

n=0:

3

x1=[20.50.91];

x=prod（x1）

7卷积

1、完成

与

两函数的卷积运算

其中：

在一个图形窗口中，画出

、

以及卷积结果。

要求每个坐标系有标题、坐标轴名称。

p=0.1;

t=0:

p:

10;

f1=exp（-2*t）.*u（t）;

f2=u（t）-u（t-4）;

f=conv（f1,f2）;

subplot（1,3,1）;

plot（t,f1,'r'）;

title（'f1（t）=e^-2*t*u（t）'）;

xlabel（'t（sec）'）;

ylabel（'f1（t）'）;

subplot（1,3,2）;

plot（t,f2,'g'）;

title（'f2（t）=u（t）-u（t-4）'）;

xlabel（'t（sec）'）;

ylabel（'f2（t）'）;

subplot（1,3,3）;

plot（f）;

title（'f（t）=f1（t）*f2（t）'）;

xlabel（'t（sec）'）;

ylabel（'f（t）'）;

四MATLAB在信号与系统中的应用

该实验用MATLAB中库函数，如tf2zp（b,a），ss2zp（A,B,C,D），zplane（z,p），freqz（b,a）等。

例如：

1．传递函数为

，求其零极点图。

程序如下：

num=[10.52]；分子系数，按降幂顺序排列

den=[10.41]；分母系数，按降幂顺序排列

[z,p]=tf2zp（num,den）；用tf2zp函数求出其零点z和极点p

zplane（z,p）作出零极点图

2．若给出的是滤波器的输入与输出的状态方程，如：

，

，求其零极点图。

程序如下：

A=[1，0；1，-3]；

B=[1；0]；

C=[-

，1]；

D=0；

[z,p]=ss2zp（A,B,C,D）；求出其零极点z,p

zplane（z,p）

在连续时间系统中，当极点在虚轴的右半平面时，系统不稳定，在虚轴上，

单阶极点系统稳定；若零点均处于左半平面内，则系统为最小相位系统。

在离散系统中，极点在单位圆外系统不稳定，在单位圆上，单阶极点系统稳定；零点在单位圆内，系统为最小相位系统。

对一滤波器，我们不仅要知道它的零点和极点，还要了解它的频率特性，本实验可求其频率特性。

对模拟滤波器，可用freqs函数求得其频率特性，对数字滤波器，则用freqz函数求得。

3已知模拟滤波器的传递函数为

，求其频率特性。

程序如下：

num=[0.20.31];

den=[10.41];

w=logspace（-1,1）;频率范围

freqs（num,den,w）

例4．数字滤波器

，取样点数为128点，求其频率特性。

程序如下：

num=[0.20.31];

den=[10.41];

freqz（num,den,128）

实验内容:

1．已知下列传递函数H（s）或H（z），求其零极点，并画出零极点图。

a.

b.

c.

d.

2．已知下列H（s）或H（z），求其频响。

a.       

b.       

c.       

d.       

1：

若某连续系统的输入为e（t），输出为r（t），系统的微分方程为：

①求该系统的单位冲激响应h（t）及其单位阶跃响应g（t）。

②若

求出系统的零状态响应y（t）

分析：

   ①求冲激响应及阶跃响应的MATLAB程序：

a=[1 5 6];b=[3 2];

subplot（2,1,1）,impulse（b,a,4）

subplot（2,1,2）,step（b,a,4）

运行结果如右：

   ②求零状态响应的MATLAB程序：

a=[1 5 6];b=[3 2];

p1=0.01;           %定义取样时间间隔为0.01

t1=0:

p1:

5;          %定义时间范围

x1=exp（-2*t1）;      %定义输入信号

lsim（b,a,x1,t1）,      %对取样间隔为0.01时系统响应进行仿真

holdon;            %保持图形窗口以便能在同一窗口中绘制多条曲线

p2=0.5;            %定义取样间隔为0.5

t2=0:

p2:

5;          %定义时间范围

x2=exp（-2*t2）;      %定义输入信号

lsim（b,a,x2,t2）,holdoff  %对取样间隔为0.5时系统响应进行仿真并解除保持

运行结果如下：

2 已知一个过阻尼二阶系统的状态方程和输出方程分别为：

  ，   r（t）=[0 1]X（t） 。

若系统初始状态为X（0）=[4 -5]T ，求系统在

作用下的全响应。

求全响应程序如下：

A=[0 1;-2 -3];B=[0 2]';C=[0 1];D=[0];

X0=[4 -5]';                      %定义系统初始状态

t=0:

0.01:

10;                      

E=[3*exp（-4*t）.*ones（size（t））]'; %定义系统激励信号

[r,x]=lsim（A,B,C,D,E,t,X0）;  %求出系统全响应的数值解

plot（t,r）                          %绘制系统全响应波形

运行结果如右。

3已知描述离散系统的差分方程为：

，且已知系统输入序列为

，

①求出系统的单位函数响应h（k）在-3~10离散时间范围内响应波形。

②求出系统零状态响应在0~15区间上的样值；并画出输入序列的时域波形以及系统零状态响应的波形

  分析：

①求系统的单位函数响应的MATLAB程序：

     a=[1,-0.25,0.5];     b=[1,1,0];         

     impz（b,a,-3:

10）,title（'单位响应'）          %绘出单位函数响应在-3~10区间上的波形

     运行结果如图a。

②求零状态响应的MATLAB程序：

a=[1,-0.25,0.5];b=[1,1,0]

k=0:

15;             %定义输入序列取值范围

x=（1/2）.^k;          %定义输入序列表达式

y=filter（b,a,x）        %求解零状态响应样值

subplot（2,1,1）,stem（k,x）   %绘制输入序列的波形   

title（'输入序列'）

subplot（2,1,2）,stem（k,y）%绘制零状态响应的波形

title（'输出序列'）

运行结果如下：

y=

 Columns1through10

   1.0000   1.7500   0.6875  -0.3281  -0.23830.1982   0.2156  -0.0218  -0.1015  -0.0086

 Columns11through16

0.0515   0.0187  -0.0204  -0.0141   0.00690.0088

图a.①运行结果图b.②运行结果

实验内容

1.已知描述系统的微分方程和激励信号e（t）分别如下，并用MATLAB绘出系统单位冲激响应和系统零状态响应的波形.

①

；

②

；

③

；

④如下图所示的电路中，已知

，

，且两电感上初始电流分别为

，如果以电阻

上电压

作为系统输出，请求出系统在激励

（v）作用下的全响应。

2.请用MATLAB分别求出下列差分方程所描述的离散系统，在0~20时间范围内的单位函数响应、阶跃响应和系统零状态响应的数值解，并绘出其波形。

另外，请将理论值与MATLAB仿真结果在对应点上的值作比较，并说出两者的区别和产生误差的原因。

①

；

②

；

③

；

④一带通滤波器可由下列差分方程描述：

， 其中

为系统输入，

为系统输出。

请求出当激励为

（选取适当的n值）时滤波器的稳态输出。

（四）、傅里叶变换的MATLAB实现

①在MATLAB中实现傅里叶变换的函数为：

        F=fourier（f）    对f（t）进行傅里叶变换，其结果为F（w）

        F=fourier（f,v）  对f（t）进行傅里叶变换，其结果为F（v）

        F=fourier（f,u,v） 对f（u）进行傅里叶变换，其结果为F（v）

②傅里叶反变换

        f=ifourier（F）    对F（w）进行傅里叶反变换，其结果为f（x）

        f=ifourier（F,U）  对F（w）进行傅里叶反变换，其结果为f（u）

        f=ifourier（F,v,u） 对F（v）进行傅里叶反变换，其结果为f（u）

例①求门函数

的傅里叶变换，并画出幅度频谱图

MATLAB程序如下：

symstw                              %定义两个符号变量t,w

Gt=sym（'Heaviside（t+1）-Heaviside（t-1）'）;    %产生门宽为2的门函数

Fw=fourier（Gt,t,w）;                     %对门函数作傅氏变换求F（jw）

FFw=maple（'convert',Fw,'piecewise'）;       %数据类型转换，转为分段函数，此处可以去掉

FFP=abs（FFw）;                        %求振幅频谱|F（jw）|

ezplot（FFP,[-10*pi10*pi]）;grid;           %绘制函数图形，并加网格

axis（[-10*pi10*pi02.2]）                %限定坐标轴范围

运行结果：

Fw=exp（i*w）*（pi*Dirac（w）-i/w）-exp（-i*w）*（pi*Dirac（w）-i/w）

%Dirac（w）为δ（ω），即傅立叶变换结果中含有奇异函数，故绘图前需作函数类型转换

FFw=-i*exp（i*w）/w+i*exp（-i*w）/w%FFw为复数

FFP=abs（-i*exp（i*w）/w+i*exp（-i*w）/w）%求FFw的模值

例②求函数

的傅里叶反变换f（t）

   MATLAB程序如下：

symstw                            %定义两个符号变量t,w

Fw=sym（'1/（1+w^2）'）;                %定义频谱函数F（jw）

ft=ifourier（Fw,w,t）;                  %对频谱函数F（jw）进行傅氏反变换

运行结果：

ft=

1/2*exp（-t）*Heaviside（t）+1/2*exp（t）*Heaviside（-t）

2、傅里叶变换的数值计算实现法

严格说来，如果不使用symbolic工具箱，是不能分析连续时间信号的。

采用数值计算方法实现连续时间信号的傅里叶变换，实质上只是借助于MATLAB的强大数值计算功能，特别是其强大的矩阵运算能力而进行的一种近似计算。

傅里叶变换的数值计算实现法的原理如下：

对于连续时间信号f（t），其傅里叶变换为：

其中τ为取样间隔，如果f（t）是时限信号，或者当|t|大于某个给定值时，f（t）的值已经衰减得很厉害，可以近似地看成是时限信号，则上式中的n取值就是有限的，假定为N，有：

若对频率变量ω进行取样，得：

通常取：

，其中

是要取的频率范围，或信号的频带宽度。

采用MATLAB实现上式时，其要点是要生成f（t）的N个样本值

的向量，以及向量

，两向量的内积（即两矩阵的乘积），结果即完成上式的傅里叶变换的数值计算。

注意：

时间取样间隔τ的确定，其依据是τ必须小于奈奎斯特（Nyquist）取样间隔。

如果f（t）不是严格的带限信号，则可以根据实际计算的精度要求来确定一个适当的频率

为信号的带宽。

例③用数值计算法实现上面门函数

的傅里叶变换，并画出幅度频谱图.

   分析：

该信号的频谱为

，其第一个过零点频率为π，一般将此频率认为是信号的带宽。

但考虑到

的形状（为抽样函数），假如将精度提高到该值的50倍，即取

，则据此确定的Nyquist取样间隔为：

。

MATLAB程序如下：

R=0.02;