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

1、 实验一 常见离散信号的MATLAB产生和图形显示 实验目的：加深对常用离散信号的理解；实验原理：单位抽样序列 在MATLAB中可以利用zeros()函数实现。如果在时间轴上延迟了k个单位，得到即： 2单位阶跃序列 在MATLAB中可以利用ones()函数实现。 3正弦序列在MATLAB中4复正弦序列在MATLAB中5指数序列在MATLAB中实验内容：编制程序产生上述5种信号（长度可输入确定），并绘出其图形。实验要求：讨论复指数序列的性质。实验过程：单位冲击序列： n=0:10; x1=1 zeros(1,10); x2=zeros(1,8) 1 zeros(1,8); subplot(1,2

2、,1); stem(n,x1); xlabel (时间序列n); ylabel(幅度); title(单位冲激序列（n）); subplot(1,2,2); stem(x2); xlabel(时间序列n); ylabel(幅度); title(延时了8个单位的冲激序列（n-8）); 单位阶跃序列： n=0:10; u=ones(1,11); stem(n,u); xlabel (时间序列n); ylabel(信号幅度); title(单位阶跃序列u（n）);正弦序列： n=1:30; x=2*sin(pi*n/6+pi/4); stem(n,x); xlabel (时间序列n); ylabel

3、(振幅); title(正弦函数序列x=2*sin(pi*n/6+pi/4);复指数序列： n=1:30; x=5*exp(j*3*n); stem(n,x); xlabel (时间序列n); ylabel(振幅); title(复指数序列x=5*exp(j*3*n);指数序列： n=1:30; x=1.8.n; stem(n,x); xlabel (时间序列n); ylabel(振幅); title(指数序列x=1.8.n);复指数序列的周期性讨论： 为了研究复指数序列的周期性质，我们分别作了正弦函数x1=1.5sin(0.3n)和x2=sin(0.6n); 的幅度特性图像。由下图看出：x1

4、=1.5sin(0.3n)的周期是20，而x2=sin(0.6n)是非周期的。理论计算中对第一个，N2/（0.3）20/3，为有理数。而第二个0.6不是的倍数，所以不是周期的。因此可以看出，实验结果和理论相符。对于离散复指数函数，只有当是纯虚数，且纯虚数的系数是的倍数时，才是周期的。其它情况下均不是。MATLAB实现过程： n=0:80; x1=1.5*sin(0.3*pi*n); x2=sin(0.6*n); subplot(1,2,1); stem(n,x1); xlabel (时间序列n); ylabel(振幅); title(正弦序列x1=1.5*sin(0.3*pi*n); subp

5、lot(1,2,2); stem(n,x2); xlabel (时间序列n); ylabel(振幅); title(正弦序列x1=1.5*sin(0.3*pi*n);实验2 离散系统的差分方程、冲激响应和卷积分析实验目的：加深对离散系统的差分方程、冲激响应和卷积分析方法的理解。实验原理：离散系统其输入、输出关系可用以下差分方程描述： 输入信号分解为冲激信号，。记系统单位冲激响应，则系统响应为如下的卷积计算式： 当时，hn是有限长度的（n：0，M），称系统为FIR系统；反之，称系统为IIR系统。 在MATLAB中，可以用函数y=Filter(p,d,x) 求解差分方程，也可以用函数 y=Conv

6、(x,h)计算卷积。实验内容：编制程序求解下列两个系统的单位冲激响应和阶跃响应，并绘出其图形。实验要求：给出理论计算结果和程序计算结果并讨论。实验过程：（1）单位冲激响应： a=1,0.6,0.08; b=1,-1; N=20;n=0:1:N; x=1,zeros(1,N); y=filter(b,a,x); stem(n,y); xlabel(时间序列n); ylabel(信号幅度); title(单位冲激响应h(n);（2）单位阶跃响应： a=1,0.6,0.08; b=1,-1; N=20;n=0:1:N-1; x= ones(1,N); y=filter(b,a,x); stem(n,

7、y); xlabel(时间序号); ylabel(信号幅度); title(单位阶跃响应h（n）);理论分析：由差分方程得系统函数为：利用分部分式法可得：，z反变换得： h(n)即为单位冲击响应。对于阶跃响应， ,所以阶跃响应 = = 所以，可见，实验结果与理论分析是一致的。单位冲激响应 a=1; b=0.2,0.2,0.2,0.2,0.2,0.2; N=20;n=0:1:N; x=1,zeros(1,N); y=filter(b,a,x); stem(n,y); xlabel(时间序号n); ylabel(信号幅度); title(单位冲激响应h（n）);单位阶跃响应 a=1; b=0.2,

8、0.2,0.2,0.2,0.2,0.2; N=20;n=0:1:N-1; x=ones(1,N); y=filter(b,a,x); stem(n,y); xlabel(时间序号); ylabel(信号幅度); title(单位阶跃响应h（n）);理论分析：由差分方程得系统函数为：z反变换得： h(n)即为单位冲击响应。对于阶跃响应， ,所以阶跃响应: = 所以，可见，实验结果与理论分析是一致的。实验3 离散系统的频率响应分析和零、极点分布实验目的：加深对离散系统的频率响应分析和零、极点分布的概念理解。实验原理：离散系统的时域方程为其变换域分析方法如下：频域 系统的频率响应为 Z域 系统的转移

9、函数为 分解因式 ，其中和称为零、极点。 在MATLAB中，可以用函数z，p，K=tf2zp（num，den）求得有理分式形式的系统转移函数的零、极点，用函数zplane（z，p）绘出零、极点分布图；也可以用函数zplane（num，den）直接绘出有理分式形式的系统转移函数的零、极点分布图。另外，在MATLAB中，可以用函数 r，p，k=residuez（num，den）完成部分分式展开计算；可以用函数sos=zp2sos（z，p，K）完成将高阶系统分解为2阶系统的串联。实验内容：求系统 的零、极点和幅度频率响应。实验要求：编程实现系统参数输入，绘出幅度频率响应曲线和零、极点分布图。实验过程

10、： 画出零极点图： num=0.0528,0.797,0.1295,0.1295,0.797,0.0528; den=1,-1.8107,2.4947,-1.8801,0.9537,-0.2336; z,p,K=tf2zp(num,den); zplane(b,a); title(零极点图); 绘制幅度频率响应曲线： num=0.0528,0.797,0.1295,0.1295,0.797,0.0528; den=1,-1.8107,2.4947,-1.8801,0.9537,-0.2336; H,W=freqz(num,den); M=abs(H); subplot(2,1,1); xlab

11、el(Mormalized Frequency(xrad/sample); ylabel(Magritude(dB); title(幅度频率响应曲线); grid; subplot(2,1,2); plot(W); xlabel(Mormalized Frequency(xrad/sample); ylabel(Phase(degree); grid; 实验4 离散信号的DTFT和DFT实验目的：加深对离散信号的DTFT和DFT的及其相互关系的理解。实验原理：序列xn 的DTFT定义：N点序列xn 的DFT定义：在MATLAB中，对形式为的DTDFT可以用函数H=Freqz（num，den，w

12、）计算；可以用函数U=fft（u，N）和u=ifft（U，N）计算N点序列的DFT正、反变换。实验内容： 分别计算16点序列 的16点和32点DFT，绘出幅度谱图形，并绘出该序列的DTFT图形。实验要求：讨论DTFT和DFT之间的相互关系。说明实验产生的现象的原因。实验步骤： = 1 * GB3 16点序列X(n)的16点及32点DFT: N=16; n=1:16; x= sin(5*pi*n/16); X1=fft(x,16); X11=abs(X1); subplot(2,1,1); stem(X11); xlabel(频率); ylabel(幅度); title(16点序列x（n）的16

13、点DFT); X2=fft(x,32); X22=abs(X2); subplot(2,1,2); stem(X22); xlabel(频率); ylabel(幅度); title(16点序列x（n）的32点DFT); = 2 * GB3 序列x（n）的DTFT: x= sin(5*pi*n/16); X1=fft(x); X11=abs(X1); stem(X11); xlabel(频率); ylabel(幅度); title(序列x（n）的DTFT); DTFT和DFT之间的区别和关系:1、DTFT是离散时间傅里叶变换，DFT是离散傅里叶变换。2、DTFT变换后的图形中的频率是一般连续的（

14、cos(wn)等这样的特殊函数除外，其变换后是冲击串)，而DFT是DTFT的等间隔抽样，是离散的点，其函数表示为X(k）,而DTFT的函数表示为（DFT是DTFT的等间隔抽样，DTFT变化后的频率响应一般是连续的，DFT变换后的频率响应是离散的）。3、DTFT是以为周期的。而DFT的序列X(k）是有限长的。4、DTFT是以复指数序列的加权和来表示的，而DFT是等间隔抽样，抽样间隔为(N为离散序列的长度)。5、DTFT和DFT都能表征原序列的信息。由于现在计算主要使用计算机，必需要是离散的值才能参与运算，因此在工程中DFT应用比较广泛，FFT是DFT的快速算法。实验5 FFT算法的应用实验目的：

15、加深对离散信号的DFT的理解及其FFT算法的运用。实验原理：N点序列的DFT和IDFT变换定义式如下：， 利用旋转因子具有周期性，可以得到快速算法（FFT）。 在MATLAB中，可以用函数X=fft（x，N）和x=ifft（X，N）计算N点序列的DFT正、反变换。实验内容：2N点实数序列N=64。用一个64点的复数FFT程序，一次算出，并绘出。（2）已知某序列在单位圆上的N=64等分样点的Z变换为。用N点IFFT程序计算，绘出和。实验要求：利用MATLAB编程完成计算，绘出相应图形。并与理论计算相比较，说明实验结果的原因。实验步骤：2N点实数序列x（n）的64点的复数FFT： k=0:N-1;

16、 n1=2*n; n2=2*n+1; x1=cos(2*pi*7*n1/N)+1/2*cos(2*pi*19*n1/N); x2=cos(2*pi*7*n2/N)+1/2*cos(2*pi*19*n2/N); XK1=fft(x1,64); XK2=fft(x2,64); XK11=XK1+(exp(-j*pi*k/N).*XK2; XK22=XK1-(exp(-j*pi*k/N).*XK2; XK11=zeros(1,N),XK11; XK22=XK22,zeros(1,N); Xk=XK11+XK22; k=0:1:2*N; mag=abs(Xk); stem(mag); xlabel(k

17、); ylabel(幅度|X(k)|); title(128点序列x（n）的64点的复数FFT);理论计算：由欧拉公式： 设，则其2N点的DFT变换为： 当时， =0 当时，即由此可得当k=14,38,90,114时有值其余为0） ， 与程序计算有相同的结论,可见，理论与实验结果相符合。 (2) X(k)的N点IFFT: N=64; n=0:N-1; Xk=1./(1-0.8.*exp(-2*pi*j*n/N); xn=ifft(Xk,N); stem(n,xn); xlabel(时间序列n); ylabel(序列x（n）); title(X(k)的ifft变换);理论计算：由，其Z变换为可得

18、，其N=64点的DFT为： = 得 其IDFT为： 则x0=1,x1=0.8,x2=0.64,x3=0.512与图像所示的xn相差很小，可以认为程序计算值与理论计算值相符。实验6 基于MATLAB的数字滤波器设计实验目的：加深对数字滤波器的常用指标和设计过程的理解。实验原理：低通滤波器的常用指标：通带边缘频率：，阻带边缘频率： ，通带起伏：，通带峰值起伏： ，阻带起伏：，最小阻带衰减： 。数字滤波器有IIR和FIR两种类型，它们的特点和设计方法不同。 在MATLAB中，可以用b，a=butter（N,Wn）等函数辅助设计IIR数字滤波器, 也可以用b=fir1(N,Wn,ftype) 等函数辅

19、助设计FIR数字滤波器。实验内容： 利用MATLAB编程设计一个数字带通滤波器，指标要求如下：通带边缘频率：，通带峰值起伏：。阻带边缘频率：，最小阻带衰减： 。分别用IIR和FIR两种数字滤波器类型进行设计。实验要求：给出IIR数字滤波器参数和FIR数字滤波器的冲激响应，绘出它们的幅度和相位频响曲线，讨论它们各自的实现形式和特点。实验步骤： = 1 * GB3 Butterworth滤波器的设计（IIR） wp=0.4*pi,0.6*pi; wr=0.3*pi,0.7*pi; Ap=1; Ar=40;N,Wn=buttord(wp/pi,wr/pi,Ap,Ar)N = 7Wn =0.3854

20、0.6146 b,a=butter(N,Wn,bandpass)b = Columns 1 through 12 0.0002 0 -0.0014 0 0.0042 0 -0.0071 0 0.0071 0 -0.0042 0 Columns 13 through 15 0.0014 0 -0.0002a = Columns 1 through 12 1.0000 0.0000 3.7738 0.0000 6.5614 0.0000 6.6518 0.0000 4.2030 0.0000 1.6437 0.0000 Columns 13 through 15 0.3666 0.0000 0.0

21、359 H,w=freqz(b,a); mag=abs(H); plot(w/pi,mag); xlabel(角频率(Omega); ylabel(幅度|Ha(jOmega)|); title(数字butterworth带通滤波器幅度响应|Ha(jOmega)|); phase=angle(H); plot(w/pi,phase); xlabel(角频率(Omega); ylabel(相位); title(数字butterworth带通滤波器相位响应曲线); = 2 * GB3 FIR滤波器的设计： wp1 = 0.4*pi; wp2 = 0.6*pi; ws1 = 0.3*pi; ws2 =

22、 0.7*pi; tr_width = min(wp1-ws1),(ws2-wp2)tr_width = 0.3142 M = ceil(6.2*pi/tr_width) + 1M = 63 n=0:1:M-1; wc1 = (ws1+wp1)/2; wc2 = (wp2+ws2)/2; wc=wc1/pi,wc2/pi; window= hanning(M); h1,w=freqz(window,1); figure(1); subplot(2,1,1) stem(window); axis(0 60 0 1.2); grid; xlabel(n); title(Hanning窗函数); s

23、ubplot(2,1,2) plot(w/pi,20*log(abs(h1)/abs(h1(1); axis(0 1 -350 0); grid; xlabel(w/pi); ylabel(幅度(dB); title(Hanning窗函数的频谱); hn = fir1(M-1,wc, hanning (M); h2,w=freqz(hn,1,512); figure(2); subplot(2,1,1) stem(n,hn); axis(0 60 -0.25 0.25); grid; xlabel(n); ylabel(h(n); title(Hanning窗函数的单位脉冲响应); subpl

24、ot(2,1,2) plot(w/pi,20*log(abs(h2)/abs(h2(1); grid; xlabel(w/pi); ylabel(幅度(dB); figure(3); phase=angle(h1); plot(phase); axis(1 pi -1 0); xlabel(w/pi); ylabel(线性相位); title(Hanning窗函数相位特性曲线); 实现形式及特点分析：在本例中，相同的技术指标下，IIR滤波器实现的阶数为N=7，而FIR滤波器的阶数N=63。因此，相同的技术指标，用IIR滤波器实现的阶数远远小于用FIR滤波器的阶数。这是由于IIR滤波器存在着输出

25、对输入的反馈。从实验中绘制的相位特性曲线可以看出，FIR滤波器可以得到严格的线性相位，而IIR滤波器做不到这一点。IIR滤波器的选择性愈好，其相位的非线性愈严重。因而，如果IIR滤波器要得到线性相位，又要满足滤波的技术要求，必须加全通网络进行相位校正，这同样会大大增加滤波器的阶数。因此，从相位特性考虑，FIR滤波器又优于IIR滤波器。附录资料：matlab画二次曲面一、螺旋线1.静态螺旋线a=0:0.1:20*pi;h=plot3(a.*cos(a),a.*sin(a),2.*a,b,linewidth,2);axis(-50,50,-50,50,0,150);grid onset(h,era

26、semode,none,markersize,22);xlabel(x轴);ylabel(y轴);zlabel(z轴);title(静态螺旋线); HYPERLINK /photo/g-OKj_1ed3vDU8433-RTuA=/876513077477666041.jpg2.动态螺旋线t=0:0.1:10*pi;i=1;h=plot3(sin(t(i),cos(t(i),t(i),*,erasemode,none);grid onaxis(-2 2 -2 2 0 35)for i=2:length(t) set(h,xdata,sin(t(i),ydata,cos(t(i),zdata,t(

27、i); drawnow pause(0.01)endtitle(动态螺旋线);(图略) 3.圆柱螺旋线t=0:0.1:10*pi;x=r.*cos(t);y=r.*sin(t);z=t;plot3(x,y,z,h,linewidth,2);grid onaxis(square)xlabel(x轴);ylabel(y轴);zlabel(z轴);title(圆柱螺旋线)HYPERLINK /photo/FWZbT1oFGgY3948_qnnNyg=/1700108859333209024.jpg 二、旋转抛物面b=0:0.2:2*pi;X,Y=meshgrid(-6:0.1:6);Z=(X.2+Y

28、.2)./4;meshc(X,Y,Z);axis(square)xlabel(x轴);ylabel(y轴);zlabel(z轴);shading flat;title(旋转抛物面)或直接用：ezsurfc(X.2+Y.2)./4) HYPERLINK /photo/1x3Uhu8ezzk8J811vbAlEA=/1700108859333209034.jpg三、椭圆柱面load clownezsurf(2*cos(u),4*sin(u),v,0,2*pi,0,2*pi)view(-105,40) %视角处理shading interp %灯光处理colormap(map) %颜色处理grid

29、on %添加网格线axis equal %使x,y轴比例一致xlabel(x轴);ylabel(y轴);zlabel(z轴);shading flat;title(椭圆柱面) %添加标题HYPERLINK /photo/QkSVYE1Dlxhfx215I259lQ=/1700108859333209038.jpg四、椭圆抛物面b=0:0.2:2*pi;X,Y=meshgrid(-6:0.1:6);Z=X.2./9+Y.2./4;meshc(X,Y,Z);axis(square)xlabel(x轴);ylabel(y轴);zlabel(z轴);shading flat;title(椭圆抛物面)或

30、直接用：ezsurfc(X.2./9+Y.2./4)HYPERLINK /photo/HC3gCvYgXN0I-ZjQdmky8Q=/1127870231680120697.jpg五、双叶双曲面ezsurf(8*tan(u)*cos(v),8.*tan(u)*sin(v),2.*sec(u),-pi./2,3*pi./2,0,2*pi)axis equalgrid onaxis squarexlabel(x轴);ylabel(y轴);zlabel(z轴);shading flat;title(双叶双曲面)HYPERLINK /photo/N7OK7I09fm9S4y70vjaLPA=/1127

31、870231680120734.jpg六、双曲柱面load clownezsurf(2*sec(u),2*tan(u),v,-pi/2,pi/2,-3*pi,3*pi)hold on %在原来的图上继续作图ezsurf(2*sec(u),2*tan(u),v,pi/2,3*pi/2,-3*pi,3*pi)colormap(map)shading interpview(-15,30)axis equalgrid onaxis equalxlabel(x轴);ylabel(y轴);zlabel(z轴);shading flat;title(双曲柱面)HYPERLINK /photo/cogod3o

32、KNYaofmvLfknS2A=/4254212798005282624.jpg七、双曲抛物面（马鞍面）X,Y=meshgrid(-7:0.1:7);Z=X.2./8-Y.2./6;meshc(X,Y,Z);view(85,20)axis(square)xlabel(x轴);ylabel(y轴);zlabel(z轴);shading flat;title(双曲抛物面)或直接用：ezsurfc(X.2./8-Y.2./6) HYPERLINK /photo/dlBlUM1Za9f3-UbeXV2WuA=/297519050383943404.jpg八、抛物柱面X,Y=meshgrid(-7:0.

33、1:7);Z=Y.2./8;h=mesh(Z);rotate(h,1 0 1,180) %旋转处理%axis(-8,8,-8,8,-2,6);axis(square)xlabel(x轴);ylabel(y轴);zlabel(z轴);shading flat;title(抛物柱面)或直接用：ezsurfc(Y.2./8) HYPERLINK /photo/MJVd3Cv8LW_rtuHyLkHfbA=/291045125919324104.jpg九、环面ezmesh(5+2*cos(u)*cos(v),(5+2*cos(u)*sin(v),2*sin(u),0,2*pi,0,2*pi)axis

34、equalgrid onxlabel(x轴);ylabel(y轴);zlabel(z轴);shading flat;title(环面)HYPERLINK /photo/5EdAokqXE2fmz8FYN0-EkA=/291045125919324112.jpg十、椭球ezsurfc(5*cos(u)*sin(v),(3*sin(u)*sin(v),4*cos(v),0,2*pi,0,2*pi)axis equalgrid onxlabel(x轴);ylabel(y轴);zlabel(z轴);shading flat;title(椭球)HYPERLINK /photo/-jAH-p581gJh0btZzQHuEw=/291045125919324147.jpg十一、单叶双曲面ezsurf(4*sec(u)*