Matlab实验仿真实验数据

附录

实验二

程序：

（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,l ast_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 =

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.用冲击响应不变法和双线性变换法将模拟滤波器1322

3)(2+++=s s s s H a 变换为数字滤波

器)(z H ，采样周期为s T s 1.0 。

程序：>> 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