代码.txt

本程序利用金井清功率谱--三角级数法模拟人造地震波 杜修力的是高频和低频部分的修改

clear all
s0=0.6       %% 取为0.64；ωg为 ； 为0.6
wr=2.0,
kecig=0.64    %%场地特征阻尼比，housner and Jennings建议值
omegag=15.6  %%场地特征频率，housner and Jennings建议值
N=20*pi,
d= 0.02*pi,
D=0.03
w0=2*pi*3.5/3/4,
tempw=d:d:N%%频率区间（d～20*pi）
temp=(tempw./omegag).^2;

ad1=(tempw.^4)./( ((tempw.^2)+(w0.^2.)).^2 )
ad2=1./( 1.+( ( D.*tempw).^2 )  )
s0w=(1.+4.*kecig.^2.*temp).*s0./((1.-temp).^2+4.*kecig.^2.*temp).*ad1.*ad2%%K.kanai谱 杜修力修改 
s1w=(1.+4.*kecig.^2.*temp).*s0./((1.-temp).^2+4.*kecig.^2.*temp) %%K.kanai谱

subplot(2,2,1)
plot(tempw,s0w,tempw,s1w)%%画谱
xlabel('freq');
ylabel('S')



for i=1:1:1000
H(i)=chol(s0w(i));%%Cholesky分解
end

thta=2*pi*rand(1000,1000);%%介于0和2pi之间均匀分布的随机数
t=2:2:2000;%%时间区间（0.02～20s）

for j=1:1:1000
a=abs(H);
    b=cos((tempw *j/100) +thta(j));%%因为步长取为0.02，所以j/100
    c=sum(a.*b);
v(j)=2*(d).^(1/2)*c;%%地震荷载模拟公式 
end

subplot(2,2,2)
plot(t/100,v)%%显示地震荷载
xlabel('t(s)');
ylabel('v(t)');

Y=fft(v);%%对数值解作傅立叶变换
Y(1)=[];%%去掉零频量
n=length(Y)/2;%%计算频率个数；
power=abs(Y(1:n)).^2/(length(Y).^2);%%计算功率谱
freq=50*(1:n)/length(Y);%%计算频率，因为步长为0.02，而不是1，故乘以100

subplot(2,2,3)
plot(freq(1:100),power(1:100))%%画功率谱
xlabel('freq');
ylabel('S');

subplot(2,2,4)
plot(freq(1:100),power(1:100),tempw,s0w)%%画功率谱
xlabel('freq');
ylabel('S');