matlab程序函数Word格式.docx
《matlab程序函数Word格式.docx》由会员分享，可在线阅读，更多相关《matlab程序函数Word格式.docx（38页珍藏版）》请在冰豆网上搜索。
%theFouriertransformX（w）andthefrequencyvectorw.
[n,m]=size（x）;
ifn<
m,
x=x'
Xk=fft（x）;
N=length（x）;
k=0:
k
（1）=eps;
X=（1-exp（-j*2*pi*k/N））./（j*2*pi*k/N/T）.*Xk'
w=2*pi*k/N/T;
4.hybrid.mSimulationofcontinuous-timesystemanddigitalcontroller
function[Y,UD]=hybrid（Np,Dp,Nd,Dd,T,t,U）;
%HYBRID[Y,UD]=hybrid（Np,Dp,Nd,Dd,T,t,U）;
%HYBRIDreturnstheoutputvectorYandthecontrolvectorUC
%ofahybridsystemwhichisdefinedbyNp,Dp,NdandDd.
%Thetransferfunctionofthecontinuoustimeplantisstored
%inNpandDp.NdandDdarethenumeratoranddenominators
%ofthetransferfunctionofadiscretetimecontroller.
%Tisthesamplingperiodofthedigitalcontroller.
%tisarowvectorofthetimesatwhichtheinputsinUappear
%andUisarowvectorofthesysteminput.
%TheincrementintmustbechosensothatTisanintegermultiple
%oftheincrement,i.e.,t=a:
b:
cthenb=T/mwherem>
1isaninteger.
%Notethatthesystemhasonlyoneinput.
[Ac,Bc,Cc,Dc]=tf2ss（Np,Dp）;
[Ad,Bd,Cd,Dd]=tf2ss（Nd,Dd）;
nsam=T/（t
（2）-t
（1））;
%numberofintegrationptspersample
%initialize
Y=0;
UD=0;
[ncr,ncc]=size（Ac）;
xc0=zeros（ncr,1）;
[ndr,ndc]=size（Ad）;
xdk=zeros（ndr,1）;
kmax=fix（t（length（t））/T）;
%#ofcompletesamplesint
fork=0:
kmax-1
%calculatecontrolandoutputofzoh
ek=U（k*nsam+1）-Y（k*nsam+1）;
xd=Ad*xdk+Bd*ek;
zoh=Cd*xdk+Dd*ek;
xdk=xd;
%integratecontinuous-timeplantwithinputofzohforTseconds
udi=zoh*ones（nsam+1,1）;
ti=t（k*nsam+1:
（k+1）*nsam+1）;
tint=t（1:
nsam+1）;
%requiredforMatlab5.0,shiftstime,but
%sincethisisatime-invariantsystem,theresultswon'
tchange
[yi,xi]=lsim（Ac,Bc,Cc,Dc,udi,tint,xc0）;
xc0=xi（nsam+1,:
）;
%augmentvectors
Y=[Y;
yi（2:
nsam+1）];
UD=[UD;
udi（2:
if（kmax*nsam+1<
length（t））
%computetailofsimulationfromt（kmax*nsam）tot_end
k=kmax;
%calculatecontrolandoutputofzoh
%integratecontinuous-timeplantwithinputofzoh
ti=t（k*nsam+1:
length（t））;
tint=ti-ti
（1）;
%again,neededforMATLAB5.0
udi=zoh*ones（length（ti）,1）;
%augmentvectors
length（yi））];
length（udi））];
5.recur.mComputessolutionofdifferenceequations
functiony=recur（a,b,n,x,x0,y0）;
%y=recur（a,b,n,x,x0,y0）
%solvesfory[n]from:
%y[n]+a1*y[n-1]+a2*y[n-2]...+an*y[n-N]
%=b0*x[n]+b1*x[n-1]+...+bm*x[n-M]
%a,b,n,x,x0andy0arevectors
%a=[a1a2...aN]
%b=[b0b1...bM]
%ncontainsthetimevaluesforwhichthesolutionwillbecomputed
%y0containstheinitialconditionsfory,inorder,
%i.e.,y0=[y[n0-N],y[n0-N+1],...,y[n0-1]]
%wheren0representsthefirstelementofn
%x0containstheinitialconditionsonx,inorder
%i.e.,x0=[x[n0-M],...,x[n0-1]]
%theoutput,y,haslength（n）
N=length（a）;
M=length（b）-1;
iflength（y0）~=N,
error（'
Lengthsofaandy0mustmatch'
）
iflength（x0）~=M,
Lengthofx0mustmatchlengthofb-1'
y=[y0zeros（1,length（n））];
x=[x0x]
a1=a（length（a）:
-1:
1）%reversestheelementsina
b1=b（length（b）:
1）
fori=N+1:
N+length（n）,
y（i）=-a1*y（i-N:
i-1）'
+b1*x（i-N:
i-N+M）'
y=y（N+1:
N+length（n））
以下是一些书上的程序（用空行隔开）
%Example1.3
%readsdatafileandplotsclosingpriceofQQQQ
c=csvread（'
QQQQdata1.csv'
1,4,[14104]）;
n=1:
10;
plot（n,c（n）,n,c（n）,'
o'
grid
xlabel（'
Day（n）'
ylabel（'
ClosingPrice'
title（'
Figure1.22'
%Example1.4
%closingpricec[n]plottedusingo'
s
%11-dayMAfilteroutputy[n]plottedusing*'
QQQQdata2.csv'
1,4,[14504]）;
fori=11:
50;
y（i）=（1/11）*sum（c（i-10:
i））;
end;
n=11:
n,y（n）,n,y（n）,'
*'
c[n]andy[n]'
Figure1.27'
pause
%plotofc[n]andy[n+5]
n=6:
45;
n,y（n+5）,n,y（n+5）,'
axis（[5453437.5]）
c[n]andy[n+5]'
Figure1.28'
%Figure1.10
%givesexampleofacontinuous-timesignal
t=0:
0.1:
30;
x=exp（-.1*t）.*sin（2/3*t）;
plot（t,x）
axis（[030-11]）;
x（t）'
Time（sec）'
Figure1.10'
%Figure1.12
n=-2:
6;
x=[001210-100];
stem（n,x,'
filled'
n'
x[n]'
Figure1.12'
%Figure1.14
n=0:
x=exp（-.1*n）.*sin（2/3*n）;
Figure1.14'
%ThisplotsFigure1.18a）andb）
n=-10:
OMEGA1=pi/3;
OMEGA2=1;
x1=cos（OMEGA1*n）;
x2=cos（OMEGA2*n）;
clf
subplot（221）,stem（n,x1,'
Figure1.18a'
axis（[-10,30,-1.5,1.5]）;
subplot（222）,stem（n,x2,'
Figure1.18b'
subplot（111）
%Example2.4
p=[0ones（1,10）zeros（1,5）];
%correspondston=-1ton=14
x=p;
v=p;
y=conv（x,v）;
n=-2:
25;
stem（n,y（1:
length（n）））%seethetutorialforobtainingclosedcircles
Example2.4'
x[n]*v[n]'
%Example2.5
40;
x=sin（.2*n）;
h=sin（.5*n）;
y=conv（x,h）;
stem（n,h）;
%seetutorialforinstructionsonclosedcircles
Example2.5'
h[n]'
stem（n,x）
length（n）））
y[n]'
%Example2.11
%numericalsolutionofRCcircuit
R=1;
C=1;
T=.2;
a=-（1-T/R/C）;
b=[0T/C];
y0=0;
x0=1;
n=1:
x=ones（1,length（n））;
y1=recur（a,b,n,x,x0,y0）;
%approximatesolution
%comparetoexactsolution
t=0:
.04:
8;
y2=1-exp（-t）;
%exactsolution
y1=[y0y1];
%augmentsI.C.ontovectorforplotting
%redefinesnaccordingly
plot（n*T,y1,'
t,y2,'
-'
%thefollowinginsertsalegend
holdon
plot（[22.32.6],[.6.6.6],'
[22.32.6],[.5.5.5]）
Example2.11'
y（t）'
text（3,.6,'
Approximatesolution'
text（3,.5,'
Exactsolution'
H=[];
h=[];
pubplot
print-dill\book\plotting\fig2_7
holdoff
%Example2.12,useoftheODEsolver
tspan=[08];
%vectorofinitialandfinaltimes
%initialvaluefory（t）
[t,y]=ode45（@ex2_12_func,tspan,y0）;
%ComparetoEulerapproximationfromExample2.11
b=[0T/R/C];
t2=0:
y2=1-exp（-t2）;
t,y,'
--'
t2,y2,'
legend（'
EulerApproximation'
'
Runge-KuttaApproximation'
ExactSolution'
%Example2.16
%Convolutionofexponentials
symstlambday
y=int（exp（2*lambda-t）,lambda,0,t）;
%for0<
t<
1
simplify（y）
y=int（exp（2*lambda-t）,lambda,0,1）+int（exp（2-t）,lambda,1,t）;
%for1<
2
y=int（exp（2*lambda-t）,lambda,0,1）+int（exp（2-t）,lambda,1,2）;
%for2<
4
y=int（exp（2*lambda-t）,lambda,t-4,1）+int（exp（2-t）,lambda,1,2）;
%for4<
5
y=int（exp（2-t）,lambda,t-4,2）%for5<
6
%Figure2.7
%LoanBalanceprogram
%Programcomputesloanbalancey[k]
y0=input（'
Amountofloan'
I=input（'
YearlyInterestrate'
c=input（'
Monthlyloanpayment'
%x[n]=c
y=[];
%definesyasanemptyvector
y
（1）=（1+（I/12））*y0-c;
forn=2:
360
y（n）=（1+（I/12））*y（n-1）-c;
ify（n）<
0,break,end
%Thefollowingcommandsarefordisplayingtheresults
formatbank
length（y）;
i=1;
fprintf（'
\nny[n]indollars'
）
whilei<
=length（n）
m=min（[19length（n）-i]）;
out=[n（i:
i+m）'
y（i:
]
i=i+20;
（Hitreturntocontinue）'
pause
formatshorte
y=[y0zeros（1,len
升级会员 