matlab.doc

华 北 电 力 大 学 实 验 报 告实验环境MATLAB7.0实验名称实验一 Matlab基本操作实 验 目 的1、认识Matlab软件环境，熟悉Matlab编程的特点。 2、学习信号的产生与画图，矩阵生成与运算，流程控制等操作。实 验 内 容1、编写程序产生如下信号，并按要求绘制相应的曲线。利用saveas命令把图形分别保存成.fig格式和.jpg格式的图像文件。（1）信号形式，利用plot命令绘制该信号在区间内，且采样间隔为的曲线。填加网格线，标题为“正弦信号曲线”，横坐标显示“x”，纵坐标显示“y”。（2）将上述信号左平移或右平移得到两个新信号y2和y3，利用hold命令在同一图中显示三条曲线并通过线型和颜色加以区分。（3）利用subplot命令将上述三个信号利用子图来显示；利用figure命令将上述三个信号分别在不同的窗体中同时显示。（4）z = 0:0.1:40；x = cos(z)；y = sin(z)，利用plot3命令绘制三维曲线并填加相应的坐标轴标注。（通信选做）试一试，如何添加图注和为曲线在指定位置添加注释。2、向量、矩阵生成及相关运算。（1）x=1:3; y=4:6; x*y=?; x.*y=?; x*y=?; x*y=?; x.*y=?; x*y.=?（2）创建矩阵，并进行如下操作，取出矩阵A的前两行元素，生成一个新矩阵B。生成一个单位矩阵，与矩阵A进行加减运算。练习利用save命令保存数据为mat格式，并利用load命令从外部文件读入数据。（3）利用wavread 与imread命令读取一个wav格式的语音文件和jpg格式的图像文件，并利用sound和image/imshow命令将其播放或显示出来。利用auwrite把语音文件写成.au格式的文件，利用imwrite把图像文件写成bmp格式的文件。实 验 内 容（4）利用矩阵求逆和乘法运算求解方程组想一想，进行矩阵的运算时，其维数需满足什么要求。3、流程控制语句相关操作。（1）使用for语句及while语句求解1至100整数中奇数的和。（2）求连续自然数的和，当和大于等于1000时，求最后一个自然数以及自然数的和。（3）Matlab函数是可以递归调用的，即在函数的内部调用函数自身，试用递归方式编写n!函数。（选做）用Matlab编程建立九九乘法表。实 验 结 果 及 分 析1、程序：（1）（3）x=0:pi/1000:2*pi;y=sin(2*x+pi/4);plot(x,y)grid ontitle(正弦信号曲线)xlabel(x)ylabel(y)hold onxl=x-pi/2;xr=x+pi/2;y2=sin(2*xl+pi/4);y3=sin(2*xr+pi/4);plot(xl,y2,r :)hold onplot(xr,y3,y -)subplot(1,3,1)plot(x,y)subplot(1,3,2)plot(xl,y2)subplot(1,3,3)plot(xr,y3)figure(2)plot(x,y)figure(3)plot(xl,y2)figure(4)plot(xr,y3)实 验 结 果 及 分 析程序：（4）z=0:0.1:40; x=cos(z);y=sin(z);plot3(x,y,z);grid on;xlabel(x)ylabel(y)zlabel(z)saveas(h,filename,format)结果：（1）（3）结果：（4）第 24 页实 验 结 果 及 分 析答：有关图形标注函数的调用格式为： title(图形名称) xlabel(x轴说明) ylabel(y轴说明) text(x,y,图形说明) legend(图例1,图例2,)2、（1）程序：x=1:3;y=4:6;c=dot(x,y);x.*yx*yx*yx.*yx*y.结果：c =32ans =4 10 18ans = 4 5 68 10 1212 15 18ans =32ans = 4 5 68 10 1212 15 18ans =32（2）程序： A=1 2 3;4 5 6;7 8 9B=A(1:2,:)C=eye;x=A+Cy=A-Csaveload Untitled.txtUntitled结果：A = 1 2 3 4 5 6 7 8 9B = 1 2 3 4 5 6实 验 结 果 及 分 析x = 2 3 4 5 6 7 8 9 10y = 0 1 2 3 4 5 6 7 8Saving to: matlab.matUntitled = 1 2 3 4 5 6 7 8 9（3）程序： x=wavread(1)plot(x)saveload()x=imread(G:05.jpg);image(x) imwrite(x,.bmp)结果：（4）程序：A = -1 1 2; 3 -1 1; -1 3 4;b = 2;6;4;x = inv(A)*b 结果：x = 1.0000 -1.00002.0000实 验 结 果 及 分 析答:进行矩阵的运算时，需保证两个矩阵的维数一致或其中一个矩阵为标量。3、（1）程序：i=1;sum=0;while(i100) sum=sum+i; i=i+2;enddisp(sum)结果：2500程序：sum=0;fori=1:2:99 sum=sum+i;endstr=计算结果为,num2str(sum);disp(str)结果：计算结果为：2500（2）程序：sum=0;n=0;while sum1000 n=n+1; sum=sum+n;enddisp(n)disp(sum)结果：45 1035（3）程序：(设n=100)clear alln=1;sum=1;while n100;n=n+1;sum=sum*n;end结果:n =100sum =9.3326e+157实 验 结 果 及 分 析程序：disp(九九乘法表 ), for i=1:9, for j=1:i,pi(j)=j*i;end， disp(pi);end 结果：九九乘法表 1 2 4 3 6 9 4 8 12 16 5 10 15 20 25 6 12 18 24 30 36 7 14 21 28 35 42 49 8 16 24 32 40 48 56 64 9 18 27 36 45 54 63 72 81实验分析：1、 MATLAB具备高性能的数值计算和可视化的科学工程计算功能，支持解释性语言输入，编程实现简单，具有丰富的数学函数功能支持。MATLAB允许与、Fortran语言接口，其部件Simulink可以采用图形输入的方式来搭构所研究的系统。包含丰富的工具包，在系统仿真、数字信号处理、图形图像分析、数理统计、通信及自动控制领域得到广泛的应用。2、 矩阵的建立方法有直接输入法，外部文件读入法，利用冒号表达式建立一个向量和建立大矩阵。直接输入矩阵时要以“”为其标志符号，矩阵的所有元素都在括号内。外部文件读入可利用load函数：load+文件名+参数,和直接导入数据。冒号表达式可以建立一个行向量，一般格式为e1:e2:e3。大矩阵可由方括号中的小矩阵和向量建立起来。3、 矩阵的基本数学运算“+、-、*”运算法则与普通数字的运算几乎相同。基本函数运算包括矩阵逆运算、行列式运算、秩运算、特征值运算。4、 MATLAB提供了多种控制程序流的的结构：if-else-end分支结构、switch-case选择结构、for循环结构、while循环结构等。实验名称实验二 Matlab在信号与系统中的应用实 验 目 的1、掌握信号与系统中Matlab编程的特点。 2、学习微分方程、卷积、系统函数的Matlab编程实现。实 验 内 容1、按照卷积原理步骤求如下两幅图中和信号卷积后的波形。2、给定x(n)=nR16(n)，h(n)=R8(n)，求x和h卷积y，并用函数stem画出x、h和y的图形。3、用Matlab在时域中编程求解 y(t)+y(t)=f(t), f(t)= exp(-3t)(t)的冲激响应、阶跃响应及零状态响应。在simulink仿真环境下，设计系统框图，分析系统的冲激响应、阶跃响应。实 验 结 果 及 分 析1、程序：t1_s=-2;t2_s=2;t1_e=0;t2_e=4;t_s=t1_s+t2_s;t_e=t1_e+t2_e;t0=min(t1_s,t2_s,t_s)-1;t1=max(t1_e,t2_e,t_e)+1;dt=0.005;t=t0:dt:t1;f1=2*rectpuls(t,2);f2=t.*rectpuls(t,4);L=length(t);tp=2*t(1):dt:2*t(L);fl=eval(f1);f2=eval(f2);y=dt*conv(fl,f2);y_max=max(y);y_min=min(y);dy=(y_max-y_min)/10;f1_max=max(f1);f1_min=min(f1);df1=(f1_max-f1_min)/10;f2_max=max(f2);f2_min=min(f2);df2=(f2_max-f2_min)/10;subplot(3,1,1);实 验 结 果 及 分 析plot(t,fl,linewidth,2);title(信号f1(t)的波形,color,b,FontSize,8)axis(t(1),t(L),f1_min-df1,f1_max+df1);grid;set(gca,Fontsize,8)subplot(3,1,2);plot(t,f2,linewidth,2);title(信号f2(t)的波形,color,b,Fontsize,8)axis(t(1),t(L),f2_min-df1,f2_max+df1);grid;set(gca,FontSize,8)subplot(3,1,3);% plot(tp,y,1inewidth,2); plot(tp,y);title(卷积f(t)=f1*f2的波形,color,b,FontSize,8)if y_min=0 y0=(y_max-y_min)/2;elsey0=0;endaxis(t(1),t(L),y_min-dy,y_max+dy);grid; set(gca,Ytick,y_min,y0,y_max,FontSize,8);结果：实 验 结 果 及 分 析2、程序：程序：n=0:15;x=n;subplot(3,1,1)stem(n,x)grid onn0=0:7;h=ones(1,8);subplot(3,1,2)stem(n0,h)grid ony=conv(x,h)n1=1:23;subplot(3,1,3)stem(n1,y)结果：y = Columns 1 through 18 0 1 3 6 10 15 21 28 36 44 52 60 68 76 84 92 84 75 Columns 19 through 23 65 54 42 29 153、程序：a=1 1;b=1;p=0.5;t=0:p:5; x=exp(-3*t);subplot(311)实 验 结 果 及 分 析impulse(b,a,0:p:5);title(冲激响应);subplot(312)step(b,a,0:p:5);title(阶跃响应);subplot(313)lsim(b,a,x,t);title(零状态响应);结果：仿真框图：实 验 结 果 及 分 析仿真结果：实验分析：1、 在MATLAB中连续信号可用向量或符号运算功能来实现。MATLAB为用户提供了专门用于求连续系统冲击响应及阶跃响应，并绘制其时域波形的函数impulse和step。2、 函数impulse()将绘出由a和b表示的连续系统在指定时间范围内的冲激响应h（t）的时域波形图，并求出指定时间范围内冲激响应的数值解。函数step（）将绘出由向量a和b表示的连续系统的阶跃响应g（t）在指定时间范围内的波形图，并能求出其数值解。3、 可以通过MATLAB命令来打开Simulink模型，向模型传递数据并进行仿真。Simulink仿真的数据结果也可以送回MATLAB工作空间中作进一步数值分析。具体步骤：启动Simulink。在命令窗中输入Simulink或点击工具栏中的按钮。打开编辑窗。根据具体问题建立数学模型。确定需要的仿真模块。对模块进行编辑，以方便模型的建立。连线。实验名称实验三 Matlab在数字信号处理中的应用实 验 目 的1、学习DFT和FFT的原理及Matlab编程实现方法。 2、学习FFT应用于卷积运算和谱分析时的Matlab编程实现方法。实 验 内 容1、用三种不同的DFT程序计算x(n)=R8(n)的傅里叶变换X(ejw),并比较三种程序计算机运行时间。（1）用for loop 语句的M函数文件dft1.m，用循环变量逐点计算X(k)；（2）编写用MATLAB矩阵运算的M函数文件dft2.m, 完成下列矩阵运算；（3）调用FFT库函数，直接计算X(k)；（4）分别利用上述三种不同方式编写的DFT程序计算序列x(n)的傅立叶变换X(ejw)，并画出相应的幅频和相频特性，再比较各个程序的计算机运行时间。2、研究高密度频谱与高分辨率频谱设有连续信号Xa(t)=cos(26.5103t)+cos(27103t)+cos(29103t)以采样频率fs=32kHz对信号Xa(t)采样,分析下列三种情况的幅频特性。（1） 采集数据长度N=16点，做N=16点的DFT，并画出幅频特性。（2） 采集数据长度N=16点，补零到256点，做N=256点的DFT，并画出幅频特性。（3） 采集数据长度N=256点，做N=256点的DFT，并画出幅频特性。观察三幅不同频率特性图，分析和比较它们的特点以及形成的原因。实 验 结 果 及 分 析1、程序：dft1.m:functionAm,pha=dft1(x)N=length(x);w=exp(-j*2*pi/N);for k=1:N sum=0; for n=1:N sum=sum+x(n)*w(k-1)*(n-1); end Am(k)=abs(sum); pha(k)=angle(sum);enddft2.mfunctionAm,pha=dft2(x)N=length(x);n=0:N-1;k=0:N-1;实 验 结 果 及 分 析w=exp(-j*2*pi/N);nk=n*k;wnk=w.(nk);Xk=x*wnk;Am=abs(Xk); pha=angle(Xk); dft3.mfunctionAm,pha=dft3(x)Xk=fft(x);Am=abs(Xk);pha=angle(Xk);dft1-main.mx=ones(1,8),zeros(1,248);t=cputime;Am1,pha1=dft1(x);t1=cputime-tn=0:(length(x)-1);w=(2*pi/length(x)*n;figure(1)subplot(2,1,1), plot(w,Am1,b); grid;title(Magnitude part);xlabel(frequency in radians);ylabel(|X(exp(jw)|);subplot(2,1,2), plot(w,pha1,r); grid;title(Phase Part);xlabel(frequency in radians);ylabel(argXexp(jw)/radians);dft2-main.mx=ones(1,8),zeros(1,248); t=cputime;Am2,pha2=dft2(x);t2=cputime-tn=0:(length(x)-1);w=(2*pi/length(x)*n;figure(2)subplot(2,1,1), plot(w,Am2,b); grid;title(Magnitude part);xlabel(frequency in radians);ylabel(|X(exp(jw)|);实 验 结 果 及 分 析subplot(2,1,2), plot(w,pha2,r); grid;title(Phase Part);xlabel(frequency in radians);ylabel(argXexp(jw)/radians);dft3-main.mx=ones(1,8),zeros(1,248);t=cputime;Am3,pha3=dft3(x); t3=cputime-t;n=0:(length(x)-1);w=(2*pi/length(x)*n;figure(3)subplot(2,1,1), plot(w,Am3,b); grid;title(Magnitude part);xlabel(frequency in radians);ylabel(|X(exp(jw)|);subplot(2,1,2), plot(w,pha3,r); grid;title(Phase Part);xlabel(frequency in radians);ylabel(argXexp(jw)/radians);结果：t1= 0.2813t2= 0.2031 实 验 结 果 及 分 析（4）从以上运行结果可以看出，调用FFT库函数直接计算X(k)速度最快，矩阵运算次之，用循环变量逐点计算运行速度最慢。因此FFT算法大大提高了DFT的实用性。2、程序：fs=32000; N=16; n=0:N-1; xa=cos(2*pi*6.5*103*n/fs)+cos(2*pi*7*103*n/fs)+cos(2*pi*9*103*n/fs);Am1,pha1=dft3(xa);n=0:(length(xa)-1);w=(2*pi/length(xa)*n;figure(1)plot(w,Am1,b); grid;title(Magnitude part);xlabel(frequency in radians);ylabel(|X(exp(jw)|);x=xa(1:16),zeros(1,240); Am2,pha2=dft2(x); n=0:(length(x)-1);w=(2*pi/length(x)*n;figure(2)plot(w,Am2,b); grid;title(Magnitude part);xlabel(frequency in radians);ylabel(|X(exp(jw)|);N0=256; n0=0:N0-1; xa0=cos(2*pi*6.5*103*n0/fs)+cos(2*pi*7*103*n0/fs)+cos(2*pi*9*103*n0/fs);Am3,pha3=dft3(xa0);n=0:(length(xa0)-1);w=(2*pi/length(xa0)*n;figure(3)plot(w,Am3,b); grid;实 验 结 果 及 分 析title(Magnitude part);xlabel(frequency in radians); ylabel(|X(exp(jw)|);结果：答：观察运行结果可知，N=16点时的DFT，由于点数太少，很难反映出频谱的细节特征，只能分辨出两个频率分量。将16点DFT补零到256点，做256点DFT，由结果可以看出，频率间隔缩小，得到了比较平滑的连续曲线，此为高密度频谱。它是由16点经过补零所得到的，并没有增加任何信息，频率分辨率并未提高，仍然只能看出两个频率分量。第三次做的是256点的DFT，采集点数大大增加，分辨率提高，可以清楚看到三根谱线的位置，是为高分辨率频谱。实验分析：1、 无限长的序列往往可以用有限长序列来逼近，对于有限长的序列我们可以使用离散傅立叶变换（DFT）。DFT是对序列傅立叶变换的等距采样。2、 快速傅立叶变换FFT是为了减少DFT运算次数的一种快速算法。它是对DFT变换式进行一次次的分解，使其成为若干小点数DFT的组合，从而减小运算量。常用的FFT是以2为基数，它的运算效率高，程序比较简单，使用也十分地方便。FFT的算法基本上可以分为两大类：按时间抽取（DIT）和按频率抽取（DIF）。3、 MATLAB的Symbolic Math Toolbox提供了能直接求解傅立叶变换及逆变换的函数fourier()和ifourier()。在绘制离散信号波形时，要使用专门绘制离散数据的stem命令。还提供了专门求离散系统单位响应，并绘制其时域波形的函数impz()。实验名称实验四 Matlab在通信系统中的应用实 验 目 的1、练习通信原理中各种调制方法及简单通信系统的Matlab编程实现；2、学习Matlab与Simulink的混合编程。实 验 内 容1用Matlab产生一个频率为2Hz、功率为1的正弦信源m(t)，设载波频率为10Hz，A分别等于0.5、1、2，试画出：（1）AM调制信号；（2）调制信号的功率谱密度；（3）相干解调后的信号波形。（1）通过得到的曲线，分析调制深度与失真的关系。（2）尝试用DSB-SC调制方法重做上题。（3）绘出加入高斯白噪声后的解调信号。2消息信号，采用频率调制中心频率为100Hz的载波。频偏常数为。（1）画出该调频信号的波形。（2）求出该调频信号的振幅谱。（3）用鉴频器解调该调频信号，并与输入信号比较。确定调制系数，加载白噪声后，重新绘出消息信号和解调信号的频谱及解调器输出的功率谱密度。3在Matlab与Simulink混合仿真环境下，建立简单通信系统的仿真模型，其中信源是间隔为0.0001s的随机二进制数，由Matlab编程封装成Simulink模块，采用BPSK调制与解调，经过AWGN信道。利用Matlab编程调用此Simulink文件，并绘制信噪比为0、2、4、6、8、10时的误码率曲线。试一试，观察信号眼图。实 验 结 果 及 分 析1、程序：dt = 0.001; fm = 1; fc = 10; T = 5; t = 0:dt:T;mt = sqrt(2) * cos(2 * pi * fm * t); A = 2;s_am = (A + mt).* cos(2 * pi * fc * t);实 验 结 果 及 分 析B = 2 * fm; figure(1)subplot(311)plot(t,s_am);hold on; plot(t,A + mt,r-);title(AM调制信号及其包络);xlabel(t); rt = s_am.*cos(2 * pi * fc * t); rt = rt - mean(rt);f,rf = T2F(t,rt);t,rt = lpf(f,rf,2 * fm); subplot(312)plot(t,rt);hold on;plot(t,mt/2,r-);title(相干解调后的信号波形与输入信号的比较);xlabel(t)subplot(313)f,sf=T2F(t,s_am); psf = (abs(sf).2)/T; plot(f,psf);axis( -2*fc 2*fc 0 max(psf);title(AM信号功率谱);xlabel(f); 结果：A=0.5 A=1实 验 结 果 及 分 析A=22、程序：Kf = 5;fc = 10;T = 5;dt = 0.001; t = 0:dt:T;fm = 1;mt = cos(2 * pi * fm * t) + sin(2 * pi * 0.5 * fm * t); A =sqrt(2);mti = 1/2/pi/fm * sin(2 * pi * fm * t) - 1/2/pi/0.5/fm * cos(2 * pi * 0.5 * fm* t);% mti = 1/2/pi/fm * sin(2 * pi * fm * t); st = A * cos(2 * pi * fc * t + 2 * pi * Kf * mti);figure(1)subplot(311)plot(t,st);hold on; plot(t, mt,r-);xlabel(t); ylabel(调频信号);subplot(312)f,sf = T2F(t,st);plot(f, abs(sf);axis( -25 25 0 3)xlabel(f); ylabel(调频信号幅度谱);% FM demodulationfor k = 1:length(st) - 1 rt(k) = (st(k+1) - st(k)/dt;endrt(length(st) = 0;subplot(313)实 验 结 果 及 分 析plot(t,rt);hold on;plot(t, A * 2 * pi * Kf * mt + A * 2 * pi * fc, r-);xlabel(t);ylabel(调频信号微分后包络)结果：3、仿真框图： 结果：实 验 结 果 及 分 析答：（1）调制深度不会超过1，否则包络就会出现失真，信号会出现额外的频谱，造成干扰信号。（2）FS = 128;t = 0 : 1/FS : 10-1/FS ;sgma=sqrt(0.001 0.01 0.05 0.1 0.3);noise=zeros(length(sgma),length(t);for k=1:length(sgma) noise(k,:)=normrnd(0,sgma(k),1,length(t);endx1= sawtooth(pi*t ,0.5);x1=(x1+1)/2;%plot(t,x1),axis(0 10 0 2);x2= square(pi*(t-0.1),90);x2=(x2+1)/2;%plot(t,x2),axis(0 10 0 2);x=x1.*x2;figure(1);subplot(3, 1 ,1),plot(t,x),axis(0 5 -1 1);s = sin(100*pi*t)subplot(3, 1 ,2),plot(t,s), axis(0 5 -1 1)X = x.* ssubplot(3, 1 ,3),plot(t,X), axis(0 10 -1 1)N=512;Y=fft(X,N);YY=Y(1:256);Y=Y(257:512),YY;n=-N/2+1:N/2;%f=FS*n/512;figure(2);k=(0:N-1)/512*128; subplot(2,1,1),plot(k-64,abs(Y)%axis(0 64 0 140)Y1=fft