人工地震动生成程序-

1、clearclcclose all hidden %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;%上升时间包络线线形参数(只有指数曲线需要具体参数,其均为1cu=fscanf(fid,'%f',1;ta=fs

2、canf(fid,'%f',1;%持时时间长度td=fscanf(fid,'%f',1;%下降时间长度%下降时间包络线线形(1-直线、2-抛物线、3-指数曲线id= fscanf(fid,'%f',1;%下降时间包络线线形(只有抛物线,指数曲线需要具体参数,其余为1cd=fscanf(fid,'%f',1; dp=fscanf(fid,'%f',1;%阴尼比值p=fscanf(fid,'%f',1;%概率系数(一般可取P=0.85nn=fscanf(fid,'%f',1;%迭代次数f

3、no=fscanf(fid,'%f',1;%输出数据文件名% %对目标反应谱取值x=fscanf(fid,'%f',2,inf;%反应谱频率和幅值数据% tatus=fclose(fid;%计算生成地震波的数据长度tl=tu+ta+td;%计算生成地震波的数据长度nt=round(fs*tl+1;%大于并最接近nt的2的幂次方为FFT长度nfft=2nestpow2(nt%计算频率间隔(Hzdf=fs/nfft%定义反应谱的离散频率向量f=0:df:(nfft/2-1*df%计算时间间隔(sdt=1/fs;%定义的离散时间向量t=0:dt:(nt-1*dt%生成

4、0到2PI的随机数为随机相位g=rand(1,nfft/s*2*pi;%建立时间包络线%建立与地震波长度相同元素为1的向量en=ones(1,nt;%上升时间阶段%确定上升时间段的长度l=round(tu*fs+1%产生上升时间段的包络线数组元素switch iucase 1 %直线en(1:l=linspace(0,1,1;% y = linspace(a,b,n generates a row vector y of n points linearlyspaced between and including a and b.case 2 %抛物线a=0:l-1;en(1:l=(a/(l-1

5、.2;case 3 %指数曲线a=0:l-1;en(1:l=1-exp(-cu*a/(l-1;end%持续时间阶段%确定0时刻到持续时间结束时刻时间段的长度m=round(tu+ta*fs+1;%下降时间阶段%产生下降时间段的包络线数组元素switch idcase 1 %直线en(m:nt=linspace(1,0,nt-m+1;case 2 %抛物线a=0:nt-m;en(m:nt=1-cd*(a/(nt-m.2;case 3 %指数曲线a=0:nt-m;en(m:nt=exp(-cd*a/(nt-m;end%按线性插值建立目标反应谱离散数据%按目标反应谱的长度生成元素为0的向量a0=ze

6、ros(x(1,:;%取目标反应谱数据的长度n=length(x(1,:;%四舍五入取整求反应谱最大频率对应数组元素的下标nb=round(x(1,n/df+1;for k=1:n-1%四舍五入取整求反应谱前一个频率数据对应数组元素的下标 l=round(x(1,K/df+1;%四舍五入取整求反应谱后一个频率数据对应数组元素的下标m=round(x(1,K+1/df+1;%线性插值产生前后两个频率数据间的反应谱数组元素 a0(1:m=linspace(x(2,k,x(2,k+1,m-l+1end%根据目标反应谱计算对应的近似功率谱a1=a0;s=zeros(1,nfft/2;k=nb:ne;s

7、(k=2*dp/(pi.*(a1(k.2./f(k./(-2*log(-log(p*pi/tl./f(k;%将功率谱转换成傅里叶幅值谱b1=sqrt(4*df*s*nfft/2;%定义元素为0的反谱传递函数矩阵hf=zeros(ne,nfft;%计算加速度反应谱传递函数矩阵for j-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;%通过FF

8、T变换求加速度反应谱传递函数向量hf(j+1,:=fft(h,nfft;endmm=nn%进行生成人工地震波迭代计算%100为最大迭代次数for k=1:100%将幅值谱和相位谱转化为实部和虚部c=b1.*exp(i*g;%将正负圆频率傅里叶谱向量组合成一仙向量d=c,c(nfft/2:-1:1;%IFFT变换,并取变换结果实部为生成的地震波e=ifft(d,nfft;%给生成的地震波加上强度包络线y=en.*real(e(1:nt;%计算反应谱%对生成的地震波进行FFT变换yf=fft(y,nfft;for j=1:ne%用地震波FFT变换结果和反应谱传递函数的乘积的逆变换做卷积运算 d=i

9、fft(yf.*hf(j,:,nfft;%求各频率对应地震的最大响应al(j=max(real(d(1:nt;end%如果达到指定的迭代次数显于图形if k=mmsubplot(2,1,1;%m 同时显示生成的地震波的强度包络线plot(t,y,t,en,t,-en;xlabel('时间(s'ylabel('加速度(g'grid onsublpot(2,1,2;%同时显示期望反应谱与反应谱计算谱、l=1:neplot(f(l,a0(l,':'f(l,a1(l;xlabel('频率(Hz'ylabel('加速度(g'legend('目标谱','计算谱'grid on;ig=input('继续迭代次数取值1-9,否则退出:'if ig>0&ig<10 %如果输入数字是1-9mm=mm+igelsebreak;end end c=bl %期望谱与计算谱的比值来修改傅里叶值谱 j=nb:ne; bl(j=c(j.*a0(j.