人工地震动生成程序

1、clearclccloseallhidden%fni=input(生成人工地震波-输入数据文件名(20041012):,s);fid=fopen(fni,r);fs=fscanf(fid,%f,1);%采样频率%tu=fscanf(fid,%f,1);%上升时间长度%上升时间包络线线形(1-直线、2抛物线、3-指数曲线)iu=fscanf(fid,%f,1);%上升时间包络线线形参数(只有指数曲线需要具体参数，其均为1)cu=fscanf(fid,%f,1);ta=fscanf(fid,%f,1);%持时时间长度td=fscanf(fid,%f,1);%下降时间长度%下降时间包络线线形(1-直

2、线、2抛物线、3-指数曲线)id=fscanf(fid,%f,1);%下降时间包络线线形(只有抛物线，指数曲线需要具体参数，其余为1)cd=fscanf(fid,%f,1);批注MS1:定义强度包络线%dp=fscanf(fid,%f,1);%阴尼比值p=fscanf(fid,%f,1);%概率系数(一般可取P=0.85)nn=fscanf(fid,%f,1);%迭代次数fno=fscanf(fid,%f,1);%输岀数据文件名%对目标反应谱取值x=fscanf(fid,%f,2,inf);%反应谱频率和幅值数据%tatus=fclose(fid);%计算生成地震波的数据长度tl=tu+ta+

3、td;%计算生成地震波的数据长度nt=round(fs*tl+1);%大于并最接近nt的2的幕次方为FFT长度nfft=2Anestpow2(nt)%计算频率间隔(Hz)df=fs/nfft%定义反应谱的离散频率向量f=0:df:(nfft/2-1)*df%计算时间间隔(s)dt=1/fs;%定义的离散时间向量t=0:dt:(nt-1)*dt%生成0到2P啲随机数为随机相位g=rand(1,nfft/s)*2*pi;%建立时间包络线%建立与地震波长度相同元素为1的向量en=ones(1,nt);%上升时间阶段%确定上升时间段的长度l=round(tu*fs)+1%产生上升时间段的包络线数组元素

4、switchiucase1%直线en(1:l)=linspace(0,1,1);%y=linspace(a,b,n)generatesarowvectoryofnpointslinearlyspacedbetweenandincludingaandb.case2%抛物线a=0:l-1;en(1:l)=(a/(l-1).A2;case3%指数曲线a=0:l-1;en(1:l)=1-exp(-cu*a/(l-1);end%持续时间阶段%确定0时刻到持续时间结束时刻时间段的长度m=round(tu+ta)*fs)+1;%下降时间阶段%产生下降时间段的包络线数组元素switchidcase1%直线en

5、(m:nt)=linspace(1,0,nt-m+1);case2%抛物线a=0:nt-m;en(m:nt)=1-cd*(a/(nt-m).A2;case3%指数曲线a=0:nt-m;en(m:nt)=exp(-cd*a/(nt-m);end%按线性插值建立目标反应谱离散数据%按目标反应谱的长度生成元素为0的向量a0=zeros(x(1,:);%取目标反应谱数据的长度n=length(x(1,:);%四舍五入取整求反应谱最大频率对应数组元素的下标nb=round(x(1,n)/df)+1;fork=1:n-1%四舍五入取整求反应谱前一个频率数据对应数组元素的下标l=round(x(1,K)/d

6、f)+1;%四舍五入取整求反应谱后一个频率数据对应数组元素的下标m=round(x(1,K+1)/df)+1;%线性插值产生前后两个频率数据间的反应谱数组元素a0(1:m)=linspace(x(2,k),x(2,k+1),m-l+1)end%根据目标反应谱计算对应的近似功率谱a1=a0;s=zeros(1,nfft/2);k=nb:ne;s(k)=2*dp/(pi.*(a1(k)A2)./f(k”(-2*log(-log(p)*pi/tl)./f(k);%将功率谱转换成傅里叶幅值谱b1=sqrt(4*df*s)*nfft/2;%定义元素为0的反谱传递函数矩阵hf=zeros(ne,nfft)

7、;%计算加速度反应谱传递函数矩阵forj-0:ne-1w=2*pi*df*jwd=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;%通过FFT变换求加速度反应谱传递函数向量hf(j+1,:)=fft(h,nfft);endmm=nn%进行生成人工地震波迭代计算%100为最大迭代次数fork=1:100%将幅值谱和相位谱转化为实部和虚部c=b1.*exp(i*g);%将正负圆

8、频率傅里叶谱向量组合成一仙向量d=c,c(nfft/2:-1:1);%IFFT变换，并取变换结果实部为生成的地震波e=ifft(d,nfft);%给生成的地震波加上强度包络线y=en.*real(e(1:nt);%计算反应谱%对生成的地震波进行FFT变换yf=fft(y,nfft);forj=1:ne%用地震波FFT变换结果和反应谱传递函数的乘积的逆变换做卷积运算d=ifft(yf.*hf(j,:),nfft);%求各频率对应地震的最大响应al(j)=max(real(d(1:nt);end%如果达到指定的迭代次数显于图形ifk=mmsubplot(2,1,1);%m同时显示生成的地震波的强度包络线plot(t,y,t,en,t,-en);xlabel(时间(s);ylabel(加速度(g);gridonsublpot(2,1,2);%同时显示期望反应谱与反应谱计算谱、l=1:neplot(f(l),a0(l),:f(l),a1(l);xlabel(濒率(Hz);ylabel(加速度(g);legend(目标谱,计算谱);gridon;ig=input(继续迭代次数取值1-9,否则退岀:);ifig0&ig10%如果输入数字是1-9mm=mm+igelsebreak;endendc=bl%期望谱与计算谱的比值来修改傅里叶值谱j=nb:ne;bl(j