earthquake_1.rar_人工 地震波_人工地震_地震_地震波 合成_地震波生成

clear clc close all hidden %fni=input('生成人工地震波-输入数据文件名:','s'); %fid=fopen(fni,'r'); %地震随机波的生成 eq=xlsread('C:\Users\dongjia\Desktop\matlab程序\地震波生成程序\jcsj.xlsx'); fs=eq(1,1); %fscanf(fid,'%f',1);%采样频率 tu=eq(2,1);%fscanf(fid,'%f',1);%上升时间长度 %iu=fscanf(fid,'%f',1);%上升时间包络线线型（1-直线，2-抛物线，3-指数曲线） %cu=fscanf(fid,'%f',1);%上升时间包络线线型参数(指数曲线需要相应参数，其他为1） ta=eq(3,1);%fscanf(fid,'%f',1);%持续时间长度 td=eq(4,1);%fscanf(fid,'%f',1);%下降时间长度 %id=fscanf(fid,'%f',1);%下降时间包络线线型（1-直线，2-抛物线，3-指数曲线） cd=eq(5,1);%fscanf(fid,'%f',1);%下降时间包络线线型参数,与地震波时间有关，5:1.5；10:1.15；20:0.8；30:0.64； dp=eq(6,1);%fscanf(fid,'%f',1);%阻尼比,根据工程类型确定，水闸一般7%，土石坝20%，拱坝5%，重力坝10% p=eq(7,1);%fscanf(fid,'%f',1);%概率系数（一般0.85以上） nn=eq(8,1);%fscanf(fid,'%f',1);%迭代次数 amax=eq(9,1);%最大加速度值 %fno=fscanf(fid,'%f',1);%输出数据文件名 x=xlsread('C:\Users\dongjia\Desktop\matlab程序\地震波生成程序\mbp.xlsx');%fscanf(fid,'%s',[2,inf]);%反应谱频率和幅值数据 %tatus=fclose(fid); tl=tu+ta+td;%计算地震波总时间长度 nt=fs*tl+1;%round(fs*tl+1);%计算地震波的数据长度 nfft=2^nextpow2(nt);%大于或最接近n的2的幂次方为FFT长度 df=fs/nfft;%频率间隔 %f=0:df:(nfft/2-1)*df;%反应谱的离散频率向量 频率点数 f=df:df:nfft*df; T=1./f; nf=length(f); dt=1/fs;%计算时间间隔 t=0:dt:(nt-1)*dt;%离散时间向量 g=rand(nf,1)*2*pi;%随机数的随机相位 相位点数 %建立时间包络线 %建立与地震波长度相同的单元为1的向量 en=ones(nt,1); Tg=0.4;%场地确定的地震时间 bmax=2.25;%反应谱最大代表值 %确定时间过程中的包络线 for v=1:1:nt if t(v)>=0&&t(v)<tu en(v,1)=(t(v)/tu)^2; elseif t(v)>=tu&&t(v)<tu+ta en(v,1)=1; elseif t(v)>=tu+ta %t0(v)<=t3 en(v,1)=exp((-1)*cd*(t(v)-tu-ta)); end end %规范标准反应谱 beta=zeros(nf,1); for m=1:1:nf if T(m)>0&&T(m)<0.1 beta(m)=1.0+12.5*T(m); elseif T(m)>=0.1&&T(m)<Tg beta(m)=bmax; elseif T(m)>=Tg beta(m)=bmax*(Tg/T(m))^0.6; %if beta(m)<bmax*0.2 %beta(m)=bmax*0.2; end end figure(1) subplot(2,1,1); plot(t,en(:,1),'k'); xlabel('时间/s'); ylabel('包络值'); subplot(2,1,2); plot(T,beta(:,1),'k'); xlabel('周期/s'); xlim([0,3]); ylabel('反应值代表值'); %按线性插值建立目标反应谱离散数据 %按目标反应谱长度生成元素为0的向量 %绘制标准反应谱 %t0=0.1; %Tg=0.4;%根据规范及场地类型确定 %beta=2.25;%标准反应谱设计值最大代表值 a0=zeros(nf,1); %取目标反应谱数据的长度 n=length(x(:,1)); %求反应谱最小频率对应数组元素的下标 nb=round(x(1,1)/df)+1; %求反应谱最大频率对应数组元素的下标 ne=round(x(n,1)/df)+1; for k=1:1:nf a0(k,1)=beta(k,1); end %for k=1:n-1 %l=round(x(k,1)/df)+1; %反应谱前一个频率对应数组下标 %m=round(x(k+1,1)/df)+1;%反应谱下一个频率对应数组下标 %a0(l:m)=linspace(x(k,2),x(k+1,2),m-l+1);%两个频率数据之间的反应谱数组元素 %end %根据目标反应谱计算对应的近似功率谱 a1=amax*a0; s=zeros(nf,1); c0k=zeros(nf,1); for j=1:1:nf s(j,1)=dp/(pi*2*pi*f(j))*a1(j)^2/log(-pi/(2*pi*f(j)*tl)*log(p)); c0k(j,1)=sqrt(4*s(j)*2*pi*df)*nfft; end %k=nb:ne; %s(k)=dp/pi.*(al(k).^2)./(2*pi*f(k))./(-2*log(-log(p)*pi/tl)./f(k)); %将功率谱转化为傅里叶幅值谱 %bl=sqrt(4*df*s)*nfft/2; %定义元素全为0的反应谱传递函数矩阵 %hf=zeros(nf,nfft); %计算加速度反应谱传递函数矩阵 %r=zeros(nf,1); %for j=1:nf %w=2*pi*f(j); % wd=w*sqrt(1-dp*dp); %e=exp(-t.*w*dp); % a=t.*wd; % s=sin(a).*((1-2*dp*dp)/(1-dp*dp)); %c=cos(a).*(2*dp/sqrt(1-dp*dp)); % h=wd*e.*(s+c)/fs;%计算加速度反应谱脉冲响应函数向量 % r(j)=max(abs(conv(at,h))); % hf(j,:)=fft(h,nfft);%通过FFT变换求加速度反应谱传递函数向量 %end %mm=nn; %进行生成人工地震波迭代计算 %最大迭代次数100 for k=1:1 c=c0k.*exp(1i*g);%将幅值谱和相位谱转换为实部和虚部 %d=[c,c(nfft/2:-1:1)];%正负圆频率傅里叶谱向量合成一个向量 e=ifft(c,nfft);%ifft逆变换，并取变换结果生成的实部为生成的地震波 b(:,1)=e(1:nt); y=en.*real(b);%给生成的地震波加上包络线 r=zeros(nf,1); for j=1:nf w=2*pi*f(j); wd=w*sqrt(1-dp*dp); e=exp(-t.*w*dp); a=t.*wd; s=sin(a).*((1-2*dp*dp)/(1-dp*dp)); c=cos(a).*(2*dp/sqrt(1-dp*dp)); h=wd*e.*(s+c)/fs;%计算加速度反应谱脉冲响应函数向量 r(j)=max(abs(conv(y,h)));%计算反应谱 sum(j,:)=conv(y,h); % hf(j,:)=fft(h,nfft);%通过FFT变换求加速度反应谱传递函数向量 end Amax=zeros(nf,2); for i=1:1:nf [ma,I]=max(abs(sum(i,:))); Amax(i,1)=I; Amax(i,2)=sum(i,I); end rx=zeros(n,1);%频率控制点对应的计算反应值 res=zeros(n,1);%频率控制点对应的计算反应值与目标值之差 q=zeros(n,1);%频率控制点对应的计算反应值与目标值之比 xnum=zeros(n,1);%频率控制点对应的频率数据点中下限 for i=1:1:n rx(i,1)=(r(fix(x(i,2)/df))+r(fix(x(i,2)/df)+1))/2; res(i,1)=abs((x(i,3)-rx(i,1))/x(i,3)); q(i)=x(i,3)/rx(i,1); xnum(i)=fix(x(i,2)/df); end if max(res)<0.1 break; else c=zeros(nf,1); for i=1:1:nf if i<=xnum(1) c(i,1)=q(1,1); else for j=1:1:n-1 if i>xnum(j)&i<=xnum(j+1) c(i,1)=(q(j+1,1)-q(j,1))/(xnum(j+1)-xnum(j))*(i-xnum(j))+q(j,1); end end end end c0kT=c0k; j=1:nf; c0k(j)=c(j).*c0kT;%a0(j)./a1(j); end end %mm=nn; %计算反应谱 %对生成的地震波进行FFT变换 %yf=fft(y,nfft); %yf1=yf'; %for j=1:nf % d=ifft(yf1.*hf(j,:),nfft);%逆变换进行卷积运算 % d1(:,1)=d(1:nt); %a1(j)=max(real(d1));%求各频率对应地震的最大响应 %end %如果达到指定迭代次数显示图形 figure(2) subplot(2,1,1); %同时显示生成的地震波和强度包络线 plot(t,y,t,en,t,-en); xlabel('时间（s)'); ylabel('加速度（m/s^2)'); subplot(2,1,2); %显示期望反应谱与反应计算谱 %l=1:nf; plot(1./f,a1,1./f,r); xlabel('频率（Hz）'); ylabel('加速度（m/s^2）'); legend('目标普’,’计算谱');