Matlab实验仿真实验数据
发布时间：2013-10-05 17:45:27
附录
实验二
程序：
（1）
>> clear
>> tic;
>> t=-5:0.5:5;
>> for n=1:size(t,2)
if(t(n)<0)
y(n)= 3*t(n)^2+5;
else
y(n)= -3*t(n)^2+5;
end
end
>> figure(1);
>> plot(t,y);
>> xlabel('x');
>> ylabel('y');
>> grid on;
>> toc;
（2）
>> clear
>> tic;
>> t=[-5:0.5:5];
>> b=t>=0;
>> y(b)=-3*t(b).^2 + 5;
>> y(~b)=3*t(~b).^2 + 5;
>> figure(2);
>> plot(t,y);
>> xlabel('x');
>> ylabel('y');
>> grid on;
>> toc;
结果：
（1）Elapsed time is 0.156000 seconds.
（2）Elapsed time is 0.094000 seconds.
实验三
程序：
（1）
>> clear;
>> n=input('ENTER A NUMBER:');
>> sum=0;
>> m=1;
>> while m
sum=sum+m;
m=m+2;
end
>> fprintf('The result of all odd numbers within a given number is:%d\n',sum);
（2）创建Fib.m文件
% 函数功能: 计算斐波那契数列的第 n 个斐波那契数
% 文 件 名： Fib.m
% 含有 n 个数的斐波那契数列的定义如下：
% f(1) = 1
% f(2) = 2
% f(n) = f(n-1) + f(n-2)
function y=Fib(n);
a(1)=1;a(2)=1;i=2;
while i<=n
a(i+1)=a(i-1)+a(i);
i=i+1;
end;
y=a(i);
结果：
（1）ENTER A NUMBER:6
The result of all odd numbers within a given number is:9
（2）
>> Fib(7)
ans =
21
实验四
（1）
程序：
创建myfun.m文件
% 函数功能: 计算x的双曲正弦、双曲余弦和双曲正切,并画出对应的图象。
% 文 件 名： myfun.m
function y=myfun();
x=-pi:pi/200:pi;
sinh=(exp(x)-exp(-x))/2;
cosh=(exp(x)+exp(-x))/2;
tanh=(exp(x)-exp(-x))/(exp(x)+exp(-x));
subplot(2,2,1);
plot(x,sinh);
xlabel('x');
ylabel('y');
title('双曲正弦');
subplot(2,2,2);
plot(x,cosh);
xlabel('x');
ylabel('y');
title('双曲余弦');
subplot(2,2,3);
plot(x,tanh);
xlabel('x');
ylabel('y');
title('双曲正切');
subplot(2,2,4);
plot(x,sinh,'r',x,cosh,'g',x,tanh,'b');
xlabel('x');
ylabel('y');
title('双曲正弦/余弦/正切');
legend('sinh','cosh','tanh');
end
结果：
（2）
程序：
创建maxminfun.m文件
% 函数功能: 计算函数的最值。
% 文 件 名： maxminfun.m
function [xmin,min_value,xmax,max_value]=maxminfun(first_value,num_steps,last_value,func)
x=first_value:num_steps:last_value;
y=feval(func,x); %feval--执行待定的函数
xmin=x(1);
xmax=x(1);
min_value=y(1);
max_value=y(1);
for n=1:size(y,2)
if(y(n)>max_value)
xmax=x(n);
max_value=y(n);
end
if(y(n)
xmin=x(n);
min_value=y(n);
end
end
创建func.m文件
% 函数功能: 函数func()。
% 文 件 名： func.m
function y=func(i)
y=i.^2+2;
end
结果：
>> [xmin,min_value,xmax,max_value] =maxminfun(-2,2,4,'func')
xmin =
0
min_value =
2
xmax =
4
max_value =
18
实验五
程序： 结果：
>> clear;
>> t=-pi:0.1:pi
>> v=2*sin(3*t);
>> subplot(2,2,1);
>> plot(t,v);
>> title('plot(t,v)');
>> subplot(2,2,2);
>> plot(v);
>> title('plot(v)');
>> subplot(2,2,3);
>> plot(t,v, 'g','LineWidth',3);
>> subplot(2,2,4);
>> plot(v,'*');
实验六
程序：
>> clear
>> fid=fopen('file_.dat','w');
>> x=linspace(-10,10,30);
>> fprintf(fid,'%f\t',x);
>> status=fclose(fid);
>> fid=fopen('file_.dat','r');
>> [array,count]=fscanf(fid,'%f',8)
>> status=fclose(fid);
结果：
（1）
array =
-10.0000
-9.3103
-8.6207
-7.9310
-7.2414
-6.5517
-5.8621
-5.1724
count = 8
（2）
实验七
程序：
（1）
>> clear
>> x = eye(5,6);
>> y = sparse(x);
>> disp(x);
>> disp(y);
（2）
>> clear
>> a{1,2} = {1 3 -7;2 0 6;0 5 1};
>> a{1,1} = {'This is a text string.'};
>> a{2,2} = {3+4*i -5;-10*i 3-4*i};
>> a{2,1} = [];
>> celldisp(a);
（3）
>> clear
>> student = struct('name', {'LiMing','WangXiao'},'addr1', {'QDLG'}, 'city', {'Jinan','Qingdao'});
>> disp('student1');
>> disp(student(1));
>> disp('student2');
>> disp(student(2));
结果：
（1）
1 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0
0 0 0 1 0 0
0 0 0 0 1 0
(1,1) 1
(2,2) 1
(3,3) 1
(4,4) 1
(5,5) 1
（2）a{1,1}{1} = This is a text string.
a{2,1} = []
a{1,2}{1,1} = 1 a{1,2}{2,1} = 2 a{1,2}{3,1} = 0
a{1,2}{1,2} = 3 a{1,2}{2,2} = 0 a{1,2}{3,2} = 5
a{1,2}{1,3} = -7 a{1,2}{2,3} = 6 a{1,2}{3,3} = 1
a{2,2}{1,1} = 3.0000 + 4.0000i a{2,2}{2,1} = 0 - 10.0000i
a{2,2}{1,2} = -5 a{2,2}{2,2} = 3.0000 - 4.0000i
（3）
student1
name: 'LiMing'
addr1: 'QDLG'
city: 'Jinan'
student2
name: 'WangXiao'
addr1: 'QDLG'
city: 'Qingdao'
实验一 阶跃响应与冲激响应
实验硬件布局图
1. 微分电路
程序：
>> clear;
>> C=47*10^(-10);R=1000;
>> aa1=1/R/C;
>> a=[1 aa1]; %系统左端模型系数
>> b=[1 0]; %系统右端模型系数
>> t=linspace(0,0.001);
>> e=(t+1)./(t+1); %激励信号
>> y=exp(-t/(R*C)); %输出信号
>> s1=step(b,a,t); %根据模型利用matlab函数直接求阶跃响应
>> plot(t,y,t,e,'-*',t,s1,'--');
>> xlabel('Time(s)');
>> ylabel('output voltage (V)');grid on;
结果：
2.二阶电路阶跃响应
程序：
>> clear;
>> L=0.01; %电感10mh
>> C=1*10^(-7); %电容10×104PF
>> R=2*sqrt(L/C); %电阻 取临界阻尼值
>> aa1=R/L; aa2=1/(C*L);
>> a=[1 aa1 aa2]; %微分方程左侧系数
>> b=aa2; %微分方程右侧系数
>> t=linspace(0,0.001); %在0到0.001秒之间
>> s1=step(b,a,t); %求系统的单位阶跃响应
>> R=2*sqrt(L/C)/10; %电阻取临界状态时的1/10
>> aa1=R/L;
>> a1=[1 aa1 aa2];
>> s2=step(b,a1,t);
>> R=2*sqrt(L/C)*2; %电阻取临界状态时的2倍
>> aa1=R/L;
>> a1=[1 aa1 aa2];
>> s3=step(b,a1,t);
>> plot(t,s1,t,s2,'--',t,s3,'*-');
>> legend('临界','欠阻尼','过阻尼') ;
>> grid on;
结果：
3.二阶电路冲激响应
程序：>> clear;
>> L=0.01; %电感10mh
>> C=1*10^(-7); %电容10×104PF
>> R=2*sqrt(L/C); %电阻 取临界阻尼值 (R=632.45)
>> aa1=R/L; aa2=1/(C*L);
>> a=[1 aa1 aa2]; %微分方程左侧系数
>> b=aa2; %微分方程右侧系数
>> t=linspace(0,0.001); %在0到0.001秒之间
>> s1=impulse(b,a,t); %求系统的单位阶跃响应
>> R=2*sqrt(L/C)/10; %电阻取临界状态时的1/10
>> aa1=R/L;
>> a1=[1 aa1 aa2];
>> s2=impulse(b,a1,t);
>> R=2*sqrt(L/C)*2; %电阻取临界状态时的2倍
>> aa1=R/L;
>> a1=[1 aa1 aa2];
>> s3=impulse(b,a1,t);
>> plot(t,s1,t,s2,'--',t,s3,'*-');
>> legend('临界','欠阻尼','过阻尼'); grid on;
结果：
4.SIMULINK仿真
二阶电路仿真框图模型
5.GUI仿真
过阻尼
欠阻尼
临界状态
实验二 IIR数字滤波器实现
附程：
1．巴特沃思型模拟低通滤波器设计函数：
function[b,a]=afd_butt(Wp,Ws,Rp,As)
if Wp<=0
error('通带边缘必须大于0')
end
if Ws<=0
error('阻带边缘必须大于通带边缘')
end
if (Rp<=0)|(As<0)
error('通带波动或阻带衰减必须大于0')
end
N=ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(Wp/Ws))); %求滤波器的阶次N
fprintf('\n***Butterworth Filter Order=%2.0f\n',N)
OmegaC=Wp/((10^(Rp/10)-1)^(1/(2*N))) %计算截止频率
[b,a]=u_buttap(N,OmegaC);
2．未归一化巴特沃思型模拟低通滤波器原型设计函数：
function[b,a]=u_buttap(N,Omegac)
[z,p,k]=buttap(N); %创建Buttord低通滤波器原型
p=p*Omegac;
k=k*Omegac^N;
B=real(poly(z)); %求原型滤波器系数B
bo=k;
b=k*B;
a=real(poly(p));
3．计算滤波器的幅值响应（绝对、相对）、相位响应、群延迟函数：
function[db,mag,pha,grd,w]=freqz_m(b,a);
[H,w]=freqz(b,a,1000,'whole');
H=(H(1:1:501))';w=(w(1:1:501))';
mag=abs(H);
db=20*log10((mag+eps)/max(mag)); %化为分贝值
pha=angle(H);
grd=grpdelay(b,a,w);
4．将结果绘制成图形函数
function y=ht(b,a,Wp,Ws,Attn,Ripple)
figure(1);
[db,mag,pha,grd,w]=freqz_m(b,a); %求模拟域频率响应
subplot(2,2,1);plot(w/pi,mag);title('Magnitude Response')
xlabel('frequency in pi units');ylabel('magnitude');axis([0,1,0,1.1])
set(gca,'XTickMode','manual','XTick',[0,Wp/pi,Ws/pi,1]);
set(gca,'YTickMode','manual','YTick',[0,Attn,Ripple,1]);grid
subplot(2,2,3);plot(w/pi,db);title('Magnitude in dB')
xlabel('frequency in pi units');ylabel('decibels');axis([0,1,-40,5])
set(gca,'XTickMode','manual','XTick',[0,Wp/pi,Ws/pi,1]);
set(gca,'YTickMode','manual','YTick',[-40,-15,-5,0]);grid
set(gca,'YTickLabelMode','manual','YTickLabel',['40';'15';' 5';' 0'])
subplot(2,2,2);plot(w/pi,pha/pi);title('Phase Response')
xlabel('frequency in pi units');ylabel('pi units');axis([0,1,-1,1])
set(gca,'XTickMode','manual','XTick',[0,Wp/pi,Ws/pi,1]);
set(gca,'YTickMode','manual','YTick',[-1,0,1]);grid
subplot(2,2,4);plot(w/pi,grd);title('Groupe Delay')
xlabel('frequency in pi units');ylabel('Samples');axis([0,1,0,10])
set(gca,'XTickMode','manual','XTick',[0,Wp/pi,Ws/pi,1]);
set(gca,'YTickMode','manual','YTick',[0:2:10]);grid
5.Impulse Invariance Transformation from Analog to Digital Filter
function [b,a] = imp_invr(c,d,T)
% [b,a] = imp_invr(c,d,T)
% b = Numerator polynomial in z^(-1) of the digital filter
% a = Denominator polynomial in z^(-1) of the digital filter
% c = Numerator polynomial in s of the analog filter
% d = Denominator polynomial in s of the analog filter
% T = Sampling (transformation) parameter
%
[R,p,k] = residue(c,d);
p = exp(p*T);
[b,a] = residuez(R,p,k);
b = real(b'); a = real(a');
6.用冲击响应不变法和双线性变换法将模拟滤波器变换为数字滤波器
，采样周期为
。
程序：>> clear
>> b=[3 2];
>> a=[2 3 1];
>> Ts=0.1;
>> [R,Ps,K]=residue(b,a);
>> Pz=exp(Ps*Ts);
>> disp('Use direct principle:')
>> [bz,az]=residue(R,Pz,K)
>> disp('Use function "IMPINVAR":')
>> [bz1,az1]=impinvar(b,a,1/Ts)
>> disp('Use function "BILINEAR":')
>> [bz2,az2]=bilinear(b,a,1/Ts)
结果：
冲击响应不变法：
>> Use direct principle:
bz =
1.5000 -1.4036
az =
1.0000 -1.8561 0.8607
>> Use function "IMPINVAR":
bz1 =
0.3000 -0.2807
az1 =
2.0000 -3.7121 1.7214
双线性变换法：
>> Use function "BILINEAR":
bz2 =
0.0720 0.0046 -0.0674
az2 =
1.0000 -1.8560 0.8606
7.用冲击响应不变法和双线性变换法设计巴特沃思原型，使满足：
ωp=0.2π Rp=1dB
ωs=0.3π As=15dB
T=1
（1）冲击响应不变法：
程序：
Wp=input('请输入通带截止频率Wp=');Rp=input('请输入通带波动Rp=');
Ws=input('请输入阻带截止频率Ws=');As=input('请输入阻带波动As=');
Wp=Wp*pi;Ws=Ws*pi;
T=1;
Wp1=Wp*T;Ws1=Ws*T;
ep=sqrt(10^(Rp/10)-1);
Ripple=sqrt(1/(1+ep*ep));Attn=10^(-As/20);
[cs,ds]=afd_butt(Wp1,Ws1,Rp,As);
[b,a]=imp_invr(cs,ds,T);
[C,B,A]=dir2par(b,a)
y=ht(b,a,Wp,Ws,Attn,Ripple)
结果：
请输入通带截止频率Wp=0.2
请输入通带波动Rp=1
请输入阻带截止频率Ws=0.3
请输入阻带波动As=15
***Butterworth Filter Order= 6
C =
[]
B =
1.8557 -0.6304
-2.1428 1.1454
0.2871 -0.4466
A =
1.0000 -0.9973 0.2570
1.0000 -1.0691 0.3699
1.0000 -1.2972 0.6949
（2）双线性变换法：
程序：
Wp=input('请输入通带截止频率Wp=');Rp=input('请输入通带波动Rp=');
Ws=input('请输入阻带截止频率Ws=');As=input('请输入阻带波动As=');
Wp=Wp*pi;Ws=Ws*pi;
T=1;Fs=1/T；
Wp1=(2/T)*tan(Wp/2);Ws1=(2/T)*tan(Ws/2);
ep=sqrt(10^(Rp/10)-1);
Ripple=sqrt(1/(1+ep*ep));Attn=10^(-As/20);
[cs,ds]=afd_butt(Wp1,Ws1,Rp,As);
[b,a]=bilinear(cs,ds,T);
[C,B,A]=dir2par(b,a)
y=ht(b,a,Wp,Ws,Attn,Ripple)
结果：
请输入通带截止频率Wp=0.2
请输入通带波动Rp=1
请输入阻带截止频率Ws=0.3
请输入阻带波动As=15
***Butterworth Filter Order= 6
C =
0.0092
B =
1.8767 -0.4881
-2.1781 1.0580
0.2927 -0.4286
A =
1.0000 -0.9459 0.2342
1.0000 -1.0541 0.3753
1.0000 -1.3143 0.7149
相关推荐
小学生劳动心得体会3篇
最好的我们心情说说
巨明：夯实基础务实发展
游泳国家二级运动员标准
C试卷
成都移动营业厅点
纪录片家园观后感5篇精华版
竹里馆相关试题
人教版四年级下册语文教学计划[最新版]
有关爱国卫生月的作文
