demultiple.m

著名的seismiclab的代码 是地震学研究人员必备的工具

load dataclean.mat;

data = Data(:,10:20);
hw = h(10:20);
h = hw;
[nt,nh] = size(data);

flow = 1.; fhigh = 60;

mu = 0.001;

 dq = 4*dt;
qmin = -4*dq;
qmax = .4;
q = qmin:dq:qmax;
nq = length(q);

mu = 0.01;
[m1,M1f,Q1] = inv_taup(data,dt,h,q,2,flow,fhigh,mu,'ls');

mu = 0.01;
[m2,M2f,Q2] = inv_taup(data,dt,h,q,2,flow,fhigh,mu,'hr');

% Filtering ..

qd = 0.1
icut = floor((qd -qmin)/dq)+1
m1f= m1;  m1f(:,1:icut) = 0;
m2f= m2;  m2f(:,1:icut) = 0;


[d1p] = for_taup(m1f,dt,h,q,2,flow,fhigh);

[d2p] = for_taup(m2f,dt,h,q,2,flow,fhigh);

t = 0:dt:(nt-1)*dt;

figure(1);clf;
figure(1);  subplot(131); imagesc(h,t,data); title('D');ylabel('t (s)')
xlabel('h (m)');
figure(1);  subplot(132); imagesc(h,t,d1p); title('M')
figure(1);  subplot(133); imagesc(h,t,data-d1p); title('P')

figure(2);clf;
figure(2);  subplot(131); imagesc(h,t,data); title('D');ylabel('t (s)')
xlabel('h (m)');
figure(2);  subplot(132); imagesc(h,t,d2p); title('M')
figure(2);  subplot(133); imagesc(h,t,data-d2p); title('P')


colormap(seismic);
figure(3);  clf;
figure(3);  subplot(121); imagesc(q,t,[m1]); title('LS PRT')     ;xlabel('q (s)')
ylabel('\tau (s)')
figure(3);  subplot(122); imagesc(q,t,[m2]); title('NI-HR PRT');xlabel('q (s)')
ylabel('\tau (s)')
colormap(seismic);


qcut = q(icut);
subplot(121); hold on; plot( [qcut,qcut],[0,max(t)],'-');
subplot(122); hold on; plot( [qcut,qcut],[0,max(t)],'-');

figure(4);clf; plot(data(:,20)); hold on; 
               plot(d1p(:,20)+1);  
               plot(d2p(:,20)+2);