两段选取地震波的代码，这里收集在这里，方便后续可能会用到的一些东西：
方法一：NewMark

function Res = Newmark(Acc,DT,T,kesi,IND)
    % parameters
    NT = max(size(T,1),size(T,2));
    NPTS = max(size(Acc,1),size(Acc,2));
    DDY = Acc;
    %% Maxima of input motion
    Emax(1) =max( abs(DDY));
    Emax(2) = 0.0;
    Emax(3) = 0.0;
    DDYF = DDY(1);
    DYF = 0.0;
    YF = 0.0;
    %% Computation of the maximun of a,v,d of the ground motion
%     for i = 2:NPTS
%         DDYM = DDY(i);
%         DY = DYF + (DDYF+DDYM)*DT/2.0;
%         Y = YF + DYF*DT + (DDYF/3.+ DDYM/6.)*DT*DT;
%         Emax(1) = max(Emax(1),abs(DDYM));
%         Emax(2) = max(Emax(2),abs(DY));
%         Emax(3) = max(Emax(3),abs(Y));
%         DDYF  = DDYM;
%         DYF = DY;
%         YF = Y;
%     end
    %Qmax = Emax(IND);%
    %% Response Computation
    H =kesi;
    Res = zeros(1,NT);
    for i = 1:NT
        if (T(i) ==0.0 & IND == 1)
                Res(i) = Emax(1);
        elseif(T(i) ==0.0 &  IND ~= 1)
                Res(i) =0.0;
        else
           W = pi*2/T(i);
           W2 = W*W;
           HW = H*W;
           WD = W*sqrt(1-H.^2);
           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);
           SS = -HW*SWDT - WD*CWDT;
           CC = -HW*CWDT + WD*SWDT;
           S1 = (E*SS + WD)/W2;
           C1 = (E*CC + HW)/W2;
           S2 = (E*DT*SS + HW*S1 + WD*C1)/W2;
           C2 = (E*DT*CC + HW*C1 - WD*S1)/W2;
           S3 = DT*S1-S2;
           C3 = DT*C1-C2;
           B11 = -S2/WDT;
           B12 = -S3/WDT;
           B21 = (HW*S2 - WD*C2)/WDT;
           B22 = (HW*S3 - WD*C3)/WDT;
           Rmax(1) = 2.*HW*abs(DDY(1))*DT;
           Rmax(2) = abs(DDY(1))*DT;
           Rmax(3) = 0.0;
           DXF = -DDY(1)*DT;
           XF = 0.0;
           for k = 2:NPTS
               DDYM = DDY(k);
               DDYF = DDY(k-1);
               X = A12*DXF + A11*XF + B12*DDYM + B11*DDYF;
               DX = A22*DXF + A21*XF + B22*DDYM + B21*DDYF;
               DDX = -2.*HW*DX - W2*X;
               Rmax(1) = max(Rmax(1),abs(DDX));
               Rmax(2) = max(Rmax(2),abs(DX));
               Rmax(3) = max(Rmax(3),abs(X));
               DXF = DX;
               XF = X;
           end
           Res(i)=Rmax(IND);
        end
    end
end
%% Newmark-Method To calculate the Response-SPEC of ground motion;
%   Acc: the Ground Motion To be calculate;
%   DT:  the time increment of the ground motion
%   T:   the Perids of the Response_SPEC
%   kesi:damping_ratio
% IND: 1 Sa; 2 Sv; 3 Sd

方法二：谱加速度法

clear;
clc;
S1= [];
dirs=dir('F:\1study\2earthquake\1seismic\*.AT2');
dircell=struct2cell(dirs)' ; 
filenames=dircell(:,1);
filenames1=strcat('F:\1study\2earthquake\1seismic\',filenames);
n=length(filenames1);
PGA=0.07;  %%目标PGA
%%
% Tt1=0.831;   %%结构基本周期
% Tt2=0.263;   %%结第二周期
% Tt3=0.146;   %%结构第三周期
%%
Tt1=1.578;   %%结构基本周期
Tt2=0.499;   %%结第二周期
Tt3=0.275;   %%结构第三周期
%%
S=zeros(301,1);
T=zeros(301,1);
abc=zeros(n,5);
%% fname为输入文件名，从PEER网站下载的地震加速度记录
for j=1:n
fname=char(filenames1(j));
fin=fopen(fname,'r');
  for m=1:4
    str = fgetl(fin);%the first 3 loops are meant be skiped.
  end
Par=textscan(str, '%*s %f %*s %f %*s','delimiter','=');
NPTS=Par{1};DT=Par{2};
abc(j,2)=NPTS;
abc(j,3)=DT;
LINES = (floor(NPTS/5)+1);
Acc = [];
    for k = 1:LINES
     str = fgetl(fin);
         if(str ~= -1) % 判断读取结果是否有效
         P = sscanf(str,'%f');
         Acc = [Acc;P];
         end
    end    
[aa,bb]=max(abs(Acc));
abc(j,4)=aa;     %%每条地震波原始的PGA
abc(j,5)=PGA/aa;
Acc=PGA/aa*Acc;  %%每条地震波PGA调幅
%%
     for l=0:1:300;  %%0-6s共301个点
        i=l*0.02;
        T(l+1,1) = i;
        S(l+1,1) = Newmark(Acc,DT,i,0.03,1); 
     end
 S1= [S1;S];
sta=fclose(fin);
end
S2=reshape(S1,301,n);   %%每条地震波的加速度反应谱
%%
aaa=num2cell(abc);
for q=1:n
    aaa(q,1)=filenames(q);   %%输入地震波的信息
end
%%
fname = 'F:\1study\2earthquake\spectrum0.16.txt';  %%读取设计反应谱
fin=fopen(fname,'r');
A1 = importdata(fname);
sta=fclose(fin);
%%
%%控制段
T0=0.1;
Tg=0.35;
deltaT1=0.2;
deltaT2=0.5;
kz1=round((Tg-T0)/0.02)+1;
kz2=round((deltaT1+deltaT2)/0.02)+1;
kz3=round((Tt1-0.2)/0.02)+1;
kz=kz1+kz2;
rsme=zeros(1,n);
for i=1:n
    kzz1=0;
    kzz2=0;
    for j=6:(kz1+5)
        kzz1=kzz1+(S2(j,i)-A1(j))^2;
    end
    for j=kz3:(kz3+kz2-1)
        kzz2=kzz2+(S2(j,i)-A1(j))^2;
    end
    rsme(i)=sqrt((kzz1+kzz2)/kz);         %%均方根误差
end
%%
[rsme1,index1]=sort(rsme);        %%从小到大排序
ts=10;              %%所选地震波条数
output=cell(ts+1,6);     %%最终选波结果
output(1,1)=cellstr('地震名称');
output(1,2)=cellstr('NPTS值');
output(1,3)=cellstr('NT值');
output(1,4)=cellstr('原PGA值');
output(1,5)=cellstr('调幅因子');
output(1,6)=cellstr('RMSE值');
%%
rsme2=num2cell(rsme);
S3=zeros(301,ts);
for i=1:ts
    S3(:,i)=S2(:,index1(i));
    output(i+1,1)=aaa(index1(i),1);   %%地震名称
    output(i+1,2)=aaa(index1(i),2);   %%NPTS值
    output(i+1,3)=aaa(index1(i),3);   %%NT值
    output(i+1,4)=aaa(index1(i),4);   %%原PGA值
    output(i+1,5)=aaa(index1(i),5);   %%调幅因子
    output(i+1,6)=rsme2(index1(i));   %%RMSE值
end
%%
% plot(T,S2)  %%画所有反应谱曲线
plot(T,S3)  %%画选择后的反应谱曲线
hold on
plot(T,A1,'k','LineWidth',2)  %%画设计反应谱曲线