数字信号处理实验讲义-长师

前言

数字信号处理的理论和技术一出现就受到人们的极大关注，发展非常迅速。20世纪60年代以来，随着信息科学和计算技术的迅速发展，数字信号处理的理论与应用飞跃发展，它的重要性日益在通信、信息技术、图像处理、遥感、声呐、雷达、生物医学、地震、语音处理等领域的应用中表现出来，并发展成为一门极其重要的学科。

数字信号处理是一门理论与实践紧密联系的课程，适当的上机实验有助于深入理解和巩固难基本理论知识，锻炼初学者用计算机和MA TLAB语言及其工具箱函数解决数字处理算法的仿真和滤波器设计问题的能力。

由于信号处理涉及大量的运算，可以说离开了计算机及相应的软件，就不可能解决任何稍微复杂的实际应用问题。Matlab是1984年美国MathWorks公司的产品，MATLAB语言具备高效、可视化及推理能力强等特点，它的推出得到了各个领域专家学者的广泛关注，其强大的扩展功能为各个领域的应用提供了基础，是目前工程界流行最广的科学计算语言。早在20世纪90年代中期，MATLAB就己成为国际公认的信号处理的标准软件和开发平台。从1996年后，美国新出版的信号处理教材就没有一本是不用MA TLAB的。

本学期实验结合数字信号处理的基本理论和基本内容，设计了五个上机实验，每个实验对应一至两个主题内容。包括常见离散信号的MA TLAB产生和图形显示、离散时间信号的Z变换、快速傅立叶变换FFT 及其应用、数字滤波器的MATLAB实现、用窗函数法设计FIR数字滤波器。

每个实验中均给出了实验方法和步骤，及部分的MATLAB程序，此外，还有习题和上机实验。通过实验可以使学生掌握数字信号处理的基本原理和方法。用科学计算语言MATLAB实现数字信号处理的方法和实践，通过实验用所学理论来分析解释程序的运行结果，进一步验证、理解和巩固学到的理论知识，从而达到掌握数字信号处理的基本原理和方法的目的。

实验一 常见离散信号的MATLAB 产生和图形显示

一、实验目的

１、学会用MA TLAB 在时域中产生一些基本的离散时间信号。 ２、了解信号的各种运算。

二、实验原理

（一）序列的产生

由于MA TLAB 数值计算的特点，用它来分析离散时间信号与系统是很方便的。离散信号是数字信号处理的最基础的内容，由于内存有限，MATLAB 无法表示无限序列。

在MATLAB 中，可以用一个列向量来表示一个有限长度的序列，但是这种表示方法没有包含采样位置的信息，要完全表示x(n)，要用x 和n 两个向量，例如

x(n)={2,1,0,2,3,-1,2,3}，在MATLAB 中表示为n=[-4,-3,-2,-1,0,1,2,3];x=[2,1,0,2,3,-1,2,3];当序列从n=0开始，则不需要采样位置信息，这时可以只用x 来表示

1．正弦序列

离散正弦序列的MATLAB 表示与连续信号类似，只不过是用stem 函数而不是用plot 函数来画出序列的波形。下面就是正弦序列k ??

?

??6sin π的MATLAB 源程序。程序运行结果如图1.19所示。

%正弦序列实现程序 k=0:39;

fk=sin(pi/6*k); stem(k,fk)

图1.19正弦序列波形

2．指数序列

离散指数序列的一般形式为k

ca ，可用MATLAB 中的数组幂运算（即点幂运算）c*k a .^来实现。下面为用MA TLAB 编写绘制离散时间实指数序列波形的函数。

function dszsu(c,a,k1,k2) %c ：指数序列的幅度 %a ：指数序列的底数

%k1：绘制序列的起始序号 %k2：绘制序列的终止序号 k=k1:k2; x=c*(a.^k); stem(k,x,'filled') hold on

plot([k1,k2],[0,0]) hold off

利用上述函数，实现实指数波形MA TLAB 程序如下（其中a 值分别为4

3454345--，，，）。 %离散时间实指数序列实现程序

subplot 221;

dszsu(1,5/4,0,20); xlabel('k'); title('f1[k]'); subplot 222

dszsu(1,3/4,0,20); xlabel('k'); title('f2[k]'); subplot 223;

dszsu(1,-5/4,0,20); xlabel('k'); title('f3[k]'); subplot 224;

dszsu(1,-3/4,0,20); xlabel('k');

title('f4[k]'); 图1.20不同底数的实指数序列

程序运行结果如图1.20所示。如图可知，对于离散时间实指数序列k

ca ，当a 的绝对值大于1时，序列为随时间发散的序列，当a 的绝对值小于1时，序列为随时间收敛的序列。同时可见，当a 的值小于零时，

其波形在增长或衰减的同时，还交替地改变序列值的符号。

对于离散时间虚指数序列，可用通过调用下列绘制虚指数序列时域波形的MA TLAB 函数。 function []=dxzsu(n1,n2,w)

%n1:绘制波形的虚指数序列的起始时间序号 %n2：绘制波形的虚指数序列的终止时间序号 %w ：虚指数序列的角频率 k=n1:n2;

f=exp(i*w*k); Xr=real(f) Xi=imag(f) Xa=abs(f) Xn=angle(f)

subplot(2,2,1),stem(k,Xr,'filled'),title('实部'); subplot(2,2,3),stem(k,Xi,'filled'),title('虚部'); subplot(2,2,2),stem(k,Xa,'filled'),title('模'); subplot(2,2,4),stem(k,Xn,'filled'),title('相角');

利用上述函数，实现虚指数波形MA TLAB 程序如下（其中虚指数分别为k

j k j ，e

e 24

π

）

%离散时间虚指数实现程序 figure(1);

dxzsu(0,20,pi/4); figure(2); dxzsu(0,20,2);

程序运行结果如图1.21（a ）、（b ）所示。由图可见，只有当虚指数序列的角频率满足ωπ2为有理数

时，信号的实部和虚部和相角都为周期序列，否则为非周期序列。

（a ）4

πk j e

波形(b)k

j e

2波形

图1.21虚指数序列波形

对于复指数序列，其一般形式为

[]k j k e r k f ω=

可以通过调用下面绘制复指数序列时域波形的MA TLAB 函数。

function dfzsu(n1,n2,r,w)

%n1:绘制波形的虚指数序列的起始时间序号 %n2：绘制波形的虚指数序列的终止时间序号 %w ：虚指数序列的角频率 %r:指数序列的底数 k=n1:n2;

f=(r*exp(i*w)).^k; Xr=real(f); Xi=imag(f); Xa=abs(f); Xn=angle(f);

subplot(2,2,1),stem(k,Xr,'filled'),title('实部'); subplot(2,2,3),stem(k,Xi,'filled'),title('虚部'); subplot(2,2,2),stem(k,Xa,'filled'),title('模'); subplot(2,2,4),stem(k,Xn,'filled'),title('相角');

利用上述函数，实现复指数序列波形MA TLAB 程序如下。 %复指数序列实现程序(r>1) figure(1);

dfzsu(0,20,1.2,pi/4);

%复指数序列实现程序(0

dfzsu(0,20,0.8,pi/4);

%复指数序列实现程序(r=1) figure(3);

dfzsu(0,20,1,pi/4);

其运行结果如图1.22（a ）、（b ）、（c ）所示。如图可见，当r>1时，复指数序列的实部和虚部分别为幅度按指数增长的正弦序列;当0

3．单位抽样序列

可以通过借助MA TLAB 中的零矩阵函数zeros 表示。全零矩阵zeros(1,N)产生一个由N 个零组成的列向

量，对于有限区间的[]kδ可以通过以下MATLAB程序表示%单位抽样序列实现程序

k=-30:30;

delta=[zeros(1,30),1,zeros(1,30)];

stem(k,delta)

程序运行结果如图1.23所示。

(a)r>1(b)0

(c)r=1

图1.22复指数序列波形

图1.23单位抽样序列波形图1.24单位阶跃序列波形

4．单位阶跃序列

可以通过借助MATLAB 中的单位矩阵函数ones 表示。单位矩阵ones(1,N)产生一个由N 个1组成的列向量，对于有限区间的[]k u 可以通过以下MATLAB 程序表示

%单位阶跃序列实现程序 k=-30:30;

uk=[zeros(1,30),ones(1,31)]; stem(k,uk)

程序运行结果如图1.24所示。

例1.7编写程序来产生下列基本脉冲序列。

（1）单位抽样序列，起点n s =0，终点n f =10，在n 0=3处有一单位脉冲。

（2）单位阶跃序列，起点n s =0，终点n f =10，在n 0=3前为0，在n s =3后为1。 （3）复指数序列，σ=-0.2，ω0=0.5。 解：程序清单如下。 clear,n0=3;ns=0;nf=10;

n1=[ns:nf];x1=[zeros(1,n0-ns),1,zeros(1,nf-n0)];%单位抽样序列的产生 n2=[ns:nf];x2=[zeros(1,n0-ns),ones(1,nf-n0+1)];%单位阶跃序列的产生 n3=[ns:nf];x3=exp((-0.2+0.5j)*n3);%复指数序列的产生

subplot(2,2,1),stem(n1,x1);title('单位脉冲序列δ(n-3)');%画图表示单位抽样序列 subplot(2,2,2),stem(n2,x2);title('单位阶跃序列u(n-3)');%画图表示单位阶跃序列 subplot(2,2,3),stem(n3,real(x3));line([0,10],[0,0])

title('复指数序列');ylabel('实部');%画图表示复指数序列的实部 subplot(2,2,4),stem(n3,imag(x3));line([0,10],[0,0]) title('复指数序列');ylabel('虚部'); 程序运行结果如图1.25所示。

图1.25几种基本脉冲序列

在数字信号处理的基本理论和MA TLAB 信号处理工具箱函数的基础上，可以自己编写一些子程序以便调用。

（1）单位抽样序列)(0n n -δ的生成函数impseq.m

function [x,n]=impseq(n0,ns,nf)

n=[ns:nf];x=[(n-n0)==0];%序列的起点为ns ，终点为nf ，在n=n0点处生成一个单位脉冲。

（2)单位阶跃序列)(0n n u -的生成函数stepseq.m function [x,n]=stepseq(n0,ns,nf)

n=[ns:nf];x=[(n-n0)>=0];%序列的起点为ns ，终点为nf ，在n=n0点处生成单位阶跃。 （3）两个信号相加的生成函数sigadd.m function [y,n]=sigadd(x1,n1,x2,n2)

n=min(min(n1),min(n2)):max(max(n1),max(n2)) y1=zeros(1,length(n));y2=y1;

y1=(find((n>=min(n1))&(n<=max(n1))==1))=x1; y2=(find((n>=min(n2))&(n<=max(n2))==1))=x2; y=y1+y2;

（4）两个信号相乘的生成函数sigmult.m function [y,n]=sigmult(x1,n1,x2,n2)

n=min(min(n1),min(n2)):max(max(n1),max(n2)) y1=zeros(1,length(n));y2=y1;

y1=(find((n>=min(n1))&(n<=max(n1))==1))=x1; y2=(find((n>=min(n2))&(n<=max(n2))==1))=x2; y=y1.*y2;

（5）序列移位)()(0n n x n y -=的生成函数sigshift.m function [y,n]=sigshift(x,m,n0) n=m+n0;y=x

（6）序列翻折)()(n x n y -=的生成函数sigfold.m function [y,n]=sigfold(x,n) y=fliplr(x);n=-fliplr(n)

（7）奇偶综合函数evenodd.m

此函数可以将任一给定的序列x(n)分解为x e (n)和x o (n)两部分。 function [xe,xo,m]=evenodd(x,n) if(imag(x)~=0)

error(‘xisnotarealsequence ’)； end

m=-fliplr(n);m1=min([m,n]);m2=max([m,n]);m=m1:m2; nm=n(1)-m(1);n1=1:length(n); x1=zeros(1,length(m)); x1(n1+nm)=x;x=x1;

xe=0.5*(x+fliplr(x));xo=0.5*(x-fliplr(x)); (8）求卷积和

求卷积和直接采用MATLAB 中的函数conv ，即y=conv(x,h)；它默认序列从n=0开始。但是如果序列是从一负值开始，即

{}{}

21:)(21:)(nh n nh n h nx n nx n x ≤

≤

≤≤

其中nx1<0或nh1<0，或两者同时为负，这样就不能直接采用conv 函数。其卷积结果序列为{}2211:)(nh nx n nh nx n y +≤≤+，这样就可构成一个新的卷积函数conv_m ，求出带下标的序列卷积。

function [y,ny]=conv_m(x,nx,h,nh)

ny1=nx(1)+nh(1);ny2=nx(length(x))+nh(length(h)); ny=[ny1:ny2]; y=conv(x,h)

例1.8产生系列序列，并绘出离散图。 （1）x1(n)=3δ(n-2)-δ(n+4)-5≤n ≤5

（2）x2(n)=n[u(n)-u(n-5)]+10e -0.2(n-5)[u(n-10)-u(n-20)]0≤n ≤20 （3）x3(n)=cos(0.04πn)+0.2w(n)0≤n ≤50

其中w(n)是均值为0，方差为1的白噪声序列。 （4）x4(n)=]1111[)(~x x x x n x =,其中x1=[10123]

图1.26产生系列序列的示意图

解：程序清单如下。

n=[-5:5];x1=3*impseq(2,-5,5)-impseq(-4,-5,5);

subplot(2,2,1);stem(n,x1);xlabel('n');ylabel('x1(n)');axis([-5,5,-2,4]); n=[0:20];x21=n.*(stepseq(0,0,20)-stepseq(10,0,20));

x22=10*exp(-0.2*(n-5)).*(stepseq(10,0,20)-stepseq(20,0,20));x2=x21+x22; subplot(2,2,2);stem(n,x2);xlabel('n');ylabel('x2(n)');axis([0,20,-1,11]) n=[0:50];x3=cos(0.04*pi*n)+0.2*randn(size(n));

subplot(2,2,3);stem(n,x3);xlabel('n');ylabel('x3(n)');axis([0,50,-2,2]) n=[-10:9];x4=[1,0,1,2,3];xtilde=x4'*ones(1,4);xtilde=(xtilde(:))';

subplot(2,2,4);stem(n,xtilde);xlabel('n');ylabel('x4(n)');axis([-10,9,-1,4]) 程序运行结果如图1.26所示。

（二）离散序列的基本运算与波形变换的MATLAB 实现

1．加法

对于离散序列来说，序列相加是将信号对应时间序号的值逐项相加，在这里不能象连续时间信号那样用符号运算来实现，而必须用向量表示的方法，即在MATLAB 中离散序列的相加需表示成两个向量的相加，因而参加运算的两序列向量必须具有相同的维数。

实现离散序列相加的MA TLAB 实用子程序如下

function [f,k]=lsxj(f1,f2,k1,k2)

%实现f(k)=f1(k)+f2(k),f1,f2,k1,k2是参加运算的二离散序列及其对应的时间序列向量,f 和k 为返回的和序列及其对应的时间序列向量

k=min(min(k1),min(k2)):max(max(k1),max(k2));%构造和序列长度 s1=zeros(1,length(k));s2=s1;%初始化新向量

s1(find((k>=min(k1))&(k<=max(k1))==1))=f1;%将f1中在和序列范围内但又无定义的点赋值为零 s2(find((k>=min(k2))&(k<=max(k2))==1))=f2;%将f2中在和序列范围内但又无定义的点赋值为零 f=s1+s2;%两长度相等序列求和 stem(k,f,'filled')

axis([(min(min(k1),min(k2))-1),(max(max(k1),max(k2))+1),(min(f)-0.5),(max(f)+0.5)])%坐标轴显示范围

图1.27离散序列相加波形图

例1.9已知两离散序列分别为

[]{}

[]{}1,1,12,1,0,1,221=--=k f k f

试用MATLAB 绘出它们的波形及[][]k f k f 21+的波形。

解：MA TLAB 的程序如下

%求两离散序列之和实现程序 f1=-2:2;k1=-2:2; f2=[111];k2=-1:1; subplot 221;

stem(k1,f1),axis([-33-2.52.5]); title('f1[k]'); subplot 222;

stem(k2,f2),axis([-33-2.52.5]); title('f2[k]'); subplot 223;

[f,k]=lsxj(f1,f2,k1,k2); title('f[k]=f1[k]+f2(k)');

程序运行的结果如图1.27所示。 2．乘法

与离散序列加法相似，这里参加运算的两序列向量必须具有相同的维数。实现离散时间信号相乘的MATLAB 实用子程序如下

function [f,k]=lsxc(f1,f2,k1,k2)

%实现f(k)=f1(k)+f2(k),f1,f2,k1,k2是参加运算的二离散序列及其对应的时间序列向量,f 和k 为返回的和序列及其对应的时间序列向量

k=min(min(k1),min(k2)):max(max(k1),max(k2));%构造和序列长度 s1=zeros(1,length(k));s2=s1;%初始化新向量

s1(find((k>=min(k1))&(k<=max(k1))==1))=f1;%将f1中在和序列范围内但又无定义的点赋值为零 s2(find((k>=min(k2))&(k<=max(k2))==1))=f2;%将f2中在和序列范围内但又无定义的点赋值为零 f=s1.*s2;%两长度相等序列求和 stem(k,f,'filled')

axis([(min(min(k1),min(k2))-1),(max(max(k1),max(k2))+1),(min(f)-0.5),(max(f)+0.5)])%坐标轴显示范围

图1.28离散序列相乘波形图

例1.10试用MATLAB 绘出例1.9中两离散序列乘法[][]k f k f 21?的波形。 解：MATLAB 的程序如下 %求两离散序列之积实现程序 f1=-2:2;k1=-2:2; f2=[111];k2=-1:1; subplot 221;

stem(k1,f1),axis([-33-2.52.5]); title('f1[k]'); subplot 222;

stem(k2,f2),axis([-33-2.52.5]); title('f2[k]'); subplot 223;

[f,k]=lsxc(f1,f2,k1,k2); title('f[k]=f1[k]*f2(k)');

程序运行的结果如图1.28所示。

3．离散序列的差分与求和

离散序列的差分[][][]1--=?k f k f k f ,在MATLAB 中用diff 函数来实现，其调用形式为

()f diff y =

离散序列的求和

[]∑=2

1

k k k k f 与信号相加运算不同，求和运算是把1

k 和2

k

之间的所有样本[]k f 加起来，

在MA TLAB中可利用sum函数来实现，其调用形式为

()2

()

y=

sum

k

f

:1k

4．离散序列的时移、反折、尺度变换

离散序列的时移、反折、尺度变换与连续时间信号相似，在此举一例来说明其MA TLAB实现过程。其MATLAB源程序如下

%离散序列图形变换实现程序

clear;

k=-12:12;

k1=2.*k+4;

f=-[stepfun(k,-3)-stepfun(k,-1)]+...

4.*[stepfun(k,-1)-stepfun(k,0)]+...

0.5*k.*[stepfun(k,0)-stepfun(k,11)];

f1=-[stepfun(k1,-3)-stepfun(k1,-1)]+...

4.*[stepfun(k1,-1)-stepfun(k1,0)]+...

0.5*k1.*[stepfun(k1,0)-stepfun(k1,11)];

subplot221;

stem(k,f);

axis([-1212-16]);

grid on;

text(-8,3,'f[k]')

subplot222;

stem(k+1,f);

axis([-1212-16]);

grid on;

text(-9.5,3,'f[k-1]')

subplot223;

stem(k,f1);

axis([-1212-16]);

grid on;

text(-8,3,'f[2k+4]')

subplot224;

stem(2-k,f);

axis([-1212-16]);

grid on;

text(5.5,3,'f[2-k]')

程序运行结果如图1.29所示。图1.29离散序列的波形变换图形

三、实验内容与步骤

1、从讲义中选择两个已有程序进行编写验证。

2、从习题与上机实验中选择两题进行编程。

四、实验仪器设备

计算机

MATLAB软件

五、实验注意事项

预先了解MATLAB的一些基础知识。

课前预先阅读并理解实验程序。

根据要求独立编程设计。

六、思考题

１.讨论复指数序列的哪个参数控制该序列的增长或衰减率？哪个参数控制该序列的振幅？

２.讨论正弦序列x3(n)的哪个参数控制该序列的相位？哪个参数控制该序列的振幅？

3.讨论算术运算符*和.*之间的区别是什么?

实验二 离散时间信号的Z 变换及快速傅立叶变换FFT 及其应用 一、实验目的

１.运用MA TLAB 理解Z 变换及其绘制H(z)的零极点图。 ２.运用MA TLAB 计算逆Z 变换。

3.利用MA TLAB 的快速傅立叶变换来计算信号的离散傅立叶变换。

4.利用MA TLAB 程序，理解进一步离散傅立叶变换的物理意义。

5.利用MA TLAB 程序，理解快速卷积算法。

二、实验原理

（一）离散时间信号的Z 变换 1．residuez:求极点留数分解 前面已经介绍，离散系统的输出

Y(z)=X(z)H(z)=B(z)/A(z)

他必然是z 的有理分式

+++-++-+-=+++++++++=

=------------1

1

11)1(1)1(1)2()1()(1)()2(1)2()1(1)1()1()()2()1()1()()2()1()()()(z k k z N p N r z p r z p r z N A z N A z A A z N B z N B z B B z A z B z Y N

N M

M

从而得出其时域信号为

+-+++++=)1()2()()1()()()()()2()2()()1()1()(n k n k n u N p N r n u p r n u p r n y n n n δδ

其中k 是当M ≥N 时的直接项，即有限序列，而其余的都是无限序列。其中的r,p,k 向量可以由下列语句求得

[r,p,k]=residuez(B,A)

2．impz ：求H(z)的Z 反变换h(n) [h,T]=impz(B,A,N)

h 为存放h(n)的列向量，时间变量N 存放在列向量T 中，当N 为标量时，表示T=[0：N-1]’，计算h(n)，n=0，1，2，…，N-1；当N 为向量时，T=N ，仅计算N 指定的整数点上的h(n)。

3．freqz ：求数字滤波器H(z)的频率响应函数 H=freqz[B,A,w]

计算由向量w 指定的数字频率点上数字滤波器H(z)的频率响应，结果存在H 向量中。

[H,w]=freqz[B,A,M]

计算出M 个频率点上的频率响应，存放在H 向量中，M 个频率存放在w 向量中，freqz 函数自动将这M 个频率点均匀设置在频率范围[0，π]之间。

若缺省w 和M 时，函数自动选取512个频率点计算。不带输出向量的freqz 函数将自动绘制幅频和相频曲线。

4．zplane:绘制H(z)的零极点图 zplane(z,p)

绘制出列向量z 中的零点（以符号“о”表示）和列向量p 中的极点（以符号“×”表示）以及参考单位圆，并在多阶零点和极点的右上角标出其阶数。如果z 和p 为矩阵，则会以不同颜色绘出z 和p 各列中的零点和极点。

zplane(B,A)

绘制出系统函数H(z)的零极点图。 （二）快速傅立叶变换FFT 及其应用

MATLAB 为计算数据的离散快速傅时叶变换，提供了一系列丰富的数学函数，主要有fft 、ifft 、fft2、

ifft2和czt 等。当所处理的数据的长度为2的幂次时，采用基-2算法进行计算，计算速度会显著增加。所以，要尽可能使所要处理的数据长度为2幂次或者用添零的方式来添补数据使之成为2的幂次。 1．fft 和ifft 函数 调用格式是： （1）()X fft Y =

如果X 是向量，则采用傅时叶变换来求解X 的离散傅里叶变换；如果X 是矩阵，则计算该矩阵每一列的离散傅里叶变换；如果X 是()D N *维数组，则是对第一个非单元素的维进行离散傅里叶变换。 （2）()N X fft Y ,=

N 是进行离散傅里叶变换的X 的数据长度，可以通过对X 进行补零或截取来实现。 （3）[]()dim ,,X fft Y =或()dim ,,N X fft Y =

在参数dim 指定的维上进行离散傅里叶变换；当X 为矩阵时，dim 用来指定变换的实施方向：dim=1，表明变换按列进行；dim=2，表明变换按行进行。

函数ifft 的参数应用与函数fft 完全相同。 2．fft2和ifft2函数

调用格式是：

（1）()X fft Y 2=

如果X 是向量，则此傅里叶变换即变成一维傅里叶变换fft ；如果X 是矩阵，则是计算该矩阵的二维快速傅里叶变换；数据二维傅里叶变换fft 2（X ）相当于()()''X fft fft ，即先对X 的列做一维傅里叶变换，然后再对变换结果的行做一维傅里叶变换。

（2）()N M X fft Y ,,2= 通过对X 进行补零或截断，使得X 成为()N M *的矩阵。 函数ifft2的参数应用与函数fft2完全相同

fftn 、ifftn 是对数据进行多维快速傅立变换，其应用与fft2、ifft2类似；在此，不再叙述。 （三）MATLAB 程序

对于离散系统其系统函数可由式（2.46）表示，则系统函数的零点和极点可以用MATLAB 的多项式求根函数roots()来实现，调用该函数的命令格式为：

p=roots(A)

其中A 为待求根的多项式的系数构成的行向量，返回向量p 则包含该多项式所有的根位置列向量。值得注意的是：求系统函数零极点时，离散系统的系统函数可能有两种形式，一种是分子和分母多项式按Z 的降幂次序排列，另一种是分子和分母多项式按1

-z 的升幂次序排列。若()z H 是以Z 的降幂次序排列，则系

数向量一定要由多项式的最高幂次开始，一直到常数项，缺项要用0补齐; 若()z H 是以按1

-z 的升幂次序

排列，则分子和分母多项式系数向量的维数一定要相同，不足的要用0补齐，不则0=z 的零点或极点就可能被漏掉。

下面是求系统函数零极点，并绘制其零极点图的MATLAB 实用函数ljdt(),该函数在绘制系统零极点图的同时，还绘出了Z 平面的单位圆。

function ljdt(A,B)

% The function to draw the pole-zero diagram for discrete system p=roots(A); %求系统极点 q=roots(B); %求系统零点

p=p'; %将极点列向量转置为行向量 q=q'; %将零点列向量转置为行向量 x=max(abs([p q 1])); %确定纵坐标范围 x=x+0.1;

y=x; %确定横坐标范围 clf hold on

axis([-x x -y y]) %确定坐标轴显示范围 w=0:pi/300:2*pi; t=exp(i*w);

plot(t) %画单位园 axis('square')

plot([-x x],[0 0]) %画横坐标轴 plot([0 0],[-y y]) %画纵坐标轴 text(0.1,x,'jIm[z]') text(y,1/10,'Re[z]')

plot(real(p),imag(p),'x') %画极点 plot(real(q),imag(q),'o') %画零点 title('pole-zero diagram for discrete system') %标注标题 hold off

例2.24 已知某离散系统的系统函数为

()1

31

45+-+=

z z z z H

试用MATLAB 求出该系统的零极点，并画出零极点分布图，判断系统是否稳定。

图2.10 例2.24系统零极点图

解 利用ljdt()函数，MATLAB 程序为 % 绘制零极点分布图的实现程序 a=[3 -1 0 0 0 1]; b=[1 1]; ljdt(a,b) p=roots(a) q=roots(b) pa=abs(p)

程序运行如果为： p =

0.7255 + 0.4633i 0.7255 - 0.4633i -0.1861 + 0.7541i -0.1861 - 0.7541i -0.7455 q =

-1 pa =

0.8608 0.8608 0.7768 0.7768 0.7455

绘制的系统零极点如图2.10所示。由图可知，该系统的所有极点均位于Z 平面的单位圆内，故该系统为稳定系统。

例2.25 已知某离散系统的系统函数为

()3

.0005.05.01

2232+--++=z z z z z z H

试用MATLAB 求出该系统的零极点，并画出零极点分布图，求系统的单位冲激响应和幅频响应，并判断系统的是否稳定。

解 MATLAB 实现程序为

% 由系统函数求解系统脉冲响应,频率响应实现程序 b=[0 1 2 1];a=[1 -0.5 -0.005 0.3]; subplot 311

zplane(b,a);xlabel('实部');ylabel('虚部'); num=[0 1 2 1];den=[1 -0.5 -0.005 0.3]; h=impz(num,den); subplot 312

stem(h);xlabel('k');title('单位脉冲响应'); [H,w]=freqz(num,den); subplot 313

plot(w/pi,abs(H)); xlabel('频率\omega'); title('频率响应')

程序运行结果如图2.11所示。该系统是稳定的。

图2.11 零极点分布图，求系统的单位冲激响应和幅频响应

利用MATLAB 分析离散系统的零极图分布与系统单位响应时域特性的关系

例2.26 已知离散系统的零极分布分别如图2.12所示，其中虚线表示单位圆，试用MATLAB 分析系统单位响应()n h 的时域特性。

解 因系统的零极点分布图已知，则系统的系统函数()z H 就可知，故可以利用MATLAB 函数impz()函数可求出系统单位响应()n h 。

对于图2.12所示的系统，其系统函数分别为

()11-=

z z H 、()8.01-=z z H 、()2

1-=z z H ()2

28.04

cos

8.021

+**-=

z z z H π

、()1

8

cos

8.021

2

+**-=

z z z H 、

()2

22.14

cos 2.121

+**-=

z z z H

则其单位响应的时域波形MATLAB 程序为

% 零极点分布与单位响应的关系实现程序 a=[1,-1]; b=[1];

subplot 321 impz(b,a); a1=[1,-0.8]; b1=[1]; subplot 322 impz(b1,a1,10); a2=[1,-2]; b2=[1]; subplot 323 impz(b2,a2,10);

a3=[1,-2*0.8*cos(pi/4),0.8^2]; b3=[1]; subplot 324 impz(b3,a3,20);

a4=[1,-2*0.8*cos(pi/8),1]; b4=[1]; subplot 325 impz(b4,a4,20);

a5=[1,-2*1.2*cos(pi/4),1.2^2]; b5=[1]; subplot 326 impz(b5,a5,20);

程序运行结果如图2.13所示。

图2.12 离散系统的零极点分布图

由此可知：离散系统的单位响应()n h 的时域特性完全由系统函数()z H 的极点位置决定，()z H 位于Z 平面单位圆内的极点决定了()n h 随时间衰减的信号分量，位于Z 平面单位圆上的极点决定了单位响应的稳定信号分量，位于Z 平面单位圆外的极点决定了单位响应随时间增长的信号分量。

图2.13 系统单位响应的时域波形图

例2.27 求一因果线性移不变系统y(n)=0.81y(n-2)+x(n)-x(n-2)的单位抽样响应h(n)，单位阶跃响

应g(n)，及绘出H(e j ω

)的幅频和相频特性。

解:程序清单如下

b=[1,0,-1];a=[1,0,-0.81]; figure(1)

subplot(2,1,1);dimpulse(b,a,50);ylabel(‘h(n)’); subplot(2,1,2);dstep(b,a,50);ylabel(‘g(n)’); figure(2)

w=[0:1:500]*pi/500; freqz(b,a,w)

程序运行结果如图2.14所示。

图2.14 系统的抽样响应、阶跃响应及幅频、相频曲线

例2.28 梳状滤波器零极点和幅频特性。梳状滤波器系统函数有如下两种类型：

FIR 型：N

z z H --=1)(1 IIR 型：N

N N

z a z z H ----=11)(2

分别令N=8，a=0.8,0.9,0.98计算并图示)(1z H 和)(2z H 的零、极点图及幅频特性，说明极点位置的影响。

解:程序清单如下

b=[1,0,0,0,0,0,0,0,-1]; %H1(z)和H2(z)的分子多项式系数向量

a0=1; %H1(z)分母多项式系数向量

a1=[1,0,0,0,0,0,0,0,-(0.8)^8]; %H2(z)分母多项式系数向量（a=0.8） a2=[1,0,0,0,0,0,0,0,-(0.9)^8]; %H2(z)分母多项式系数向量（a=0.9） a3=[1,0,0,0,0,0,0,0,-(0.98)^8]; %H2(z)分母多项式系数向量（a=0.98） [H,w]=freqz(b,a0); %H1(z)的频响函数 [H1,w1]=freqz(b,a1); [H2,w2]=freqz(b,a2); [H3,w3]=freqz(b,a3);

subplot(4,2,1);zplane(b,a0); xlabel('实部');ylabel('虚部');title('FIR 梳状滤波器零点图') subplot(4,2,2);zplane(b,a1);xlabel('实部');ylabel('虚部'); title('IIR 梳状滤波器零极点图 a=0.8')

subplot(4,2,3);plot(w/pi,abs(H));title('FIR 梳状滤波器幅频响应曲线')

subplot(4,2,4); plot(w/pi,abs(H1));title('IIR 梳状滤波器幅频响应曲线 a=0.8')

程序运行结果如图2.15所示。

由图看出，阶数相同的时候，IIR 滤波器具有更平坦的通带特性和更窄的过渡带，极点距离单位圆越近，这个特性就越明显。

图2.15 八阶梳状滤波器零、极点图及幅频响应曲线

利用系统函数求解离散系统差分方程的MATLAB

例2.29 求解差分方程()()()()()()214.045.0245.014.0---+=----n x n x n x n y n y n y ，其中，

()()n n x n ε8.0=，初始状态()01=-y ，()12=-y ，()11=-x ，()22=-x 。

解 将方程两边进行z 变换得

数字信号处理实验二报告

实验二 IIR数字滤波器设计及软件实现 1．实验目的 （1）熟悉用双线性变换法设计IIR数字滤波器的原理与方法； （2）学会调用MATLAB信号处理工具箱中滤波器设计函数（或滤波器设计分析工具fdatool）设计各种IIR数字滤波器，学会根据滤波需求确定滤波器指标参数。 （3）掌握IIR数字滤波器的MATLAB实现方法。 （3）通过观察滤波器输入输出信号的时域波形及其频谱，建立数字滤波的概念。 2．实验原理 设计IIR数字滤波器一般采用间接法（脉冲响应不变法和双线性变换法），应用最广泛的是双线性变换法。基本设计过程是：①先将给定的数字滤波器的指标转换成过渡模拟滤波器的指标；②设计过渡模拟滤波器；③将过渡模拟滤波器系统函数转换成数字滤波器的系统函数。MATLAB信号处理工具箱中的各种IIR数字滤波器设计函数都是采用双线性变换法。第六章介绍的滤波器设计函数butter、cheby1 、cheby2 和ellip可以分别被调用来直接设计巴特沃斯、切比雪夫1、切比雪夫2和椭圆模拟和数字滤波器。本实验要求读者调用如上函数直接设计IIR数字滤波器。 本实验的数字滤波器的MATLAB实现是指调用MATLAB信号处理工具箱函数filter对给定的输入信号x(n)进行滤波，得到滤波后的输出信号y(n）。 3. 实验内容及步骤 （1）调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st，该函数还会自动绘图显示st的时域波形和幅频特性曲线，如图1所示。由图可见，三路信号时域混叠无法在时域分离。但频域是分离的，所以可以通过滤波的方法在频域分离，这就是本实验的目的。 图1 三路调幅信号st的时域波形和幅频特性曲线 （2）要求将st中三路调幅信号分离，通过观察st的幅频特性曲线，分别确定可以分离st中三路抑制载波单频调幅信号的三个滤波器（低通滤波器、带通滤波器、高通滤波器）的通带截止频率和阻带截止频率。要求滤波器的通带最大衰减为0.1dB,阻带最小衰减为

数字信号处理实验一

一、实验目的 1. 通过本次实验回忆并熟悉MATLAB这个软件。 2. 通过本次实验学会如何利用MATLAB进行序列的简单运算。 3. 通过本次实验深刻理解理论课上的数字信号处理的一个常见方法——对时刻n的样本附近的一些样本求平均，产生所需的输出信号。 3. 通过振幅调制信号的产生来理解载波信号与调制信号之间的关系。 二、实验内容 1. 编写程序在MATLAB中实现从被加性噪声污染的信号中移除噪声的算法，本次试验采用三点滑动平均算法，可直接输入程序P1.5。 2. 通过运行程序得出的结果回答习题Q1.31-Q1.33的问题，加深对算法思想的理解。 3. 编写程序在MATLAB中实现振幅调制信号产生的算法，可直接输入程序P1.6。 4. 通过运行程序得出的结果回答习题Q1.34-Q1.35的问题，加深对算法思想的理解。 三、主要算法与程序 1. 三点滑动平均算法的核心程序： %程序P1.5 %通过平均的信号平滑 clf; R=51; d=0.8*(rand(R,1)-0.5);%产生随噪声 m=0:R-1; s=2*m.*(0.9.^m);%产生为污染的信号 x=s+d';%产生被噪音污染的信号 subplot(2,1,1); plot(m,d','r-',m,s,'g--',m,x,'b-.');

xlabel('时间序号n');ylabel('振幅'); legend('d[n]','s[n]','x[n]'); x1=[0 0 x];x2=[0 x 0];x3=[x 0 0]; y=(x1+x2+x3)/3; subplot(2,1,2); plot(m,y(2:R+1),'r-',m,s,'g--'); legend('y[n]','s[n]'); xlabel('时间序号n');ylabel('振幅'); 2. 振幅调制信号的产生核心程序：（由于要几个结果，因此利用subplot函数画图） %程序P1.6 %振幅调制信号的产生 n=0:100; m=0.1;fH=0.1;fL=0.01; m1=0.3;fH1=0.3;fL1=0.03; xH=sin(2*pi*fH*n); xL=sin(2*pi*fL*n); y=(1+m*xL).*xH; xH1=sin(2*pi*fH1*n); xL1=sin(2*pi*fL1*n); y1=(1+m1*xL).*xH; y2=(1+m*xL).*xH1; y3=(1+m*xL1).*xH; subplot(2,2,1); stem(n,y); grid; xlabel('时间序号n');ylabel('振幅');title('m=0.1;fH=0.1;fL=0.01;'); subplot(2,2,2); stem(n,y1); grid; xlabel('时间序号n');ylabel('振幅');title('m=0.3;fH=0.1;fL=0.01;'); subplot(2,2,3); stem(n,y2); grid; xlabel('时间序号n');ylabel('振幅');title('m=0.3;fH=0.3;fL=0.01;'); subplot(2,2,4); stem(n,y3); grid;

数字信号处理基础实验指导书

《数字信号处理》实验指导书 光电工程学院二○○九年十月

实验一离散时间信号分析 一、实验目的 1．掌握各种常用的序列，理解其数学表达式和波形表示。 2．掌握在计算机中生成及绘制数字信号波形的方法。 3．掌握序列的相加、相乘、移位、反转等基本运算及计算机实现与作用。 4．掌握线性卷积软件实现的方法。 5．掌握计算机的使用方法和常用系统软件及应用软件的使用。 6．通过编程，上机调试程序，进一步增强使用计算机解决问题的能力。 二、实验原理 1．序列的基本概念 离散时间信号在数学上可用时间序列来表示，其中代表序列的第n个数字，n代表时间的序列，n的取值范围为的整数，n取其它值没有意义。离散时间信号可以是由模拟信号通过采样得到，例如对模拟信号进行等间隔采样，采样间隔为T，得到一个有序的数字序列就是离散时间信号，简称序列。 2．常用序列 常用序列有：单位脉冲序列（单位抽样）、单位阶跃序列、矩形序列、实指数序列、复指数序列、正弦型序列等。 3．序列的基本运算 序列的运算包括移位、反转、和、积、标乘、累加、差分运算等。 4．序列的卷积运算 上式的运算关系称为卷积运算，式中代表两个序列卷积运算。两个序列的卷积是一个序列与另一个序列反褶后逐次移位乘积之和，故称为离散卷积，也称两序列的线性卷积。其计算的过程包括以下4个步骤。 （1）反褶：先将和的变量换成，变成和，再将以纵轴为对称轴反褶成。 （2）移位：将移位，得。当为正数时，右移位；当为负数时，左

移位。 （3）相乘：将和的对应点值相乘。 （4）求和：将以上所有对应点的乘积累加起来，即得。 三、主要实验仪器及材料 微型计算机、Matlab软件6.5或更高版本。 四、实验内容 1．知识准备 认真复习以上基础理论，理解本实验所用到的实验原理。 2．离散时间信号（序列）的产生 利用MATLAB或C语言编程产生和绘制下列有限长序列： （1）单位脉冲序列 （2）单位阶跃序列 （3）矩形序列 （4）正弦型序列 （5）任意序列 3．序列的运算 利用MATLAB编程完成上述两序列的移位、反转、加法、乘法等运算，并绘制运算后序列的波形。 4．卷积运算 利用MATLAB编制一个计算两个序列线性卷积的通用程序，计算上述两序列，并绘制卷积后序列的波形。 5．上机调试并打印或记录实验结果。 6．完成实验报告。 五、实验报告要求 1. 简述实验原理及目的。 2. 给出上述序列的实验结果。 3. 列出计算卷积的公式，画出程序框图，并列出实验程序清单 （可略）（包括必要的程序说明）。 4. 记录调试运行情况及所遇问题的解决方法。 5. 给出实验结果，并对结果做出分析。 6. 简要回答思考题。 1 如何产生方波信号序列和锯齿波信号序列? 2 实验中所产生的正弦序列的频率是多少?是否是周期序列?

数字信号处理实验程序2.

2.1 clc close all; n=0:15; p=8;q=2; x=exp(-(n-p.^2/q; figure(1; subplot(3,1,1; stem(n,x; title('exp(-(n-p^2/q,p=8,q=2'; xk1=fft(x,16; q=4; x=exp(-(n-p.^2/q; subplot(3,1,2; xk2=fft(x,16; stem(n,x; title('exp(-(n-p^2/q,p=8,q=4'; q=8; x=exp(-(n-p.^2/q;

xk3=fft(x,16; subplot(3,1,3; stem(n,x; title('exp(-(n-p^2/q,p=8,q=8';%时域特性figure(2; subplot(3,1,1; stem(n,abs(xk1; title('exp(-(n-p^2/q,p=8,q=2'; subplot(3,1,2; stem(n,abs(xk2; title('exp(-(n-p^2/q,p=8,q=4'; subplot(3,1,3; stem(n,abs(xk3; title('exp(-(n-p^2/q,p=8,q=8';%频域特性%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%% p=8;q=8; figure(3; subplot(3,1,1; stem(n,x; title('exp(-(n-p^2/q,p=8,q=8';

xk1=fft(x,16; p=13; x=exp(-(n-p.^2/q; subplot(3,1,2; xk2=fft(x,16; stem(n,x; title('exp(-(n-p^2/q,p=13,q=8'; p=14; x=exp(-(n-p.^2/q; xk3=fft(x,16; subplot(3,1,3; stem(n,x; title('exp(-(n-p^2/q,p=14,q=8';%时域特性figure(4; subplot(3,1,1; stem(n,abs(xk1; title('exp(-(n-p^2/q,p=8,q=8'; subplot(3,1,2; stem(n,abs(xk2; title('exp(-(n-p^2/q,p=13,q=8'; subplot(3,1,3;

数字信号处理实验作业

实验6 数字滤波器的网络结构 一、实验目的： 1、加深对数字滤波器分类与结构的了解。 2、明确数字滤波器的基本结构及其相互间的转换方法。 3、掌握用MA TLAB 语言进行数字滤波器结构间相互转换的子函数及程序编写方法。 二、实验原理： 1、数字滤波器的分类 离散LSI 系统对信号的响应过程实际上就是对信号进行滤波的过程。因此，离散LSI 系统又称为数字滤波器。 数字滤波器从滤波功能上可以分为低通、高通、带通、带阻以及全通滤波器；根据单位脉冲响应的特性，又可以分为有限长单位脉冲响应滤波器（FIR ）和无限长单位脉冲响应滤波器（IIR ）。 一个离散LSI 系统可以用系统函数来表示： M -m -1-2-m m m=0 012m N -1-2-k -k 12k k k=1 b z b +b z +b z ++b z Y(z)b(z)H(z)=== =X(z)a(z) 1+a z +a z ++a z 1+a z ∑∑ 也可以用差分方程来表示： N M k m k=1 m=0 y(n)+a y(n-k)=b x(n-m)∑∑ 以上两个公式中，当a k 至少有一个不为0时，则在有限Z 平面上存在极点，表达的是以一个IIR 数字滤波器；当a k 全都为0时，系统不存在极点，表达的是一个FIR 数字滤波器。FIR 数字滤波器可以看成是IIR 数字滤波器的a k 全都为0时的一个特例。 IIR 数字滤波器的基本结构分为直接Ⅰ型、直接Ⅱ型、直接Ⅲ型、级联型和并联型。 FIR 数字滤波器的基本结构分为横截型（又称直接型或卷积型）、级联型、线性相位型及频率采样型等。本实验对线性相位型及频率采样型不做讨论，见实验10、12。 另外，滤波器的一种新型结构——格型结构也逐步投入应用，有全零点FIR 系统格型结构、全极点IIR 系统格型结构以及全零极点IIR 系统格型结构。 2、IIR 数字滤波器的基本结构与实现 （1）直接型与级联型、并联型的转换 例6-1 已知一个系统的传递函数为 -1-2-3 -1-2-3 8-4z +11z -2z H(z)=1-1.25z +0.75z -0.125z 将其从直接型（其信号流图如图6-1所示）转换为级联型和并联型。

数字信号处理实验1,2,3,4

实验一 连续时间系统的时域和频域分析相关MATLAB 函数1.设描述连续时间系统的微分方程为：)()()()()()()()(01)1(1)(01)1(1)(t f b t f b t f b t f b t y a t y a t y a t y a m m m m n n n n +'+++=+'+++---- 则可用向量和表示该系统，即 a b ] ,,,,[011a a a a a n n -=],,,,[011b b b b b m m -=注意，向量和的元素一定要以微分方程时间求导的降幂次序排列，且缺项要用0补齐。a b 如微分方程)()()(2)(3)(t f t f t y t y t y +''=+'+''表示该系统的向量为 ]2 3 1[=a ]1 0 1[=b （1）求解冲激响应：impulse()函数impulse()函数有以下四种调用格式： ① impulse(b,a) 该调用格式以默认方式绘制由向量和定义的连续时间系统的冲激响应的时域波形。a b ② impulse(b,a,t)该调用格式绘制由向量和定义的连续时间系统在时间范围内的冲激响应的时a b t ~0域波形。③ impulse(b,a, t1:p:t2)该调用格式绘制由向量和定义的连续时间系统在时间范围内，且以时间间a b 21~t t 隔均匀抽样的冲激响应的时域波形。p ④ y=impulse(b,a,t1:p:t2)该调用格式并不绘制系统冲激响应的波形，而是求出由向量和定义的连续时间系a b 统在时间范围内以时间间隔均匀抽样的系统冲激响应的数值解。21~t t p （2）求解阶跃响应：step()函数 step()函数也有四种调用格式：① step(b,a) ② step(b,a,t) ③ step(b,a, t1:p:t2) ④ y=step(b,a,t1:p:t2) 上述调用格式的功能与impulse()函数完全相同。

数字信号处理实验作业

实验5 抽样定理 一、实验目的： 1、了解用MA TLAB 语言进行时域、频域抽样及信号重建的方法。 2、进一步加深对时域、频域抽样定理的基本原理的理解。 3、观察信号抽样与恢复的图形，掌握采样频率的确定方法和插公式的编程方法。 二、实验原理： 1、时域抽样与信号的重建 （1）对连续信号进行采样 例5-1 已知一个连续时间信号sin sin(),1Hz 3 ππ=0001f(t)=(2f t)+6f t f ，取最高有限带宽频率f m =5f 0，分别显示原连续时间信号波形和F s >2f m 、F s =2f m 、F s <2f m 三情况下抽样信号的波形。 程序清单如下： %分别取Fs=fm ，Fs=2fm ，Fs=3fm 来研究问题 dt=0.1; f0=1; T0=1/f0; m=5*f0; Tm=1/fm; t=-2:dt:2; f=sin(2*pi*f0*t)+1/3*sin(6*pi*f0*t); subplot(4,1,1); plot(t,f); axis([min(t),max(t),1.1*min(f),1.1*max(f)]); title('原连续信号和抽样信号'); for i=1:3; fs=i*fm;Ts=1/fs; n=-2:Ts:2; f=sin(2*pi*f0*n)+1/3*sin(6*pi*f0*n); subplot(4,1,i+1);stem(n,f,'filled'); axis([min(n),max(n),1.1*min(f),1.1*max(f)]); end 程序运行结果如图5-1所示：

原连续信号和抽样信号 图5-1 （2）连续信号和抽样信号的频谱 由理论分析可知，信号的频谱图可以很直观地反映出抽样信号能否恢复原模拟信号。因此，我们对上述三种情况下的时域信号求幅度谱，来进一步分析和验证时域抽样定理。 例5-2编程求解例5-1中连续信号及其三种抽样频率（F s>2f m、F s=2f m、F s<2f m）下的抽样信号的幅度谱。 程序清单如下： dt=0.1;f0=1;T0=1/f0;fm=5*f0;Tm=1/fm; t=-2:dt:2;N=length(t); f=sin(2*pi*f0*t)+1/3*sin(6*pi*f0*t); wm=2*pi*fm;k=0:N-1;w1=k*wm/N; F1=f*exp(-j*t'*w1)*dt;subplot(4,1,1);plot(w1/(2*pi),abs(F1)); axis([0,max(4*fm),1.1*min(abs(F1)),1.1*max(abs(F1))]); for i=1:3; if i<=2 c=0;else c=1;end fs=(i+c)*fm;Ts=1/fs; n=-2:Ts:2;N=length(n); f=sin(2*pi*f0*n)+1/3*sin(6*pi*f0*n); wm=2*pi*fs;k=0:N-1; w=k*wm/N;F=f*exp(-j*n'*w)*Ts; subplot(4,1,i+1);plot(w/(2*pi),abs(F)); axis([0,max(4*fm),1.1*min(abs(F)),1.1*max(abs(F))]); end 程序运行结果如图5-2所示。 由图可见，当满足F s≥2f m条件时，抽样信号的频谱没有混叠现象；当不满足F s≥2f m 条件时，抽样信号的频谱发生了混叠，即图5-2的第二行F s＜2f m的频谱图，，在f m=5f0的围，频谱出现了镜像对称的部分。

数字信号处理实验报告(实验1_4)

实验一 MATLAB 仿真软件的基本操作命令和使用方法 实验容 1、帮助命令 使用 help 命令，查找 sqrt （开方）函数的使用方法； 2、MATLAB 命令窗口 （1）在MATLAB 命令窗口直接输入命令行计算3 1)5.0sin(21+=πy 的值； （2）求多项式 p(x) = x3 + 2x+ 4的根； 3、矩阵运算 （1）矩阵的乘法 已知 A=[1 2;3 4]， B=[5 5;7 8]，求 A^2*B

（2）矩阵的行列式 已知A=[1 2 3;4 5 6;7 8 9]，求A （3）矩阵的转置及共轭转置 已知A=[1 2 3;4 5 6;7 8 9]，求A' 已知B=[5+i,2-i,1;6*i,4,9-i], 求B.' , B' （4）特征值、特征向量、特征多项式 已知A=[1.2 3 5 0.9;5 1.7 5 6;3 9 0 1;1 2 3 4] ，求矩阵A的特征值、特征向量、特征多项式；

（5）使用冒号选出指定元素 已知：A=[1 2 3;4 5 6;7 8 9]；求A 中第3 列前2 个元素；A 中所有列第2，3 行的元素； 4、Matlab 基本编程方法 （1）编写命令文件：计算1+2+…+n<2000 时的最大n 值；

（2）编写函数文件：分别用for 和while 循环结构编写程序，求 2 的0 到15 次幂的和。

5、MATLAB基本绘图命令 （1）绘制余弦曲线 y=cos(t)，t∈[0，2π]

（2）在同一坐标系中绘制余弦曲线 y=cos(t-0.25)和正弦曲线 y=sin(t-0.5)， t∈[0，2π] （3）绘制[0，4π]区间上的 x1=10sint 曲线，并要求： （a）线形为点划线、颜色为红色、数据点标记为加号； （b）坐标轴控制：显示围、刻度线、比例、网络线 （c）标注控制：坐标轴名称、标题、相应文本； >> clear;

数字信号处理实验报告一

武汉工程大学 数字信号处理实验报告 姓名：周权 学号：1204140228 班级：通信工程02

一、实验设备 计算机，MATLAB语言环境。 二、实验基础理论 1.序列的相关概念 2.常见序列 3.序列的基本运算 4.离散傅里叶变换的相关概念 5.Z变换的相关概念 三、实验内容与步骤 1.离散时间信号（序列）的产生 利用MATLAB语言编程产生和绘制单位样值信号、单位阶跃序列、指数序列、正弦序列及随机离散信号的波形表示。 四实验目的 认识常用的各种信号，理解其数字表达式和波形表示，掌握在计算机中生成及绘制数字信号波形的方法，掌握序列的简单运算及计算机实现与作用，理解离散时间傅里叶变换，Z变换及它们的性质和信号的频域分

实验一离散时间信号（序列）的产生 代码一 单位样值 x=2; y=1; stem(x,y); title('单位样值 ') 单位阶跃序列 n0=0; n1=-10; n2=10; n=[n1:n2]; x=[(n-n0)>=0]; stem(n,x); xlabel('n'); ylabel('x{n}'); title('单位阶跃序列');

实指数序列 n=[0:10]; x=(0.5).^n; stem(n,x); xlabel('n'); ylabel('x{n}'); title('实指数序列');

正弦序列 n=[-100:100]; x=2*sin(0.05*pi*n); stem(n,x); xlabel('n'); ylabel('x{n}'); title('正弦序列');

随机序列 n=[1:10]; x=rand(1,10); subplot(221); stem(n,x); xlabel('n'); ylabel('x{n}'); title('随机序列');

数字信号处理实验1认识实验

实验1认识实验-MATLAB语言上机操作实践 一、实验目的 ㈠了解MATLAB语言的主要特点、作用。 ㈡学会MATLAB主界面简单的操作使用方法。 ㈢学习简单的数组赋值、运算、绘图、流程控制编程。 二、实验原理 ㈠简单的数组赋值方法 MATLAB中的变量和常量都可以是数组（或矩阵）,且每个元素都可以是复数。 在MATLAB指令窗口输入数组A=[1 2 3；4 5 6；7 8 9]，观察输出结果。然后，键入：A（4，2）= 11 键入：A (5，：) = [-13 -14 -15] 键入：A（4，3）= abs (A(5，1)) 键入：A ([2，5]，：) = [ ] 键入：A/2 键入：A (4，：) = [sqrt(3) (4+5)/6*2 –7] 观察以上各输出结果。将A式中分号改为空格或逗号，情况又如何？请在每式的后面标注其含义。 2．在MATLAB指令窗口输入B=[1+2i，3+4i；5+6i ，7+8i], 观察输出结果。 键入：C=[1，3；5，7]+[2，4；6，8]*i，观察输出结果。 如果C式中i前的*号省略，结果如何？ 键入：D = sqrt (2+3i) 键入：D*D 键入：E = C’, F = conj(C), G = conj(C)’ 观察以上各输出结果, 请在每式的后面标注其含义。 3．在MATLAB指令窗口输入H1=ones(3,2)，H2=zeros(2,3)，H3=eye(4)，观察输出结果。 ㈡、数组的基本运算 1．输入A=[1 3 5]，B= [2 4 6]，求C=A+B，D=A-2，E=B-A 2．求F1=A*3，F2=A.*B，F3=A./B，F4=A.\B, F5=B.\A, F6=B.^A, F7=2./B, F8=B.\2 *3．求B'，Z1=A*B’，Z2=B’*A 观察以上各输出结果,比较各种运算的区别，理解其含义。 ㈢、常用函数及相应的信号波形显示 例1：显示曲线f(t)=2sin(2πt)，（t>0） ⅰ点击空白文档图标（New M-file），打开文本编辑器。 ⅱ键入：t=0:0.01:3; (1) f=2*sin(2*pi*t); (2) plot(t,f); title(‘f（t）-t曲线’)； xlabel（‘t’），ylabel（‘f（t）’）；

数字信号处理实验答案完整版

数字信号处理实验答案 HEN system office room 【HEN16H-HENS2AHENS8Q8-HENH1688】

实验一熟悉Matlab环境 一、实验目的 1.熟悉MATLAB的主要操作命令。 2.学会简单的矩阵输入和数据读写。 3.掌握简单的绘图命令。 4.用MATLAB编程并学会创建函数。 5.观察离散系统的频率响应。 二、实验内容 认真阅读本章附录，在MATLAB环境下重新做一遍附录中的例子，体会各条命令的含义。在熟悉了MATLAB基本命令的基础上，完成以下实验。 上机实验内容： （1）数组的加、减、乘、除和乘方运算。输入A=[1 2 3 4]，B=[3 4 5 6]，求C=A+B，D=A-B，E=A.*B，F=A./B，G=A.^B并用stem语句画出A、B、C、D、E、F、G。 clear all; a=[1 2 3 4]; b=[3 4 5 6]; c=a+b; d=a-b; e=a.*b; f=a./b; g=a.^b; n=1:4; subplot(4,2,1);stem(n,a); xlabel('n');xlim([0 5]);ylabel('A'); subplot(4,2,2);stem(n,b); xlabel('n');xlim([0 5]);ylabel('B'); subplot(4,2,3);stem(n,c); xlabel('n');xlim([0 5]);ylabel('C'); subplot(4,2,4);stem(n,d); xlabel('n');xlim([0 5]);ylabel('D'); subplot(4,2,5);stem(n,e); xlabel('n');xlim([0 5]);ylabel('E'); subplot(4,2,6);stem(n,f); xlabel('n');xlim([0 5]);ylabel('F'); subplot(4,2,7);stem(n,g); xlabel('n');xlim([0 5]);ylabel('G'); （2）用MATLAB实现下列序列： a) x(n)= 0≤n≤15 b) x(n)=e+3j)n 0≤n≤15 c) x(n)=3cosπn+π)+2sinπn+π) 0≤n≤15 d) 将c)中的x(n)扩展为以16为周期的函数x(n)=x(n+16),绘出四个周期。

数字信号处理实验及参考程序

数字信号处理实验实验一离散时间信号与系统及MA TLAB实现 1.单位冲激信号： n = -5:5; x = (n==0); subplot(122); stem(n, x); 2.单位阶跃信号： x=zeros(1,11); n0=0; n1=-5; n2=5; n = n1:n2; x(:,n+6) = ((n-n0)>=0); stem(n,x); 3.正弦序列： n = 0:1/3200:1/100; x=3*sin(200*pi*n+1.2); stem(n,x); 4.指数序列 n = 0:1/2:10; x1= 3*(0.7.^n); x2=3*exp((0.7+j*314)*n); subplot(221); stem(n,x1); subplot(222); stem(n,x2); 5.信号延迟 n=0:20; Y1=sin(100*n); Y2=sin(100*(n-3)); subplot(221); stem(n,Y1); subplot(222); stem(n,Y2);

6.信号相加 X1=[2 0.5 0.9 1 0 0 0 0]; X2=[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7]; X=X1+X2; stem(X); 7.信号翻转 X1=[2 0.5 0.9 1]; n=1:4; X2=X1(5-n); subplot(221); stem(n,X1); subplot(222); stem(n,X2); 8.用MATLAB计算序列{-2 0 1 –1 3}和序列{1 2 0 -1}的离散卷积。a=[-2 0 1 -1 3]; b=[1 2 0 -1]; c=conv(a,b); M=length(c)-1; n=0:1:M; stem(n,c); xlabel('n'); ylabel('幅度'); 9．用MA TLAB计算差分方程 当输入序列为时的输出结果。 N=41; a=[0.8 -0.44 0.36 0.22]; b=[1 0.7 -0.45 -0.6]; x=[1 zeros(1,N-1)]; k=0:1:N-1; y=filter(a,b,x); stem(k,y) xlabel('n'); ylabel('幅度') 10.冲激响应impz N=64; a=[0.8 -0.44 0.36 0.22];

数字信号处理实验一

实验一 离散时间信号分析 班级 信息131班 学号 201312030103 姓名 陈娇 日期 一、实验目的 掌握两个序列的相加、相乘、移位、反褶、卷积等基本运算。 二、实验原理 1．序列的基本概念 离散时间信号在数学上可用时间序列)}({n x 来表示，其中)(n x 代表序列的第n 个数字，n 代表时间的序列，n 的取值范围为+∞<<∞-n 的整数，n 取其它值)(n x 没有意义。离散时间信号可以是由模拟信号通过采样得到，例如对模拟信号)(t x a 进行等间隔采样，采样间隔为T ，得到)}({nT x a 一个有序的数字序列就是离散时间信号，简称序列。 2．常用序列 常用序列有：单位脉冲序列（单位抽样）） （n δ、单位阶跃序列)(n u 、矩形序列)(n R N 、实指数序列、复指数序列、正弦型序列等。 3．序列的基本运算 序列的运算包括移位、反褶、和、积、标乘、累加、差分运算等。 4．序列的卷积运算 ∑∞ -∞==-= m n h n x m n h m x n y )(*)()()()( 上式的运算关系称为卷积运算，式中代表两个序列卷积运算。两个序列的卷积是一个序列与另一个序列反褶后逐次移位乘积之和，故称为离散卷积，也称两序列的线性卷积。其计算的过程包括以下4个步骤。 （1）反褶：先将)(n x 和)(n h 的变量n 换成m ，变成)(m x 和)(m h ，再将)(m h 以纵轴为对称轴反褶成)(m h -。

（2）移位：将)(m h -移位n ，得)(m n h -。当n 为正数时，右移n 位；当n 为负数时，左移n 位。 （3）相乘：将)(m n h -和)(m x 的对应点值相乘。 （4）求和：将以上所有对应点的乘积累加起来，即得)(n y 。 三、主要实验仪器及材料 微型计算机、Matlab6.5 教学版、TC 编程环境。 四、实验内容 （1）用Matlab 或C 语言编制两个序列的相加、相乘、移位、反褶、卷积等的程序； （2）画出两个序列运算以后的图形； （3）对结果进行分析； （4）完成实验报告。 五、实验结果 六、实验总结

实验一 基于Matlab的数字信号处理基本

实验一 基于Matlab 的数字信号处理基本操作 一、 实验目的：学会运用MA TLAB 表示的常用离散时间信号；学会运用MA TLAB 实现离 散时间信号的基本运算。 二、 实验仪器：电脑一台，MATLAB6.5或更高级版本软件一套。 三、 实验内容： （一） 离散时间信号在MATLAB 中的表示 离散时间信号是指在离散时刻才有定义的信号，简称离散信号，或者序列。离散序列通常用)(n x 来表示，自变量必须是整数。 离散时间信号的波形绘制在MATLAB 中一般用stem 函数。stem 函数的基本用法和plot 函数一样，它绘制的波形图的每个样本点上有一个小圆圈，默认是空心的。如果要实心，需使用参数“fill ”、“filled ”，或者参数“.”。由于MATLAB 中矩阵元素的个数有限，所以MA TLAB 只能表示一定时间范围内有限长度的序列；而对于无限序列，也只能在一定时间范围内表示出来。类似于连续时间信号，离散时间信号也有一些典型的离散时间信号。 1. 单位取样序列 单位取样序列)(n δ，也称为单位冲激序列，定义为 ) 0() 0(0 1)(≠=?? ?=n n n δ 要注意，单位冲激序列不是单位冲激函数的简单离散抽样，它在n =0处是取确定的值1。在MATLAB 中，冲激序列可以通过编写以下的impDT .m 文件来实现，即 function y=impDT(n) y=(n==0); %当参数为0时冲激为1,否则为0 调用该函数时n 必须为整数或整数向量。 【实例1-1】 利用MATLAB 的impDT 函数绘出单位冲激序列的波形图。 解：MATLAB 源程序为 >>n=-3:3; >>x=impDT(n); >>stem(n,x,'fill'),xlabel('n'),grid on >>title('单位冲激序列') >>axis([-3 3 -0.1 1.1]) 程序运行结果如图1-1所示。 图1-1 单位冲激序列

数字信号处理基础实验报告_

本科生实验报告 实验课程数字信号处理基础 学院名称地球物理学院 专业名称地球物理学 学生姓名 学生学号 指导教师王山山 实验地点5417 实验成绩 二〇一四年十一月二〇一四年十二月

填写说明 1、适用于本科生所有的实验报告（印制实验报告册除外）； 2、专业填写为专业全称，有专业方向的用小括号标明； 3、格式要求： ①用A4纸双面打印（封面双面打印）或在A4大小纸上用蓝黑色水笔书写。 ②打印排版：正文用宋体小四号，1.5倍行距，页边距采取默认形式（上下2.54cm， 左右2.54cm，页眉1.5cm，页脚1.75cm）。字符间距为默认值（缩放100%，间距：标准）；页码用小五号字底端居中。 ③具体要求： 题目（二号黑体居中）； 摘要（“摘要”二字用小二号黑体居中，隔行书写摘要的文字部分，小4号宋体）； 关键词（隔行顶格书写“关键词”三字，提炼3-5个关键词，用分号隔开，小4号黑体)； 正文部分采用三级标题； 第1章××(小二号黑体居中，段前0.5行) 1.1 ×××××小三号黑体×××××（段前、段后0.5行） 1.1.1小四号黑体（段前、段后0.5行） 参考文献（黑体小二号居中，段前0.5行），参考文献用五号宋体，参照《参考文献著录规则（GB/T 7714－2005）》。

实验一生成离散信号并计算其振幅谱 并将信号进行奇偶分解 一、实验原理 单位脉冲响应h(t)=exp(-a*t*t)*sin(2*3.14*f*t)进行离散抽样，分别得到t=0.002s，0.009s，0.011s采样的结果。用Excel软件绘图显示计算结果。并将信号进行奇偶分解，分别得到奇对称信号h(n)-h(-n)与偶对称信号h(n)+h(-n)。用Excel 软件绘图显示计算结果。 二、实验程序代码 （1）离散抽样 double a,t; a=2*f*f*log(m); int i; for(i=0;i

数字信号处理上机实验代码

文件名：tstem.m（实验一、二需要） 程序： f unction tstem(xn,yn) %时域序列绘图函数 %xn:被绘图的信号数据序列，yn:绘图信号的纵坐标名称（字符串）n=0:length(xn)-1; stem(n,xn,'.'); xlabel('n');ylabel('yn'); axis([0,n(end),min(xn),1.2*max(xn)]); 文件名：tplot.m（实验一、四需要） 程序： function tplot(xn,T,yn) %时域序列连续曲线绘图函数 %xn:信号数据序列,yn:绘图信号的纵坐标名称（字符串） %T为采样间隔 n=0;length(xn)-1;t=n*T; plot(t,xn); xlabel('t/s');ylabel(yn); axis([0,t(end),min(xn),1.2*max(xn)]); 文件名：myplot.m（实验一、四需要）

%(1)myplot;计算时域离散系统损耗函数并绘制曲线图。function myplot(B,A) %B为系统函数分子多项式系数向量 %A为系统函数分母多项式系数向量 [H,W]=freqz(B,A,1000) m=abs(H); plot(W/pi,20*log10(m/max(m)));grid on; xlabel('\omega/\pi');ylabel('幅度(dB)') axis([0,1,-80,5]);title('损耗函数曲线'); 文件名：mstem.m（实验一、三需要） 程序： function mstem(Xk) %mstem(Xk)绘制频域采样序列向量Xk的幅频特性图 M=length(Xk); k=0:M-1;wk=2*k/M;%产生M点DFT对应的采样点频率（关于pi归一化值） stem(wk,abs(Xk),'.');box on;%绘制M点DFT的幅频特性图xlabel('w/\pi');ylabel('幅度'); axis([0,2,0,1.2*max(abs(Xk))]); 文件名：mpplot.m（实验一需要）

数字信号处理第二章上机题作业

数字信号处理作业实验题报告 第一章16.(1) 实验目的： 求解差分方程所描述的系统的单位脉冲响应和单位阶跃响应。 实验要求： 运用matlab求出y(n)=0.6y(n-1)-0.08y(n-2)+x(n)的单位脉冲响应和单位阶跃响应的示意图。 源程序： B1=1;A1=[1, -0.6, 0.08]; ys=2; %设差分方程 xn=[1, zeros(1, 20)]; %xn=单位脉冲序列，长度N=31 xi=filtic(B1, A1, ys); hn1=filter(B1, A1, xn, xi); %求系统输出信号hn1 n=0:length(hn1)-1; subplot(2, 1, 1);stem(n, hn1, '.') title('单位脉冲响应'); xlabel('n');ylabel('h(n)') xn=ones(1, 20); sn1=filter(B1, A1, xn, xi); %求系统输出信号sn1 n=0:length(sn1)-1; Subplot(2, 1, 2); stem(n, sn1, '.') title('单位阶跃响应'); xlabel('n'); ylabel('s(n)')

运行结果： 实验分析： 单位脉冲响应逐渐趋于0，阶跃响应保持不变，由此可见，是个稳定系统。

第二章31题 实验目的： 用matlab判断系统是否稳定。 实验要求： 用matlab画出系统的极，零点分布图，输入单位阶跃序列u(n)检查系统是否稳定。 源程序： A=[2, -2.98, 0.17, 2.3418, -1.5147]; B=[0, 0, 1, 5, -50]; subplot(2,1,1); zplane(B,A); %求H(z)的极点 p=roots(A); %求H(z)的模 pm=abs(p); if max(pm)<1 disp('系统因果稳定'), else,disp('系统因果不稳定'),end un=ones(1,800); sn=filter(B, A, un); n=0:length(sn)-1; subplot(2, 1, 2);plot(n, sn) xlabel('n');ylabel('s(n)')

数字信号处理基础实验报告 (2)

成都理工大学 《信号处理基础》实验 开设时间：2013—2014学年第2学期

题目1：信号的产生和显示 一、实验目的： 认识基本信号 通过使用MATLAB 设计简单程序, 掌握对MATLAB 的基本使用方法 二、实验原理： 找出下列表达式的信号与：正弦信号、最小相位信号、最大相位信号、零相位信号的对应关系。 1、sin60t 2、e-60t sin60t 3、(1- e-60t)sin60t 4、e60t sin60t 三、实验内容： 产生上述信号的信号并显示 （1）t=[-pi/30:0.001:pi/30]; f=sin(60*t); plot(t,f) 产生图形如下：

（2）t=[0:0.001:pi/30]; f=exp(-60*t).*sin(60*t); plot(t,f) 产生图形如下:

（3）t=[-5*pi/30:0.001:5*pi/30]; f=(1-exp(-60*t)).*sin(60*t); plot(t,f) 产生图形如下： (4) t=[-pi/30:0.001:pi/30]; f=exp(6*t).*sin(60*t); plot(t,f) 产生如下波形:

四、实验结果与讨论： 讨论上述信号的特点 从第一个波形图可以看出，它的波形与正弦函数sin（t）的相像，只是相位上有改变，是一个正弦信号。最大相位信号的能量集中在后面，最小相位能量集中在前面，所以第二个是一个最小相位，第四个是一个最大相位信号。第三个由于波形在t>0时没有，所以是一个零相位信号。 题目2：频谱分析与显示 一、实验目的 初步认识频谱分析

数字信号处理实验

子程序： function myplot(B,A) %myplot(B,A) %时域离散系统损耗函数绘图 %B为系统函数分子多项式系数向量 %A为系统函数分母多项式系数向量 [H,W]=freqz(B,A,1000); m=abs(H); plot(W/pi,20*log10(m/max(m)));grid on; xlabel('\omega/\pi');ylabel('幅度(dB)') axis([0,1,-80,5]);title('损耗函数曲线'); function tplot(xn,T,yn) %时域序列连续曲线绘图函数 % xn:信号数据序列，yn:绘图信号的纵坐标名称（字符串） % T为采样间隔 n=0:length(xn)-1;t=n*T; plot(t,xn); xlabel('t/s');ylabel(yn); axis([0,t(end),min(xn),1.2*max(xn)]) 程序： %实验4程序exp4.m % IIR数字滤波器设计及软件实现 clear all;close all Fs=10000;T=1/Fs; %采样频率 %调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st st=mstg; %低通滤波器设计与实现========================================= fp=280;fs=450; wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60; %DF指标（低通滤波器的通、阻带边界频）[N,wp]=ellipord(wp,ws,rp,rs); %调用ellipord计算椭圆DF阶数N和通带截止频率wp [B,A]=ellip(N,rp,rs,wp); %调用ellip计算椭圆带通DF系统函数系数向量B和A y1t=filter(B,A,st); %滤波器软件实现 % 低通滤波器设计与实现绘图部分 figure(2);subplot(3,1,1); myplot(B,A); %调用绘图函数myplot绘制损耗函数曲线 yt='y_1(t)'; subplot(3,1,2);tplot(y1t,T,yt); %调用绘图函数tplot绘制滤波器输出波形 %带通滤波器设计与实现==================================================== fpl=440;fpu=560;fsl=275;fsu=900; wp=[2*fpl/Fs,2*fpu/Fs];ws=[2*fsl/Fs,2*fsu/Fs];rp=0.1;rs=60;