冲击响应谱计算的matlab程序.docx
《冲击响应谱计算的matlab程序.docx》由会员分享，可在线阅读，更多相关《冲击响应谱计算的matlab程序.docx（10页珍藏版）》请在冰豆网上搜索。
冲击响应谱计算的matlab程序
disp（''）
disp（'verJuly3,2006'）
disp（'byTomIrvineEmail'）
disp（''）
disp（'Thisprogramcalculatestheshockresponsespectrum'）
disp（'ofanaccelerationtimehistory,whichispre-loadedintoMatlab.'）
disp（'Thetimehistorymusthavetwocolumns:
time（sec）&acceleration'）
disp（''）
%
cleart;
cleary;
clearyy;
clearn;
clearfn;
cleara1;
cleara2
clearb1;
clearb2;
clearjnum;
clearTHM;
clearresp;
clearx_pos;
clearx_neg;
%
iunit=input（'Enteraccelerationunit:
1=G2=m/sec^2'）;
%
disp（''）
disp（'Selectfileinputmethod'）;
disp（'1=externalASCIIfile'）;
disp（'2=filepreloadedintoMatlab'）;
file_choice=input（''）;
%
if（file_choice==1）
[filename,pathname]=uigetfile（'*.*'）;
filename=fullfile（pathname,filename）;
%
fid=fopen（filename,'r'）;
THM=fscanf（fid,'%g%g',[2inf]）;
THM=THM';
else
THM=input（'Enterthematrixname:
'）;
end
%
t=double（THM（:
1））;
y=double（THM（:
2））;
%
tmx=max（t）;
tmi=min（t）;
n=length（y）;
%
out1=sprintf（'\n%dsamples\n',n）;
disp（out1）
%
dt=（tmx-tmi）/（n-1）;
sr=1./dt;
%
out1=sprintf（'SR=%gsamples/secdt=%gsec\n',sr,dt）;
disp（out1）
%
fn
（1）=input（'Enterthestartingfrequency（Hz）'）;
iffn
（1）>sr/30.
fn
（1）=sr/30.;
end
%
idamp=input（'Enterdampingformat:
1=dampingratio2=Q'）;
%
disp（''）
if（idamp==1）
damp=input（'Enterdampingratio（typically'）;
else
Q=input（'Entertheamplificationfactor（typicallyQ=10）'）;
damp=1./（2.*Q）;
end
%
disp（''）
disp（'Selectalgorithm:
'）
disp（'1=Kelly-Richman2=Smallwood'）;
ialgorithm=input（''）;
%
tmax=（tmx-tmi）+1./fn
（1）;
limit=round（tmax/dt）;
n=limit;
yy=zeros（1,limit）;
fori=1:
length（y）
yy（i）=y（i）;
end
%
disp（''）
disp（'Calculatingresponse.....'）
%
%SRSengine
%
forj=1:
1000
%
omega=2.*pi*fn（j）;
omegad=omega*sqrt（damp^2））;
cosd=cos（omegad*dt）;
sind=sin（omegad*dt）;
domegadt=damp*omega*dt;
%
if（ialgorithm==1）
a1（j）=2.*exp（-domegadt）*cosd;
a2（j）=-exp（-2.*domegadt）;
b1（j）=2.*domegadt;
b2（j）=omega*dt*exp（-domegadt）;
b2（j）=b2（j）*（（omega/omegad）*.*（damp^2））*sind-2.*damp*cosd）;
b3（j）=0;
%
else
E=exp（-damp*omega*dt）;
K=omegad*dt;
C=E*cos（K）;
S=E*sin（K）;
Sp=S/K;
%
a1（j）=2*C;
a2（j）=-E^2;
b1（j）=;
b2（j）=2.*（Sp-C）;
b3（j）=E^2-Sp;
end
forward=[b1（j）,b2（j）,b3（j）];
back=[1,-a1（j）,-a2（j）];
%
resp=filter（forward,back,yy）;
%
x_pos（j）=max（resp）;
x_neg（j）=min（resp）;
%
jnum=j;
iffn（j）>sr/8.
break
end
fn（j+1）=fn
（1）*（2.^（j*（1./12.）））;
end
%
%Outputoptions
%
disp（''）
disp（'Selectoutputoption'）;
choice=input（'1=plotonly2=plot&outputtextfile'）;
disp（''）
%
ifchoice==2
%%
[writefname,writepname]=uiputfile（'*','SaveSRSdataas'）;
writepfname=fullfile（writepname,writefname）;
writedata=[fn'x_pos'（abs（x_neg））'];
fid=fopen（writepfname,'w'）;
fprintf（fid,'%g%g%g\n',writedata'）;
fclose（fid）;
%%
%disp（'Enteroutputfilename'）;
%SRS_filename=input（'','s'）;
%
%fid=fopen（SRS_filename,'w'）;
%forj=1:
jnum
%fprintf（fid,'%%%\n',fn（j）,x_pos（j）,abs（x_neg（j）））;
%end
%fclose（fid）;
end
%
%PlotSRS
%
disp（''）
disp（'Plottingoutput.....'）
%
%Findlimitsforplot
%
srs_max=max（x_pos）;
ifmax（abs（x_neg））>srs_max
srs_max=max（abs（x_neg））;
end
srs_min=min（x_pos）;
ifmin（abs（x_neg））srs_min=min（abs（x_neg））;
end
%
figure
（1）;
plot（fn,x_pos,fn,abs（x_neg）,'-.'）;
%
ifiunit==1
ylabel（'PeakAccel（G）'）;
else
ylabel（'PeakAccel（m/sec^2）'）;
end
xlabel（'NaturalFrequency（Hz）'）;
Q=1./（2.*damp）;
out5=sprintf（'AccelerationShockResponseSpectrumQ=%g',Q）;
title（out5）;
grid;
set（gca,'MinorGridLineStyle','none','GridLineStyle',':
','XScale','log','YScale','log'）;
legend（'positive','negative',2）;
%
ymax=10^（round（log10（srs_max）+）;
ymin=10^（round（log10（srs_min））;
%
fmax=max（fn）;
fmin=fmax/10.;
%
fmax=10^（round（log10（fmax）+）;
%
iffn
（1）>=
fmin=;
end
iffn
（1）>=1
fmin=1;
end
iffn
（1）>=10
fmin=10;
end
iffn
（1）>=100
fmin=100;
end
axis（[fmin,fmax,ymin,ymax]）;
%
disp（''）
disp（'Plotpseudovelocity'）;
vchoice=input（'1=yes2=no'）;
if（vchoice==1）
figure
（2）;
%
%Converttopseudovelocity
%
forj=1:
jnum
ifiunit==1
x_pos（j）=386.*x_pos（j）/（2.*pi*fn（j））;
x_neg（j）=386.*x_neg（j）/（2.*pi*fn（j））;
else
x_pos（j）=x_pos（j）/（2.*pi*fn（j））;
x_neg（j）=x_neg（j）/（2.*pi*fn（j））;
end
end
%
srs_max=max（x_pos）;
ifmax（abs（x_neg））>srs_max
srs_max=max（abs（x_neg））;
end
srs_min=min（x_pos）;
ifmin（abs（x_neg））srs_min=min（abs（x_neg））;
end
%
plot（fn,x_pos,fn,abs（x_neg）,'-.'）;
%
ifiunit==1
ylabel（'Velocity（in/sec）'）;
else
ylabel（'Velocity（m/sec）'）;
end
xlabel（'NaturalFrequency（Hz）'）;
Q=1./（2.*damp）;
out5=sprintf（'PseudoVelocityShockResponseSpectrumQ=%g',Q）;
title（out5）;
grid;
set（gca,'MinorGridLineStyle','none','GridLineStyle',':
','XScale','log','YScale','log'）;
legend（'positive','negative',2）;
%
ymax=10^（round（log10（srs_max）+）;
ymin=10^（round（log10（srs_min））;
%
axis（[fmin,fmax,ymin,ymax]）;
end
升级会员 