multipointmanualprocess.rar_地震 小波_滤波 地震_经验小波_经验小波分解_谐波 小波

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)