数字信号处理MATLAB仿真.doc

实验一 数字信号处理的Matlab仿真一、实验目的1、掌握连续信号及其MATLAB实现方法；2、掌握离散信号及其MATLAB实现方法3、掌握离散信号的基本运算方法，以及MATLAB实现4、了解离散傅里叶变换的MATLAB实现5、了解IIR数字滤波器设计6、了解FIR数字滤波器设计1二、实验设备计算机，Matlab软件三、实验内容（一）、 连续信号及其MATLAB实现1、 单位冲击信号例1.1：t=1/A=50时，单位脉冲序列的MATLAB实现程序如下：clear all;t1=-0.5:0.001:0;A=50;A1=1/A;n1=length(t1);u1=zeros(1,n1);t2=0:0.001:A1;t0=0;u2=A*stepfun(t2,t0);t3=A1:0.001:1;n3=length(t3);u3=zeros(1,n3);t=t1 t2 t3;u=u1 u2 u3;plot(t,u)axis(-0.5 1 0 A+2)2、 任意函数例1.2：用MATLAB画出如下表达式的脉冲序列clear all;t=-2:1:3;N=length(t);x=zeros(1,N);x(1)=0.4;x(2)=0.8x(3)=1.2;x(4)=1.5;x(5)=1.0;x(6)=0.7;stem(t,x);axis(-2.2 3.2 0 1.7)3、 单位阶跃函数例1.3：用MATLAB实现单位阶跃函数clear all;t=-0.5:0.001:1;t0=0;u=stepfun(t,t0);plot(t,u)axis(-0.5 1 -0.2 1.2)4、 斜坡函数例1.4：用MATLAB实现g(t)=3(t-1)clear all;t=0:0.01:3;B=3;t0=1;u=stepfun(t,t0);n=length(t);for i=1:n u(i)=B*u(i)*(t(i)-t0);endplot(t,u)axis(-0.2 3.1 -0.2 6.2)5、 实指数函数例1.5：用MATLAB实现clear all;t=0:0.001:3;A=3;a=0.5;u=A*exp(a*t);plot(t,u)axis(-0.2 3.1 -0.2 14)6、 正弦函数例1.6：用MATLAB实现正弦函数f(t)=3cos(10t+1)clear all;t=-0.5:0.001:1;A=3;f=5;fai=1;u=A*sin(2*pi*f*t+fai);plot(t,u)axis(-0.5 1 -3.2 3.2)（二）、离散信号及其MATLAB实现1、 单位冲激序列例2.1：用MATLAB产生64点的单位冲激序列clear all;N=64;x=zeros(1,N);x(1)=1;xn=0:N-1;stem(xn,x)axis(-1 65 0 1.1)2、 任意序列例2.2：用MATLAB画出如下表达式的脉冲序列clear all;N=8;x=zeros(1,N);x(1)=8.0;x(2)=3.4x(3)=1.8;x(4)=5.6;x(5)=2.9;x(6)=0.7;xn=0:N-1;stem(xn,x)axis(-1 8 0 8.2)3、 单位阶跃序列例2.3：用MATLAB实现单位阶跃函数clear all;N=32;x=ones(1,N);xn=0:N-1;stem(xn,x)axis(-1 32 0 1.1)4、 斜坡序列例2.4：用MATLAB实现g(n)=3(n-4)点数为32的斜坡序列clear all;N=32;k=4B=3;t0=1;x=zeros(1,k) ones(1,N-k);for i=1:N x(i)=B*x(i)*(i-k);endxn=0:N-1;stem(xn,x)axis(-1 32 0 90)5、 正弦序列例2.5：用MATLAB实现幅度A=3，频率f=100，初始相位=1.2，点数为32的正弦信号clear all;N=32;A=3;f=100;fai=1.2;xn=0:N-1;x=A*sin(2*pi*f*(xn/N)+fai);stem(xn,x)axis(-1 32 -3.2 3.2)6、 实指数序列例2.6：用MATLAB实现，点数为32的实指数序列clear all;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的实指数序列clear all;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的均匀分布的随机序列与高斯随机序列clear all;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程序如下：clear all;N=32;w=100;k=3;x1=zeros(1,k);xn=0:N-1;x2=sin(100*xn);figure(1)stem(xn,x2)x=x1 x2;axis(-1 N -1.1 1.1)N=N+k;xn=0:N-1;figure(2)stem(xn,x)axis(-1 N -1.1 1.1)2、 信号相加若信号，值得注意的是当序列和的长度不相等或者位置不对应时，首先应该使两者的位置对齐，然后通过zeros函数左右补零使其长度相等后再相加例3.2：用MATLAB实现两序列相加clear all;n1=0:3x1=2 0.5 0.9 1;figure(1)stem(n1,x1)axis(-1 8 0 2.1 )n2=0:7x2= 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7;figure(2)stem(n2,x2)axis(-1 8 0 0.8 )n=0:7;x1=x1 zeros(1,8-length(n1);x2= zeros(1,8-length(n2),x2;x=x1+x2;figure(3)stem(n,x)axis(-1 8 0 2.1)3、 信号相乘信号序列和相乘所得信号的表达式为：这是样本与样本之间的点乘运算，在MATLAB中可采用“.*”来实现，但是在信号序列相乘之前，应对其做与相加运算一样的操作。例3.3：用MATLAB实现上例中两序列相乘clear all;n1=0:3x1=2 0.5 0.9 1;figure(1)stem(n1,x1)axis(-1 8 0 2.1 )n2=0:7x2= 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7;figure(2)stem(n2,x2)axis(-1 8 0 0.8 )n=0:7;x1=x1 zeros(1,8-length(n1);x2= zeros(1,8-length(n2),x2;x=x1.*x2;figure(3)stem(n,x)axis(-1 8 0 0.35)4、 信号翻转信号翻转的表达式为：y(n)=x(-n)，在MATLAB中可以用fliplr函数实现此操作例3.4：用MATLAB实现“信号相加”中的序列翻转clear all;n=0:3x1=2 0.5 0.9 1;x=fliplr(x1);stem(n,x)axis(-1 4 0 2.1 )5、 信号和对于N点信号，其和的定义为：例3.5：用MATLAB实现“信号相加”中的序列和clear all;n=0:3x1=2 0.5 0.9 1;x=sum(x1)6、 信号积对于N点信号，其积的定义为：例3.5：用MATLAB实现“信号相加”中的序列积clear all;n=0:3x1=2 0.5 0.9 1;x=prod(x1)（四）、离散傅里叶变换的MATLAB实现例：若是一个N=32的有限序列，利用MATLAB计算它的DFT并画出图形。N=32;n=0:N-1;xn=cos(pi*n/6);k=0:N-1;WN=exp(-j*2*pi/N);nk=n*k;WNnk=WN.nk;Xk=xn*WNnk;figure(1)stem(n,xn)figure(2)stem(k,abs(Xk)在MATLAB中，可以直接利用内部函数fft来实现FFT算法，该函数是机器语言，而不是MATLAB指令写成的，执行速度很快。常用格式为：y=fft(x)y=fft(x,N)（五）、IIR数字滤波器设计1、 基于巴特沃斯法直接设计IIR数字滤波器例.1：设计一个10阶的带通巴特沃斯数字滤波器，带通频率为100Hz到200Hz，采样频率为1000Hz，绘出该滤波器的幅频于相频特性，以及其冲击响应图clear all;N=10;Wn=100 200/500;b,a=butter(N,Wn,bandpass);freqz(b,a,128,1000)figure(2)y,t=impz(b,a,101);stem(t,y)2、 基于切比雪夫法直接设计IIR数字滤波器例5.2：设计一个切比雪夫型数字低通滤波器，要求：Ws=200Hz,Wp=100Hz,Rp=3dB,Rs=30dB,Fs=1000Hzclear all;Wp=100;Rp=3;Ws=200;Rs=30;Fs=1000;N,Wn=cheb1ord(Wp/(Fs/2),Ws/(Fs/2),Rp,Rs);b,a=cheby1(N,Rp,Wn);freqz(b,a,512,1000);例5.3：设计一个切比雪夫型数字带通滤波器，要求带通范围100-250Hz，带阻上限为300Hz，下限为50Hz，通带内纹波小于3dB，阻带纹波为30 dB，抽样频率为1000 Hz，并利用最小的阶次实现。clear all;Wpl=100;Wph=250;Wp=Wpl,Wph;Rp=3;Wsl=50;Wsh=300;Ws=Wsl,Wsh;Rs=30;Fs=1000;N,Wn=cheb2ord(Wp/(Fs/2),Ws/(Fs/2),Rp,Rs);b,a=cheby2(N,Rp,Wn);freqz(b,a,512,1000);（六）、FIR数字滤波器设计1、在MATLAB 中产生窗函数十分简单：（1）矩形窗（Rectangle Window）调用格式：w=boxcar(n)，根据长度n 产生一个矩形窗w。（2）三角窗（Triangular Window）调用格式：w=triang(n) ，根据长度n 产生一个三角窗w。（3）汉宁窗（Hanning Window）调用格式：w=hanning(n) ，根据长度n 产生一个汉宁窗w。（4）海明窗（Hamming Window）调用格式：w=hamming(n) ，根据长度n 产生一个海明窗w。（5）布拉克曼窗（Blackman Window）调用格式：w=blackman(n) ，根据长度n 产生一个布拉克曼窗w。（6）恺撒窗（Kaiser Window）调用格式：w=kaiser(n,beta) ，根据长度n 和影响窗函数旁瓣的参数产生一个恺撒窗w。2、基于窗函数的FIR 滤波器设计利用MATLAB 提供的函数firl 来实现调用格式：firl (n,Wn,ftype,Window)，n 为阶数、Wn 是截止频率（如果输入是形如W1 W2的矢量时，本函数将设计带通滤波器，其通带为W1=n0; %生成离散序列fstem(n,f,filled);title(单位阶跃序列)axis(n1 n2 0 1.1*max(f);%限定横坐标和纵坐标的显示范围ylabel(幅度f(n);n=0:16;f=exp(0.1+1.6*pi*j)*n); %建立信号fstem(n,f,fill),grid onxlabel(n),title(f(n);n=0:20; f=3*sin(pi/4*n); %建立信号fstem(n,f,fill),xlabel(n),grid ontitle(f(n)axis(0,20,-3.5,3.5); %限定横坐标和纵坐标的显示范围n=-20:20; f=sinc(n/5); %建立信号fstem(n,f,fill),xlabel(n),grid ontitle(f(n)(3)f=50;Um=2;nt=2; %显示周期个数N=16; %采样点数T=1/f; %信号周期dt=T/N; %采样时间间隔n=0:nt*N-1; %建立离散时间的时间序列tn=n*dt; %确定时间序列样点在时间轴上的位置x=1+0.5*Um*sawtooth(2*f*pi*tn);%产生时域信号subplot(2,1,1);plot(tn,x); %显示原连续信号axis(0 nt*T 0 2); %限定横坐标和纵坐标的显示范围ylabel(x(t);subplot(2,1,2);stem(tn,x); %显示经采样的信号axis(0 nt*T 0 2);ylabel(x(n);(4)f=200;nt=2; %显示周期个数Fs=4000;N=Fs/f; %求采样点数T=1/f; %求信号周期dt=T/N; %采样时间间隔n=0:nt*N-1; %建立离散时间的时间序列tn=n*dt; %确定时间序列样点在时间轴上的位置x=square(2*f*pi*tn); %产生时域信号subplot(2,1,1);plot(tn,x); %显示原连续信号axis(0 nt*T -1 1); %限定横坐标和纵坐标的显示范围ylabel(x(t);subplot(2,1,2);stem(tn,x); %显示经采样的信号axis(0 nt*T -1 1);ylabel(x(n);实验三:离散序列的基本运算（2）n1=-5;n2=5;n0=0;n3=-3;n4=4; n=n1:n2; %生成离散信号的时间序列f1=n=n0; %建立f1信号f2=n=n3; %建立f2信号f3=n=n4; %建立f3信号f=f2+2*f3; %建立f信号stem(n,f,filled); %画f信号的波形图axis(n1 n2 0 1.1*max(f); %限定横坐标和纵坐标的显示范围ylabel(f(n);n1=-5;n2=5;n0=0;n3=2;n4=-2;n=n1:n2;f1=n=n0; %建立f1信号f2=n=n3; %建立f2信号f3=n=n4; %建立f3信号f=f2+f3; %建立f信号stem(n,f,filled); %画f信号的波形图axis(n1 n2 0 1.1*max(f); %限定横坐标和纵坐标的显示范围ylabel(f(n);n1=0;n2=20;n3=3;n4=-3;n=n1:n2;x1=3*cos(2*pi*n/10); %建立原信号x(n）x2=3*cos(2*pi*(n-n3)/10); %建立x(n-3)信号x3=3*cos(2*pi*(n-n4)/10); %建立x(n+3)信号subplot(3,1,1);stem(n,x1,filled);%画x1的波形图ylabel(x(n);subplot(3,1,2);stem(n,x2,filled);%画x(n-3)的波形图ylabel(x(n-3);subplot(3,1,3);stem(n,x3,filled);%画x(n+3)的波形图ylabel(x(n+3);n1=0;n2=24;n=n1:n2;x1=exp(-n)/16); %建立x1(n)信号x2=5*sin(2*pi*n/10); %建立x2(n)信号x=x1.*x2; %建立x1(n)*x2(n)信号subplot(3,1,1);stem(n,x1,filled);%画x1(n)的波形图ylabel(x1(n);subplot(3,1,2);stem(n,x2,filled);%画x2(n)的波形图ylabel(x1(n);subplot(3,1,3);stem(n,x,filled); %画x1(n)*x2(n)的波形图ylabel(x(n);（3）n1=0;n2=20;n=n1:n2;x0=n.*sin(n); %建立原信号x(n）y1=(n-3).*sin(n-3); %建立信号x(n-3）y2=(-n).*sin(-n); %建立信号x(-n）y3=-(n.*sin(n); %建立信号-x(n）y4=(-n+3).*sin(-n+3); %建立信号x(-n+3）y5=n/2.*sin(n/2); %建立信号x(n/2）subplot(5,1,1);stem(n,y1,filled);ylabel(y1(n)=x(n-3);subplot(3,2,2);stem(n,y2,filled);ylabel(y2(n)=x(-n);subplot(3,2,3);stem(n,y3,filled);ylabel(y3(n)=-x(n);subplot(5,1,4);stem(n,y4,filled);ylabel(y4(n)=x(-n+3);subplot(3,2,5);stem(n,y5,filled);ylabel(y5(n)=x(n/2);（4）n1=-4:-1;x1=2*n1+5; %建立信号x1(n）n2=0:4;x2=6; %建立信号x2(n）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=min(n1)&(n=min(n2)&(n=max(n2)=x2; %为y2赋值x=y1+y2; %建立信号xy3=2*x %建立信号y3y4=2*x; %建立信号y4y5=x; %建立信号y5subplot(3,1,1);stem(n3,y3,filled);title(x1(n);axis(min(n3),max(n3),1.1*min(y3),1.1*max(y3);subplot(3,1,2);stem(n4,y4,filled);title(x2(n);axis(min(n4),max(n4),1.1*min(y4),1.1*max(y4);subplot(3,1,3);stem(n5,y5,filled);title(x3(n);axis(min(n5),max(n5),1.1*min(y5),1.1*max(y5);附录一 信号、系统和系统响应 1、理想采样信号序列 （1）首先产生信号x(n),0=n=50 n=0:50; %定义序列的长度是50 A=444.128; %设置信号有关的参数 a=50*sqrt(2.0)*pi; T=0.001; %采样率 w0=50*sqrt(2.0)*pi; x=A*exp(-a*n*T).*sin(w0*n*T); %pi是MATLAB定义的，信号乘可采用“.*” close all %清除已经绘制的x(n)图形 subplot(3,1,1);stem(x); %绘制x(n)的图形 title(理想采样信号序列); （2）绘制信号x(n)的幅度谱和相位谱 k=-25:25; W=(pi/12.5)