数字信号处理指导书matlab版.docx
《数字信号处理指导书matlab版.docx》由会员分享，可在线阅读，更多相关《数字信号处理指导书matlab版.docx（46页珍藏版）》请在冰点文库上搜索。
数字信号处理指导书matlab版
实验1时域离散信号的产生
一、实验目的
●学会运用MATLAB产生常用离散时间信号。
二、实验涉及的matlab子函数
1、square
功能：
产生矩形波
调用格式：
x=square（t）；类似于sin（t），产生周期为2*pi，幅值为+—1的方波。
x=square（t，duty）；产生制定周期的矩形波，其中duty用于指定脉冲宽度与整个周期的比例。
2、rand
功能：
产生rand随机信号。
调用格式：
x=rand（n，m）；用于产生一组具有n行m列的随机信号。
三、实验原理
在时间轴的离散点上取值的信号，称为离散时间信号。
通常，离散时间信号用x（n）表示，其幅度可以在某一范围内连续取值。
由于信号处理所用的设备主要是计算机或专用的信号处理芯片，均以有限的位数来表示信号的幅度，因此，信号的幅度也必须“量化”，即取离散值。
我们把时间和幅度上均取离散值的信号称为时域离散信号或数字信号。
在MATLAB中，时域离散信号可以通过编写程序直接生成，也可以通过对连续信号的等间隔抽样获得。
下面介绍常用的时域离散信号及其程序。
1、单位抽样序列
MATLAB源程序为
1）function[x,n]=impuls（n0,n1,n2）
%Generatesx（n）=delta（n-n0）;n=n0处建立一个单位抽样序列
%[x,n]=impuls（n0,n1,n2）
if（（n0n2）|（n1>n2））
error（'argumentsmustsatisfyn1<=n0<=n2'）
end
n=[n1:
n2];
x=[zeros（1,（n0-n1））,1,zeros（1,（n2-n0））];
将上述文件存为：
impuls.m，在命令窗口输入
n0=0,n1=-10,n2=11;
[x,n]=impuls（n0,n1,n2）;stem（n,x,’filled’）
2）n1=-5;n2=5;n0=0;
n=n1:
n2;
x=[n==n0];
stem（n,x,'filled','k'）;
axis（[n1,n2,1.1*min（x）,1.1*max（x）]）;
title（'单位脉冲序列'）;
xlabel（'时间（n）'）;
ylabel（'幅度x（n）'）;
3）n1=-5;n2=5;k=0;
n=n1:
n2;
nt=length（n）;%求n点的个数
nk=abs（k-n1）+1;%确定k在n序列中的位置
x=zeros（1,nt）;%对所有样点置0
x（nk）=1;%对抽样点置1
stem（n,x,'filled','k'）;
axis（[n1,n2,0,1.1*max（x）]）;
title（'单位脉冲序列'）;
xlabel（'时间（n）'）;
Ylabel（'幅度x（n）'）;
2、单位阶跃序列
MATLAB源程序为：
1）n1=-2;n2=8;n0=0;
n=n1:
n2;%生成离散信号的时间序列
x=[n>=n0];%生成离散信号x（n）
stem（n,x,'filled','k'）;%绘制脉冲杆图，且圆点处用实芯圆表示
axis（[n1,n2,0,1.1*max（x）]）;
title（'单位阶跃序列'）;
xlabel（'时间（n）'）;
Ylabel（'幅度x（n）'）;
2）n1=-2;n2=8;k=0;
n=n1:
n2;
nt=length（n）;%求n点的个数
nk=abs（k-n1）+1;%确定k在n序列中的位置
x=[zeros（1,nk-1）,ones（1,nt-nk+1）];%对所有样点置0
stem（n,x,'filled','k'）;
axis（[n1,n2,0,1.1*max（x）]）;
title（'单位阶跃序列'）;
xlabel（'时间（n）'）;
ylabel（'幅度x（n）'）;
3、正弦序列
x（n）=Um
例、已知一时域周期性正弦信号的频率为1HZ，振幅幅度为1V，在窗口中显示两个周期的信号波形，并对该信号的一个周期进行32点采样获得离散信号。
显示连续信号和采样获得离散信号波形。
MATLAB源程序为：
f=1;Um=1;nt=2;;%输入信号频率、振幅和显示周期数
N=32;T=1/f;%N为采样点数，T为窗口显示时间
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）'）;
4、矩形序列
将square表示式中的t换成n，且n取整数，则可以获得矩形序列。
例、一个连续的周期性矩形信号频率为5kHZ，信号幅度为0-2V之间，脉冲宽度与周期的比例为1:
4,且要求在窗口上显示其2个周期的信号波形，并对信号的一个周期进行16点采样来获得离散信号，显示原连续信号与采样获得的离散信号。
MATLAB源程序为：
f=5000;nt=2;
N=16;T=1/f;
dt=T/N;
n=0:
nt*N-1;
tn=n*dt;
x=square（2*f*pi*tn,25）+1;%产生时域信号,且幅度在0～2V之间
subplot（2,1,1）;stairs（tn,x,'k'）;
axis（[0nt*T1.1*min（x）1.1*max（x）]）;
ylabel（'x（t）'）;
subplot（2,1,2）;stem（tn,x,'filled','k'）;
axis（[0nt*T1.1*min（x）1.1*max（x）]）;
ylabel（'x（n）'）;
注意：
直接用square子函数产生的信号波形，其幅度为-1~1之间。
5、rand函数
在实际系统的研究和处理中，常常需要产生随机信号，MATLAB提供的rand函数可以为我们生成随机信号。
例、生成一组41点构成的连续随机信号和与之相应的随机序列。
MATLAB源程序为：
tn=0:
40;
N=length（tn）;
x=rand（1,N）;
subplot（1,2,1）,plot（tn,x,'k'）;
subplot（1,2,2）,stem（tn,x,'filled','k'）;
四、实验任务
1、产生离散序列：
（1）f（n）=
（显示-3（2）f（n）=u（n-1）（显示-52、一个连续的周期性正弦信号频率为50HZ，信号幅度在0~2V之间，在窗口上显示2个周期的信号波形，对信号的一个周期进行16点采样来获得离散信号，显示原连续信号和采样获得的离散信号波形。
3、一个连续的周期性方波信号频率为200HZ,信号幅度在0~2V之间，在窗口上显示两个周期的信号波形，用4Khz的频率对连续信号进行采样，显示原连续信号和采样获得的离散信号波形。
实验2离散序列的基本运算
一、实验目的
●学会运用MATLAB进行离散序列的运算，并掌握程序的编写方法。
二、实验涉及的matlab子函数
1、find
功能：
寻找非零元素的索引号。
调用格式：
find（（n>=min（n1））&（n<=max（n1）））；在符合关系运算条件的范围内寻找非零元素的索引号。
2、fliplr
功能：
对矩阵行元素进行左右翻转。
调用格式：
x1=fliplr（x）；将x的行元素左右翻转，赋给变量x1。
三、实验原理
离散序列的时域运算包括信号的相加、相乘，信号的时域变换包括信号的移位、反折、尺度变换等。
在MATLAB中，离散序列的相加、相乘等运算是两个向量之间的运算，因此参加运算的两个序列向量必须具有相同的维数，否则应进行相应的处理。
下面介绍各种离散序列的时域运算和时域变换的性质。
1、序列的移位
x1（n）=x（n-m）
例、x0=u（n），（显示-10x1=u（n+6）（显示-10x2=u（n-4）（显示-10MATLAB源程序为：
n1=-10;n2=10;k0=0;k1=-6;k2=4;
n=n1:
n2;%生成离散信号的时间序列
x0=[n>=k0];%生成离散信号x0（n）
x1=[（n-k1）>=0];%生成离散信号x1（n）
x2=[（n-k2）>=0];%生成离散信号x2（n）
subplot（3,1,1）,stem（n,x0,'filled','k'）;
axis（[n1,n2,1.1*min（x0）,1.1*max（x0）]）;
ylabel（'u（n）'）;
subplot（3,1,2）,stem（n,x1,'filled','k'）;
axis（[n1,n2,1.1*min（x1）,1.1*max（x1）]）;
ylabel（'u（n+6）'）;
subplot（3,1,3）,stem（n,x2,'filled','k'）;
axis（[n1,n2,1.1*min（x2）,1.1*max（x2）]）;
ylabel（'u（n-4）'）;
例、已知x（n）=
，求x（n-2）和x（n+2）在-2~10区间的波形。
MATLAB源程序为：
n=-2:
10;n0=2;n1=-2;
x=2*sin（2*pi*n/10）;%建立原信号x（n）
x1=2*sin（2*pi*（n-n0）/10）;%建立x（n-2）信号
x2=2*sin（2*pi*（n-n1）/10）;%建立x（n+2）信号
subplot（3,1,1）,stem（n,x,'filled','k'）;
ylabel（'x（n）'）;
subplot（3,1,2）,stem（n,x1,'filled','k'）;
ylabel（'x（n-2）'）;
subplot（3,1,3）,stem（n,x2,'filled','k'）;
ylabel（'x（n+2）'）;
2、序列的相加
情况1两序列具有相同的维数
例、求x（n）=
（0MATLAB源程序为：
n1=0;n2=10;n01=2;n02=4;
n=n1:
n2;
x1=[（n-n01）==0];
x2=[（n-n02）==0];
x3=x1+x2;
subplot（3,1,1）;stem（n,x1,'filled','k'）;
axis（[n1,n2,1.1*min（x1）,1.1*max（x1）]）;
ylabel（'δ（n-2）'）;
subplot（3,1,2）;stem（n,x2,'filled','k'）;
axis（[n1,n2,1.1*min（x2）,1.1*max（x2）]）;
ylabel（'δ（n-4）'）;
subplot（3,1,3）;stem（n,x3,'filled','k'）;
axis（[n1,n2,1.1*min（x3）,1.1*max（x3）]）;
ylabel（'δ（n-2）+δ（n-4）'）;
情况2两序列具有不同的维数
例、x1（n）=u（n+2）（-4X2（n）=u（n-4）（-5求x（n）=x1（n）+X2（n）=u（n-4）
MATLAB源程序为：
n1=-4:
6;n01=-2;
x1=[（n1-n01）>=0];%建立x1信号
n2=-5:
8;n02=4;
x2=[（n2-n02）>=0];%建立x2信号
n=min（[n1,n2]）:
max（[n1,n2]）;%为x信号建立时间序列n
N=length（n）;%求时间序列n的点数N
y1=zeros（1,N）;%新建一维N列的y1全0数组
y2=zeros（1,N）;%新建一维N列的y2全0数组
y1（find（（n>=min（n1））&（n<=max（n1））））=x1;%为y1赋值
y2（find（（n>=min（n2））&（n<=max（n2））））=x2;%为y2赋值
x=y1+y2;
subplot（3,1,1）,stem（n1,x1,'filled','k'）;
ylabel（'x1（n）'）;
axis（[min（n）,max（n）,1.1*min（x）,1.1*max（x）]）;
subplot（3,1,2）,stem（n2,x2,'filled','k'）;
ylabel（'x2（n）'）;
axis（[min（n）,max（n）,1.1*min（x）,1.1*max（x）]）;
subplot（3,1,3）,stem（n,x,'filled','k'）;
ylabel（'x（n）'）;
axis（[min（n）,max（n）,1.1*min（x）,1.1*max（x）]）;
3、序列相乘
例、x1（n）=
（-4X2（n）=u（n+1）（-2求x（n）=x1（n）xX2（n）
MATLAB源程序为：
n1=-4:
10;
x1=3*exp（-0.25*n）;%建立x1信号
n2=-2:
6;n02=-1;
x2=[（n2-n02）>=0];%建立x2信号
n=min（[n1,n2]）:
max（[n1,n2]）;%为x信号建立时间序列n
N=length（n）;%求时间序列n的点数N
y1=zeros（1,N）;%新建一维N列的y1全0数组
y2=zeros（1,N）;%新建一维N列的y2全0数组
y1（find（（n>=min（n1））&（n<=max（n1））））=x1;%为y1赋值
y2（find（（n>=min（n2））&（n<=max（n2））））=x2;%为y2赋值
x=y1.*y2;
subplot（3,1,1）,stem（n1,x1,'filled','k'）;
ylabel（'x1（n）'）;
axis（[min（n）,max（n）,1.1*min（x1）,1.1*max（x1）]）;
subplot（3,1,2）,stem（n2,x2,'filled','k'）;
ylabel（'x2（n）'）;
axis（[min（n）,max（n）,1.1*min（x2）,1.1*max（x2）]）;
subplot（3,1,3）,stem（n,x,'filled','k'）;
ylabel（'x（n）'）;
axis（[min（n）,max（n）,1.1*min（x）,1.1*max（x）]）;
4、序列反折
例、x（n）=
（-4求x（-n）。
MATLAB源程序为：
n=-4:
4;
x=exp（-0.3*n）;
x1=fliplr（x）;
n1=-fliplr（n）;
subplot（1,2,1）,stem（n,x,'filled','k'）;
title（'x（n）'）;
subplot（1,2,2）,stem（n1,x1,'filled','k'）;
title（'x（-n）'）;
5、序列的尺度变换
例、x（n）=sin（2*pi*n），求x（2n）和x（n/2）
MATLAB源程序为：
n=0:
40;tn=n./20;%每周期取20个点
x=sin（2*pi*tn）;%建立原信号x（n）
x1=sin（2*pi*tn*2）;%建立x（2n）信号
x2=sin（2*pi*tn/2）;%建立x（n/2）信号
subplot（3,1,1）,stem（tn,x,'filled','k'）;
ylabel（'x（n）'）;
axis（[0,2,1.1*min（x）,1.1*max（x）]）;
subplot（3,1,2）,stem（tn,x1,'filled','k'）;
ylabel（'x（2n）'）;
axis（[0,2,1.1*min（x）,1.1*max（x）]）;
subplot（3,1,3）,stem（tn,x2,'filled','k'）;
ylabel（'x（n/2）'）;
axis（[0,2,1.1*min（x）,1.1*max（x）]）;
四、实验任务
1、x（n）=
（-52、x（n）=u（n-2）+u（n+2）（-53、已知x（n）=3cos（2*pi*n/10）,显示x（n-3）和x（n+3）在0~20区间的波形。
4、已知x1=exp（-n/16），x2=5sin（2*pi*n/10）,显示x1乘以x2在0~24区间的波形。
5、已知x（n）=nsin（n），显示0~20区间的波形
y1（n）=x（n-3）,y2（n）=x（-n），y3=-x（n）,y4=x（-n+3）,y5（n）=x（n/2）.
实验3模拟原型滤波器的设计
一、实验目的
●学会运用MATLAB设计模拟低通滤波器原型的设计方法。
二、实验涉及的matlab子函数
1、buttord
功能：
确定巴特沃斯滤波器的阶数和3dB截止频率。
调用格式：
[N,Wc]=buttord（wp,ws,Rp,As,’s’），计算巴特沃斯滤波器的阶数和3dB截止频率。
Rp为通带最大衰减，As为阻带最小衰减。
2、che1ord
功能：
确定切比雪夫1型滤波器的阶数和通带截止频率。
调用格式：
[N,Wp]=buttord（wp,ws,Rp,As,’s’）
3、che2ord
功能：
确定切比雪夫2型滤波器的阶数和阻带截止频率。
调用格式：
[N,Wn]=buttord（wp,ws,Rp,As,’s’）
4、ellipord
功能：
确定椭圆滤波器的阶数和通带截止频率。
调用格式：
[N,Wn]=buttord（wp,ws,Rp,As,’s’）。
5、buttap
功能：
巴特沃斯模拟低通滤波器原型（即归一化的滤波器）。
调用格式：
[z，p，k]=buttap（n），设计巴特沃斯模拟低通滤波器原型，其系统函数为
6、cheb1ap
功能：
切比雪夫1型模拟低通滤波器原型。
调用格式：
[z，p，k]=cheb1ap（n，Rp），设计切比雪夫1型模拟低通滤波器原型，其通带最大衰减为Rp，系统函数为
7、cheb2ap
功能：
切比雪夫2型模拟低通滤波器原型。
调用格式：
[z，p，k]=cheb2ap（n，As），设计切比雪夫2型模拟低通滤波器原型，其阻带最小衰减为As，系统函数为
8、ellipap
功能：
椭圆模拟低通滤波器原型。
调用格式：
[z，p，k]=ellipap（n，Rp，As），设计椭圆模拟低通滤波器原型，其系统函数为
9、poly
功能：
求某向量指定根所对应的特征多项式。
调用格式：
P=poly（
），多项式P是一个特征多项式，
的元素是多项式P的根。
例如：
>>a=[12]；P=poly（a）；
则P=1-32
10、poly2str
Pa=poly2str（a,’s’）
例如：
>>P=[1,-3,2];PA=poly2str（P,'s'）
则PA=s^2-3s+2
11、pzmap
功能：
显示连续系统的零极点分布图。
调用格式：
pzmap（b，a）；绘制由行向量b和a构成的系统函数确定的零极点分布图
pzmap（p，z）；绘制由列向量z确定的零点和p确定的极点构成的零极点分布图。
三、实验原理
由于IIR数字滤波器是在已知的归一化的低通模拟滤波器的基础上设计的，主要包括巴特沃斯、切比雪夫、椭圆低通滤波器，因此把这些低通滤波器称为滤波器原型。
下面介绍各种滤波器的设计。
例1、进行巴特沃斯滤波器原型的设计，获得任意阶数N的归一化的系统函数公式和幅频响应。
MATLAB源程序为：
%巴特沃斯模拟滤波器原型
n=input（'N='）;%输入滤波器阶数N
%计算n阶模拟低通原型，得到左半平面零极点
[z0,p0,k0]=buttap（n）;
b0=k0*poly（z0）%求滤波器系数b0
a0=poly（p0）%求滤波器系数a0
[h,w]=freqs（b0,a0）;%显示系统的频率特性
plot（w,abs（h）,'r'）,%画幅频特性图
axis（[0,5,0,1.1]）,
ylable（'幅度'）;xlabel（'f（HZ）'）;
Pb=poly2str（b0,'p'）%给出b0决定的关于p多项式
Pa=poly2str（a0,'p'）%给出a0决定的关于p多项式
输入阶数N=2
则运行结果为
b0=1
a0=1.00001.41421.0000
Pb=1
Pa=p^2+1.4142p+1
即归一化的2阶巴特沃斯滤波器原型的系统函数为
例2、通过模拟滤波器原型设计一个巴特沃斯模拟低通滤波器，要求通带截止频率2khz，通带最大衰减1db，阻带截止频率5khz，阻带最小衰减20db。
MATLAB源程序为：
%巴特沃斯模拟滤波器
fp=2000;wp=2*pi*fp%输入滤波器的通带截止频率
fs=5000;ws=2*pi*fs%输入滤波器的阻带截止频率
Rp=1;As=20;%输入滤波器的通阻带衰减指标
%计算滤波器的阶数和3dB截止频率
[N,wc]=buttord（wp,ws,Rp,As,'s'）
%计算n阶模拟低通原型
[z0,p0,k0]=buttap（n）%字母后加0表明这是原型滤波器的各指标，而不是所求的滤波器的
b0=k0*poly（z0）%求归一化滤波器分子系数b0
a0=poly（p0）%求归一化滤波器分母系数a0
[H,w]=freqs（b0,a0）;%求归一化系统的频率特性
dbH=20*log10（abs（H）/max（abs（H）））;%将归一化系统的幅频特性化为分贝值，注意此时的分贝值为负的，为了使分贝图和幅频响应图一致
subplot（2,2,1）,plot（w*wc/（2*pi）,abs（H））,grid%画所求滤波器的幅频响应图
axis（[0,6000,0,1.1]）;ylabel（'幅度'）;xlabel（'f（Hz）'）;
subplot（2,2,2）,plot（w*wc/（2*pi）,angle（H））,grid%画所求滤波器的相频响应图
axis（[0,6000,-4,4]）;ylabel（'相位'）;xlabel（'f（Hz）'）;
subplot（2,2,3）,plot（w*wc/（2*pi）,dbH）,grid%画所求滤波器的幅频响应分贝图
axis（[0,6000,-30,2]）;
ylabel（'幅度（dB）'）;xlabel（'f（Hz）'）;
subplot（2,2,4）,plot（p0*wc,'x'）;%画所求滤
升级会员 