单质点阻尼系的地震反应(线性加速度法-MATLAB代码)

function [acc,vel,dis]=sdof_response(h,f,dt,ddy)
%sdof_response 计算单质点阻尼系的地震响应
%   h=阻尼比
%   f=频率(Hz)
%   dt=采样间隔
%   ddy=加速度时程

% h=0.02;
% f=50;
% dt=0.01;
% ddy=rand(1,1000);

n_sample=length(ddy);
acc=zeros(1,n_sample);  %反应加速度时程
vel=zeros(1,n_sample);  %反应速度时程
dis=zeros(1,n_sample);  %反应位移时程

w=2*pi*f;
w2=w*w;
hw=h*w;
wd=w*sqrt(1-h*h);
wdt=wd*dt;
E=exp(-hw*dt);
cwdt=cos(wdt);
swdt=sin(wdt);
a11=E*(cwdt+hw*swdt/wd);
a12=E*swdt/wd;
a21=-E*w2*swdt/wd;
a22=E*(cwdt-hw*swdt/wd);

s1=2*h/w2/w/dt;
s2=(1-2*h*h)/w2/wdt;
s3=1/w2/dt;
s4=h/w/wdt;
b11=E*((1/w2+s1)*cwdt+(h/w/wd-s2)*swdt)-s1;
b12=E*(-s1*cwdt+s2*swdt)-1/w2+s1;
b21=E*(-s3*cwdt-(s4+1/wd)*swdt)+s3;
b22=E*(s3*cwdt+s4*swdt)-s3;

acc(1)=2*hw*ddy(1)*dt;
vel(1)=-ddy(1)*dt;
dis(1)=0;

for i=2:n_sample
    dis(i)=a11*dis(i-1)+a12*vel(i-1)+b11*ddy(i-1)+b12*ddy(i);
    vel(i)=a21*dis(i-1)+a22*vel(i-1)+b21*ddy(i-1)+b22*ddy(i);
    acc(i)=-2*hw*vel(i)-w2*dis(i);
end