多点人工生成地震波程序

上官楼

我编了一个多支承点拟合地震波的程序，可是计算结果很大，有些不符合实际，程序如下，请高手指点一下，谢谢！
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;                                                               
    elseif  t>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=[0 100 200 300 -100 -200 -300 ]; 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变换，是不是这个错了？

ChaChing

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

上官楼

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

daizhengyue

我也很迷惑，同上求解