反应谱 matlab,冲击响应谱计算的matlab程序.doc

冲击响应谱计算的matlab程序.doc

disp( )disp( srs.mver2.0July3,2006 )disp( byTomIrvineEmail:tomirvine@ )disp( )disp( Thisprogramcalculatestheshockresponsespectrum )disp( ofanaccelerationtimehistory,whichispre-loadedintoMatlab. )disp( Thetimehistorymusthavetwocolumns:time(sec)cleary;clearyy;clearn;clearfn;cleara1;cleara2clearb1;clearb2;clearjnum;clearTHM;clearresp;clearx_pos;clearx_neg;%iunit=( Enteraccelerationunit:1=G2=m/sec^2 );%disp( )disp( Selectfile );disp( 1=externalASCIIfile );disp( 2=filepreloadedintoMatlab );file_choice=( );%if(file_choice==1)[filename,pathname]=uigetfile( *.* );filename=fullfile(pathname,filename);%fid=(filename, r );THM=fscanf(fid, %g%g ,[2inf]);THM=THM ;elseTHM=( 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)=( Enterthestartingfrequency(Hz) );iffn(1)>sr/30.fn(1)=sr/30.;end%idamp=( Enterdampingat:1=dampingratio2=Q );%disp( )if(idamp==1)damp=( Enterdampingratio(typically0.05) );elseQ=( Entertheamplificationfactor(typicallyQ=10) );damp=1./(2.*Q);end%disp( )disp( Selectalgorithm: )disp( 1=Kelly-Richman2=Smallwood );ialgorithm=( );%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=( 1=plotonly2=plotdisp( )%ifchoice==2%%[writefname,writepname]=uiputfile( * , SaveSRSdataas );writepfname=fullfile(writepname,writefname);writedata=[fn x_pos (abs(x_neg)) ];fid=(writepfname, w );fprintf(fid, %g%g%g\n ,writedata );fclose(fid);%%%disp( Enteroutputfilename );%SRS_filename=( , s );%%fid=(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))=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=( 1=yes2=no );if(vchoic