上官楼 发表于 2009-4-18 18:17

多点人工生成地震波程序

我编了一个多支承点拟合地震波的程序,可是计算结果很大,有些不符合实际,程序如下,请高手指点一下,谢谢!
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 编辑 ]

上官楼 发表于 2009-4-18 18:19

加速度居然达到150 m/s^2,:@Q ,请高手指点一下,是不是拟合错误,在生成地震波的过程中没有用到fft变换,是不是这个错了?

ChaChing 发表于 2009-4-18 20:48

LZ这个并非编程问题! 个人水平专业有限, 移动版块不知对否?

上官楼 发表于 2009-4-19 20:05

sorry sir!
我是病急乱投医,没准这里隐藏一位博学大师,请多多见谅:loveliness:

daizhengyue 发表于 2015-3-15 10:53

我也很迷惑,同上求解
页: [1]
查看完整版本: 多点人工生成地震波程序