直接法计算冲击响应谱

%% 直接积分法：利用Duhamel积分来求解单自由度系统响应谱的方法

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear   %清楚内存中所有变量和函数

clc     %清楚工作窗口中所有显示的内容

close all hidden %关闭所有隐藏的窗口

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 激励数据读取

fki=input('激励-输入文件名：','s');

fid=fopen(fki,'r');%以只读方式打开数据文件

dt=input('采样时间间隔dt：');%fs=fscanf(fid,'%f',1);%采样频率（200Hz）

f0=input('分析频率范围f0：');%fs=fscanf(fid,'%f',1);%采样频率（200Hz）

ddx=fscanf(fid,'%f',[1,inf]);%读入时程数据存成行向量

total=zeros(1,length(ddx));

tau_max=zeros(1,f0);

tau_min=1000+zeros(1,f0);

xi=0.05;

for f=1:f0

    wn=2*pi*f;

    wd=wn*sqrt(1-xi^2);

    for i=0:length(ddx)-1

        su=zeros(1,length(ddx));

        for k=0:i

            e=exp(-xi*wn*dt*(i-k));

            s=sin(wd*(i-k)*dt);

            su(k+1)=ddx(k+1)*e*s*dt/wd;

        end       

        total(i+1)=sum(su);

    end

    tau_max(f)=wn*max(total); 

    tau_min(f)=wn*min(total);   

end

f=1:f0;

plot((1:f0),tau_max,'--o',(1:f0),tau_min); 

hold on;

fno=input('请输入保存文件名称及后缀：','s');

fid=fopen(fno,'w');%以写的方式打开文件或建立一个新文件

for k=1:f0

    fprintf(fid,'%f %f \n',f(k),tau_max(k));

end