多点人工生成地震波程序.txt

谐波小波程序变参数的积分程序多点人工生成地震波程序低通滤波程序经验模式分解降噪程序

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)