人工地震动生成程序

clear

clc

close 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);

%上升时间包络线线形参数（只有指数曲线需要具体参数，其均为1）

cu=fscanf(fid,'%f',1);

ta=fscanf(fid,'%f',1);%持时时间长度

td=fscanf(fid,'%f',1);%下降时间长度

%下降时间包络线线形(1-直线、2-抛物线、3-指数曲线)

id= fscanf(fid,'%f',1);

%下降时间包络线线形（只有抛物线，指数曲线需要具体参数，其余为1）

cd=fscanf(fid,'%f',1);

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+td;

%计算生成地震波的数据长度

nt=round(fs*tl+1);

%大于并最接近nt的2的幂次方为FFT长度

nfft=2^nestpow2(nt)

%计算频率间隔（Hz）

df=fs/nfft

%定义反应谱的离散频率向量

f=0:df:(nfft/2-1)*df

%计算时间间隔（s）

dt=1/fs;

%定义的离散时间向量

t=0:dt:(nt-1)*dt

%生成0到2PI的随机数为随机相位

g=rand(1,nfft/s)*2*pi;

%建立时间包络线

%建立与地震波长度相同元素为1的向量

en=ones(1,nt);

%上升时间阶段

%确定上升时间段的长度

l=round(tu*fs)+1

%产生上升时间段的包络线数组元素

switch iu

case 1 %直线

en(1:l)=linspace(0,1,1);% y = linspace(a,b,n) generates a row vector y of n points linearly

spaced between and including a and b.

case 2 %抛物线

a=0:l-1;

en(1:l)=(a/(l-1)).^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 id

case 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=zeros(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+1)

end

%根据目标反应谱计算对应的近似功率谱

a1=a0;

s=zeros(1,nfft/2);

k=nb:ne;

s(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-1

w=2*pi*df*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;

%通过FFT变换求加速度反应谱传递函数向量

hf(j+1,:)=fft(h,nfft);

end

mm=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=ifft(yf.*hf(j,:),nfft);

%求各频率对应地震的最大响应

al(j)=max(real(d(1:nt)));

end

%如果达到指定的迭代次数显于图形

if k==mm

subplot(2,1,1);

%m 同时显示生成的地震波的强度包络线

plot(t,y,t,en,t,-en);

xlabel('时间（s）');

ylabel('加速度（g）');

grid on

sublpot(2,1,2);

%同时显示期望反应谱与反应谱计算谱、

l=1:ne

plot(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-9

mm=mm+ig

else

break;

end

end

c=bl

%期望谱与计算谱的比值来修改傅里叶值谱

j=nb:ne;

bl(j)=c(j).*a0(j)./a1(j);

end

%打开文件输入人工地震动数据

fid=fopen(fno,'W');

for K=1:nt

%每一行输出两个实型数据,t为时间,y为人工地震动信号值 fprintf(fid,'%f%f\n',t(k),y(k));

end

status=fclose(fid);