多点人工生成地震波程序
我编了一个多支承点拟合地震波的程序，可是计算结果很大，有些不符合实际，程序如下，请高手指点一下，谢谢！clc; clear all
% history time
for ii=1:2000
t=ii*0.02; g(ii)=ii*0.02;
if t>= 0 & t <= 2, yy(ii)=(t/2)^2;
elseift>2 & t<= 10, yy(ii)=1.0;
else yy(ii)=exp(-0.115*(t-10)); % 形函数
end
end
for kkkk=1:2000
t=kkkk*0.02;
sum1=0.0; sum2=0.0; sum3=0.0; sum4=0.0;
for jjj=1:1000
u1(jjj)=0.0; u2(jjj)=0.0; u3(jjj)=0.0; u4(jjj)=0.0;
end
for n=1:1000 % 内循环开始
w=n*0.0628;
r1(n)=unifrnd(0,6.28); r2(n)=unifrnd(0,6.28);
r3(n)=unifrnd(0,6.28); r4(n)=unifrnd(0,6.28); %随即相角
%相干函数
d=; a=0.736; af=0.147;
sit=3300*(1+(w/(1.5*pi)^2))^(-0.5);
for kk=1:7
a12=cos(w*d(kk)/2500); a13=sin(w*d(kk)/2500);
a1(kk)=a12+a13*i ;
end
for mm=1:7
rr3=a*exp(-2*d(mm)*(1-a+af*a)/(af*sit))+(1-a)*exp(-2*d(mm)*(1-a+af*a)/sit);
rr1=rr3*real(a1(mm)); rr2=rr3*imag(a1(mm));
r(mm)=rr1+rr2*i;
end
%power spectral denstiy function
w1=2.5 ; kc=0.6;
s=(w1^4+4*(kc*w1*w)^2)/((w1^2-w^2)^2+4*(kc*w1*w)^2);
% cross power spectral denstiy function
ss=[ r(1) r(2) r(3) r(4)
r(5) r(1) r(2) r(3)
r(6) r(5) r(1) r(2)
r(7) r(6) r(5) r(1)];
% 计算时间历程系数
l=zeros(4,4); aa=zeros(4,4); amp=zeros(4,4);
l(1,1)=r(1);
for i=1:4
l(i,1)=ss(i,1)/l(1,1);
end
for i=2:4
for j=2:i
if i== j
s1=0.0;
for k=1:i-1
s1=s1+l(i,k)*conj(l(i,k));
end
l(i,j)=(ss(i,j)-s1)^0.5;
else
s2=0.0;
for k=1:j-1
s2=s2+l(i,k)*conj(l(j,k));
end
l(i,j)=(ss(i,j)-s2)/l(j,j);
end
end
end
for i=1:4
for j=1:i
%if i==j
% aa(i,j)=2*(s*0.0628)^0.5 ;
%amp(i,j)=angle(l(i,j));
%else
aa1=real(l(i,j)); aa2=imag(l(i,j));
aa(i,j)=2*(s*0.0628*(aa1^2+aa2^2))^0.5 ;
%amp(i,j)=atan2(imag(l(i,j))/real(l(i,j)))
amp(i,j)=angle(l(i,j));
% end
end
end
u4(n)=(aa(4,1)*cos(w*t+r1(n)+amp(4,1))+aa(4,2)*cos(w*t+r2(n)+amp(4,2))+aa(4,3)*cos(w*t+r3(n)+amp(4,3))+aa(4,4)*cos(w*t+r4(n)+amp(4,4)));
sum4=sum4+u4(n);
end
ff4(kkkk)=sum4;
f4(kkkk)=ff4(kkkk)*yy(kkkk);
end
plot (g,f4)
图形如下：
[ 本帖最后由 ChaChing 于 2009-4-18 20:31 编辑 ] 加速度居然达到150 m/s^2，:@Q ，请高手指点一下，是不是拟合错误，在生成地震波的过程中没有用到fft变换，是不是这个错了？ LZ这个并非编程问题! 个人水平专业有限, 移动版块不知对否? sorry sir！
我是病急乱投医，没准这里隐藏一位博学大师，请多多见谅:loveliness: 我也很迷惑，同上求解
页:
[1]