1、数字信号处理指导书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）表示，其幅度可以在某一范围内连
2、续取值。由于信号处理所用的设备主要是计算机或专用的信号处理芯片，均以有限的位数来表示信号的幅度，因此，信号的幅度也必须“量化”，即取离散值。我们把时间和幅度上均取离散值的信号称为时域离散信号或数字信号。在MATLAB中，时域离散信号可以通过编写程序直接生成，也可以通过对连续信号的等间隔抽样获得。 下面介绍常用的时域离散信号及其程序。1、单位抽样序列MATLAB源程序为1) function x,n = impuls (n0,n1,n2)% Generates x(n) = delta(n-n0); n=n0 处建立一个单位抽样序列% x,n = impuls (n0,n1,n2)if (n0
3、n2) | (n1 n2)error(arguments must satisfy n1 = n0 =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); %对所有样点置0stem(n,x,fil
4、led,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; %时间序列样点在时
5、间轴上的位置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个周期的信号波形，并对信号的
6、一个周期进行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; %产生时域信号,且幅度在02V之间 subplot(2,1,1);stairs(tn,x,k); axis(0 nt*T 1.1*min(x) 1.1*max(x);ylabel(x(t);subplot(2,1,2);stem(tn,x,filled,k); axis(0 nt*T 1.1*min(x) 1.1*max(x);ylabel(x(
7、n);注意：直接用square子函数产生的信号波形，其幅度为-11之间。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）= （显示-3n4）（2）f（n）=u（n-1） （显示-5n=min(n1)）&(n=max
8、(n1)）；在符合关系运算条件的范围内寻找非零元素的索引号。2、fliplr功能：对矩阵行元素进行左右翻转。调用格式：x1= fliplr（x）；将x的行元素左右翻转，赋给变量x1。三、实验原理离散序列的时域运算包括信号的相加、相乘，信号的时域变换包括信号的移位、反折、尺度变换等。在MATLAB中，离散序列的相加、相乘等运算是两个向量之间的运算，因此参加运算的两个序列向量必须具有相同的维数，否则应进行相应的处理。下面介绍各种离散序列的时域运算和时域变换的性质。1、 序列的移位x1（n）=x（n-m）例、x0=u（n），（显示-10n10）x1=u（n+6）（显示-10n10）,x2= u（n-
9、4）（显示-10n=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*mi
10、n(x2),1.1*max(x2);ylabel(u(n-4);例、已知x（n）=，求x（n-2）和x（n+2）在-210区间的波形。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);su
11、bplot(3,1,3),stem(n,x2,filled,k); ylabel(x(n+2);2、 序列的相加情况1 两序列具有相同的维数例、求x(n)= (0n10)MATLAB源程序为：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
12、),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） （-4n6）X2(n)=u(n-4) (-5n=0; %建立x1信号n2=-5:8;n02=4;x2=(n2-n02)=0; %建立x2信号n=min(n1,n2):max(n1,n2); %为x信号建立时间序列nN=length(n); %求时间序列n的点数Ny1=zeros(1,N); %新建一维N列的y1全
13、0数组y2=zeros(1,N); %新建一维N列的y2全0数组y1(find(n=min(n1)&(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
14、,k);ylabel(x(n);axis(min(n),max(n),1.1*min(x),1.1*max(x);3、 序列相乘例、x1（n）= （-4n10）X2(n)=u(n+1) (-2n=0; %建立x2信号n=min(n1,n2):max(n1,n2); %为x信号建立时间序列nN=length(n); %求时间序列n的点数Ny1=zeros(1,N); %新建一维N列的y1全0数组y2=zeros(1,N); %新建一维N列的y2全0数组y1(find(n=min(n1)&(n=min(n2)&(n=max(n2)=x2; %为y2赋值x=y1.*y2;subplot(3,1,1)
15、,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）= （-4n4）求x（-n）。MATLAB源程序为：n=-4
16、: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)信
17、号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）= （-5n5）2、x（n）=u（n-2）+u（n+2） （-5na=1 2 ；P=po
18、ly(a)；则P = 1 -3 210、poly2strPa= poly2str(a,s)例如： P= 1,-3,2; PA=poly2str(P,s)则PA = s2 - 3 s + 211、pzmap功能：显示连续系统的零极点分布图。调用格式：pzmap（b，a）；绘制由行向量b和a构成的系统函数确定的零极点分布图pzmap(p，z)；绘制由列向量z确定的零点和p确定的极点构成的零极点分布图。三、实验原理由于IIR数字滤波器是在已知的归一化的低通模拟滤波器的基础上设计的，主要包括巴特沃斯、切比雪夫、椭圆低通滤波器，因此把这些低通滤波器称为滤波器原型。下面介绍各种滤波器的设计。例1、 进行巴
19、特沃斯滤波器原型的设计，获得任意阶数N的归一化的系统函数公式和幅频响应。MATLAB源程序为：%巴特沃斯模拟滤波器原型n=input(N= ); %输入滤波器阶数N%计算n阶模拟低通原型，得到左半平面零极点z0,p0,k0=buttap(n); b0=k0*poly(z0) %求滤波器系数b0a0=poly(p0) %求滤波器系数a0h,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=po
20、ly2str(a0,p) %给出a0决定的关于p多项式输入阶数N=2则运行结果为b0 = 1a0 = 1.0000 1.4142 1.0000Pb = 1Pa = p2 + 1.4142 p + 1即归一化的2阶巴特沃斯滤波器原型的系统函数为例2、 通过模拟滤波器原型设计一个巴特沃斯模拟低通滤波器，要求通带截止频率2khz，通带最大衰减1db，阻带截止频率5khz，阻带最小衰减20db。MATLAB源程序为：%巴特沃斯模拟滤波器fp=2000;wp=2*pi*fp %输入滤波器的通带截止频率fs=5000;ws=2*pi*fs %输入滤波器的阻带截止频率Rp=1;As=20; %输入滤波器的通
21、阻带衰减指标%计算滤波器的阶数和3dB截止频率N,wc=buttord(wp,ws,Rp,As,s) %计算n阶模拟低通原型z0,p0,k0=buttap(n) %字母后加0表明这是原型滤波器的各指标，而不是所求的滤波器的b0=k0*poly(z0) %求归一化滤波器分子系数b0a0=poly(p0) %求归一化滤波器分母系数a0H,w=freqs(b0,a0); %求归一化系统的频率特性dbH=20*log10(abs(H)/max(abs(H); %将归一化系统的幅频特性化为分贝值，注意此时的分贝值为负的，为了使分贝图和幅频响应图一致subplot(2,2,1),plot(w*wc/(2*
22、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); %画所求滤
