冲击响应谱计算的matlab程序

disp('')disp('srs.mver2.0July3,2006')disp('byTomIrvineEmail:[e-mail]')disp('')disp('Thisprogramcalculatestheshockresponsespectrum')disp('ofanaccelerationtimehistory,whichispre-loadedintoMatlab.')disp('Thetimehistorymusthavetwocolumns:time(sec)&acceleration')disp('')%cleart;cleary;clearyy;clearn;clearfn;cleara1;cleara2clearb1;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';elseTHM=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(typically0.05)');elseQ=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(1.-(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)*(1.-2.*(damp^2))*sind-2.*damp*cosd);b3(j)=0;%elseE=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)=1.-Sp;b2(j)=2.*(Sp-C);b3(j)=E^2-Sp;endforward=[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.breakendfn(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,'%7.2f%10.3f%10.3f\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_maxsrs_max=max(abs(x_neg));endsrs_min=min(x_pos);ifmin(abs(x_neg))srs_minsrs_min=min(abs(x_neg));end%figure(1);plot(fn,x_pos,fn,abs(x_neg),'-.');%ifiunit==1ylabel('PeakAccel(G)');elseylabel('PeakAccel(m/sec^2)');endxlabel('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)+0.8));ymin=10^(round(log10(srs_min)-0.6));%fmax=max(fn);fmin=fmax/10.;%fmax=10^(round(log10(fmax)+0.5));%iffn(1)=0.1fmin=0.1;endiffn(1)=1fmin=1;endiffn(1)=10fmin=10;endiffn(1)=100fmin=100;endaxis([fmin,fmax,ymin,ymax]);%disp('')disp('Plotpseudovelocity?');vchoice=input('1=yes2=no');if(vchoice==1)figure(2);%%Converttopseudovelocity%forj=1:jnumifiunit==1x_pos(j)=386.*x_pos(j)/(2.*pi*fn(j));x_neg(j)=386.*x_neg(j)/(2.*pi*fn(j));elsex_pos(j)=x_pos(j)/(2.*pi*fn(j));x_neg(j)=x_neg(j)/(2.*pi*fn(j));endend%srs_max=max(x_pos);ifmax(abs(x_neg))srs_maxsrs_max=max(abs(x_neg));endsrs_min=min(x_pos);ifmin(abs(x_neg))srs_minsrs_min=min(abs(x_neg));end%plot(fn,x_pos,fn,abs(x_neg),'-.');%ifiunit==1ylabel('Velocity(in/sec)');elseylabel('Velocity(m/sec)');endxlabel('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)+0.8));ymin=10^(round(log10(srs_min)-0.6));%axis([fmin,fmax,ymin,ymax]);end