Skip to content

R
R
提交 52d4ee57 authored 作者: Heyaoyao's avatar Heyaoyao
浏览文件
操作

Initial commit

上级
全部展开 显示空白字符变更
内嵌 并排
正在显示 包含 809 行增加 和 0 行删除
H_Kappa_stacking.m 0 → 100644
浏览文件 @ 52d4ee57
差异被折叠。
Utilities/decon_waterlevel_matlab_v3_trueamp.m 0 → 100644
浏览文件 @ 52d4ee57
function [ RFt, rms, nwl ] = decon_waterlevel_matlab_v3_trueamp( numer, denom, tshift, dt, waterlevel, gaussa )
% calculate frequency
NT = length(numer);
Fs = 1/dt;
% fny = Fs/2; % nyquist
df = Fs/NT;
Q = ceil((NT+1)/2);
freq = df*(0:1:Q-1);
w = 2*pi*[freq,-fliplr(freq(2:NT-Q+1))];
w = reshape(w,[],1);
numer = reshape(numer,[],1);
denom = reshape(denom,[],1);
numerf = fft( numer );
denomf = fft( denom );
% denominator
Df = denomf.* conj(denomf);
dmax = max( real( Df ) );
% add water level correction
phi1 = waterlevel*dmax; % water level
nwl = length( find(Df<phi1) ); % number corrected
% Df( real(Df)<phi1 ) = phi1;
Df = Df + phi1;
% numerator
Nf = numerf .* conj(denomf);
% compute RF
RFf = Nf ./ Df;
rff_nofilter_hyy = abs(RFf);
% compute predicted numerator
pnumerf = RFf.*denomf;
% add gauss filter
gfilter = exp(-w.*w/(4*gaussa*gaussa))/dt*sqrt(pi)/gaussa;
RFf = RFf .* gfilter;
rff_filter_hyy = abs(RFf);
% add phase shift
RFf = RFf .* exp(-1i*w*tshift);
RFt = (1+waterlevel)*real(ifft(RFf));
% RFt = real(ifft(RFf));
% %% temp plot by hyy
% % figure
% fontsize = 32;
% tickfontsize = 28;
% t = (0:NT-1)*dt;
% figa = plot(t, RFt, 'LineWidth', 2);
% set(gca, 'Fontsize', tickfontsize)
% xlabel('Time (sec)', 'FontSize', tickfontsize)
% ylabel('Amplitude', 'FontSize', tickfontsize)
% title('RFs calculated by different ways', 'FontSize', fontsize)
% grid on
% box on
% compute the fit
pnumer=real(ifft(pnumerf));
rms = sum( (pnumer - numer).^2 )/sum(numer.^2);
% %% temp plot by hyy
% nf_hyy = abs(Nf);
% df_hyy = abs(Df);
% t = (0:NT-1)*dt;
% fontsize = 22;
% tickfontsize = 18;
% smooth_scale = 5;
% figure
% xi1 = linspace(min(freq), max(freq), numel(freq)*smooth_scale);
% yi1 = interp1(freq, nf_hyy(1:Q), xi1, 'pchip');
% semilogx(xi1, yi1, 'k', 'LineWidth', 2)
% hold on
% xi2 = linspace(min(freq), max(freq), numel(freq)*smooth_scale);
% yi2 = interp1(freq, df_hyy(1:Q), xi2, 'pchip');
% semilogx(xi2, yi2, 'r', 'LineWidth', 1)
% legend('Nf','Df')
% set(gca, 'Fontsize', tickfontsize)
% title(strcat('Nf and Df (waterlevel = ',num2str(waterlevel),')'), 'FontSize', fontsize);
% ylabel('Amplitude', 'FontSize', tickfontsize);
% xlabel('frequency (Hz)', 'FontSize', tickfontsize)
% xlim([freq(2) freq(end)]);
% grid on
% box on
%
% figure
% subplot(3,1,1)
% xi3 = linspace(min(freq), max(freq), numel(freq)*smooth_scale);
% yi3 = interp1(freq, rff_nofilter_hyy(1:Q), xi3, 'pchip');
% semilogx(xi3, yi3, 'k', 'LineWidth', 2);
% set(gca, 'Fontsize', tickfontsize)
% title(strcat('RFf (waterlevel = ', num2str(waterlevel), ')'), 'FontSize', fontsize);
% ylabel('Amplitude', 'FontSize', tickfontsize);
% xlabel('Frequency (Hz)', 'FontSize', tickfontsize)
% xlim([freq(2) freq(end)]);
% grid on
% box on
%
% subplot(3,1,2)
% xi4 = linspace(min(freq), max(freq), numel(freq)*smooth_scale);
% yi4 = interp1(freq, rff_filter_hyy(1:Q), xi4, 'pchip');
% semilogx(xi4, yi4, 'k', 'LineWidth', 2);
% set(gca, 'Fontsize', tickfontsize)
% title(strcat('RFf (With guass filter(factor = ', num2str(gaussa), '), waterlevel = ', num2str(waterlevel), ')'), 'FontSize', fontsize);
% ylabel('Amplitude', 'FontSize', tickfontsize);
% xlabel('Frequency (Hz)', 'FontSize', tickfontsize)
% xlim([freq(2) freq(end)]);
% grid on
% box on
%
% subplot(3,1,3)
% % single_wig7(t, RFt, 'k', 'LineWidth', 2);
% single_wig7(t, RFt, 0,'blue',[0.82745 0.82745 0.82745],'yes', 'yes')
% set(gca, 'Fontsize', tickfontsize)
% title(strcat('RFt (waterlevel = ', num2str(waterlevel), ')'), 'FontSize', fontsize);
% ylabel('Amplitude', 'FontSize', tickfontsize);
% xlabel('Time (sec)', 'FontSize', tickfontsize)
% axis tight
% grid on
% box on
% qqq =1;
end
Utilities/filtering.m 0 → 100644
浏览文件 @ 52d4ee57
% function [ dat ] = filtering (data, delta, lowf, highf)
% %FILTERING Summary of this function goes here
% % Detailed explanation goes here
%
% % this program is aimed to filter data using a 2-pass butterworth
% % filter. You can also use other filters by changing butter to others
%
%
% fs = 1/delta;
% nyq = fs/2;
% [B,A] = butter(order,[lowf/nyq highf/nyq]);
% dat = filtfilt(B,A,data);
% usage:
% dat = filtering(data, delta, lowf, highf, order) - bandpass filter
% dat = filtering(data, delta, lowf, highf) - bandpass filter
% dat = filtering(data, delta, lowf, 'low', order) - lowpass filter
% dat = filtering(data, delta, lowf, 'low') - lowpass filter
% change 'low' to 'high' for highpass filter
% else no filtering
function [ dat ] = filtering (varargin)
% default butterworth order
order = 2;
if nargin <=3
% fprintf('Number of input parameters less than 4!\n');
return;
else
data = varargin{1};
dat = data;
delta = varargin{2};
f1 = varargin{3};
if ~isnumeric(delta) || ~isnumeric(f1)
% fprintf('Delta or f1 not numeric!\n No filtering\n');
return;
end
fs = 1/delta;
nyq = fs/2;
if nargin == 4
if strcmpi(varargin{4},'high') || strcmpi(varargin{4},'hp')
[B,A] = butter(order,f1/nyq,'high');
elseif strcmpi(varargin{4},'low') || strcmpi(varargin{4},'lp')
[B,A] = butter(order,f1/nyq,'low');
elseif isnumeric(varargin{4})
f2 = varargin{4};
[B,A] = butter(order,[f1/nyq f2/nyq]);
else
% fprintf('Check input para!\n');
return;
end
elseif nargin >=5
if (strcmpi(varargin{4},'high') || strcmpi(varargin{4},'hp')) && isnumeric(varargin{5})
order = varargin{5};
[B,A] = butter(order,f1/nyq,'high');
elseif (strcmpi(varargin{4},'low') || strcmpi(varargin{4},'lp')) && isnumeric(varargin{5})
order = varargin{5};
[B,A] = butter(order,f1/nyq,'low');
elseif isnumeric(varargin{4}) && isnumeric(varargin{5})
f2 = varargin{4};
order = varargin{5};
[B,A] = butter(order,[f1/nyq f2/nyq]);
else
% fprintf('Check input para!\n');
return;
end
end
data = detrend(data);
% taper = tukeywin(length(data),0.1);
% data = data.*reshape(taper,size(data,1),size(data,2));
dat = filtfilt(B,A,double(data));
end
\ No newline at end of file
Utilities/imksac.m 0 → 100644
浏览文件 @ 52d4ee57
差异被折叠。
Utilities/irdsac.m 0 → 100644
浏览文件 @ 52d4ee57
function [ sachd, sacdata ] = irdsac( sacfile )
% irdsac generates a structure 'sachd' and a vector 'sacdata' from an
% exist sac file.
%
% 'sachd' contains the following elements
% -------------------------------------------------------------------------
%delta stla evla data(10) iftype dist xminimum trcLen
%b stlo evlo label(3) idep az xmaximum scale
%e stel evel iztype baz yminimum
%o stdp evdp iinst gcarc ymaximum
%a cmpaz nzyear istreg norid
%t0 cmpinc nzjday ievreg nevid
%t1 kstnm nzhour ievtyp nwfid
%t2 kcmpnm nzmin iqual nxsize
%t3 knetwk nzsec isynth nysize
%t4 nzmsec
%t5 kevnm
%t6 mag
%t7 imagtyp
%t8 imagsrc
%t9
%f
%k0
%ka
%kt1
%kt2
%kt3
%kt4
%kt5
%kt6
%kt7
%kt8
%kt9
%kf
%response is a 10-element array, and trcLen is a scalar.
%
% author: taokai@pku.edu.cn (any suggestion is appreciated)
%%
% Default byte-order
% endian = 'big' big-endian byte order (e.g., UNIX)
% = 'lil' little-endian byte order (e.g., LINUX)
sachd = []; sacdata = [];
fid_lil = fopen(sacfile,'r','ieee-le');
fid_big = fopen(sacfile,'r','ieee-be');
if fid_lil < 0 && fid_big < 0
return;
end
% check header version == 6 and the byte order
%--------------------------------------------------------------------------
% If the header version is not NVHDR == 6 then the sacfile is likely of the
% opposite byte order. This will give h(77) some ridiculously large
% number. NVHDR can also be 4 or 5. In this case it is an old SAC file
% and rsac cannot read this file in. To correct, read the SAC file into
% the newest verson of SAC and w over.
if((1+fseek(fid_lil,304,'bof'))&&(1+fseek(fid_big,304,'bof')))
nvhdr_lil = fread(fid_lil,1,'int32');
nvhdr_big = fread(fid_big,1,'int32');
else
return;
% message = 'Failed to read sac file!';
% error(message)
end
if nvhdr_lil == 6
fid = fid_lil;
% disp('little-endian byte order!')
status = fseek(fid,0,'bof');
fclose(fid_big);
elseif nvhdr_big == 6
fid = fid_big;
disp('big-endian byte order!')
status = fseek(fid,0,'bof');
fclose(fid_lil);
else
return;
% message = ['nvhdr ~= 6: nvhdr_lil = ',num2str(nvhdr_lil),', nvhdr_big = ',num2str(nvhdr_big)];
% error(message)
end
%% read in sac header
h = zeros(302,1);
% read in 70*(single precision real) header variables:
h(1:70) = fread(fid,70,'single');
% read in 35*(single precision interger) header variables:
h(71:105) = fread(fid,35,'int32');
% read in 4*logical(single precision interger) header variables
h(106:110) = fread(fid,5,'int32');
% read in 192*char header variables
h(111:302) = fread(fid,192,'char');
% add header signature for testing files for SAC format
%---------------------------------------------------------------------------
% h(303) = 84;
% h(304) = 65;
% h(305) = 79;
%% write the structured sachd
% read real header variables
%---------------------------------------------------------------------------
sachd.delta = h(1);
sachd.depmin = h(2);
sachd.depmax = h(3);
sachd.scale = h(4);
sachd.odelta = h(5);
sachd.b = h(6);
sachd.e = h(7);
sachd.o = h(8);
sachd.a = h(9);
sachd.t0 = h(11);
sachd.t1 = h(12);
sachd.t2 = h(13);
sachd.t3 = h(14);
sachd.t4 = h(15);
sachd.t5 = h(16);
sachd.t6 = h(17);
sachd.t7 = h(18);
sachd.t8 = h(19);
sachd.t9 = h(20);
sachd.f = h(21);
sachd.resp0 = h(22);
sachd.resp1 = h(23);
sachd.resp2 = h(24);
sachd.resp3 = h(25);
sachd.resp4 = h(26);
sachd.resp5 = h(27);
sachd.resp6 = h(28);
sachd.resp7 = h(29);
sachd.resp8 = h(30);
sachd.resp9 = h(31);
sachd.stla = h(32);
sachd.stlo = h(33);
sachd.stel = h(34);
sachd.stdp = h(35);
sachd.evla = h(36);
sachd.evlo = h(37);
sachd.evel = h(38);
sachd.evdp = h(39);
sachd.mag = h(40);
sachd.user0 = h(41);
sachd.user1 = h(42);
sachd.user2 = h(43);
sachd.user3 = h(44);
sachd.user4 = h(45);
sachd.user5 = h(46);
sachd.user6 = h(47);
sachd.user7 = h(48);
sachd.user8 = h(49);
sachd.user9 = h(50);
sachd.dist = h(51);
sachd.az = h(52);
sachd.baz = h(53);
sachd.gcarc = h(54);
sachd.depmen = h(57);
sachd.cmpaz = h(58);
sachd.cmpinc = h(59);
sachd.xminimum = h(60);
sachd.xmaximum = h(61);
sachd.yminimum = h(62);
sachd.ymaximum = h(63);
% read integer header variables
%---------------------------------------------------------------------------
sachd.nzyear = round(h(71));
sachd.nzjday = round(h(72));
sachd.nzhour = round(h(73));
sachd.nzmin = round(h(74));
sachd.nzsec = round(h(75));
sachd.nzmsec = round(h(76));
sachd.nvhdr = round(h(77));
sachd.norid = round(h(78));
sachd.nevid = round(h(79));
sachd.npts = round(h(80));
sachd.nwfid = round(h(82));
sachd.nxsize = round(h(83));
sachd.nysize = round(h(84));
sachd.iftype = round(h(86));
sachd.idep = round(h(87));
sachd.iztype = round(h(88));
sachd.iinst = round(h(90));
sachd.istreg = round(h(91));
sachd.ievreg = round(h(92));
sachd.ievtyp = round(h(93));
sachd.iqual = round(h(94));
sachd.isynth = round(h(95));
sachd.imagtyp = round(h(96));
sachd.imagsrc = round(h(97));
%read logical header variables
%---------------------------------------------------------------------------
sachd.leven = round(h(106));
sachd.lpspol = round(h(107));
sachd.lovrok = round(h(108));
sachd.lcalda = round(h(109));
%read character header variables
%---------------------------------------------------------------------------
sachd.kstnm = char(h(111:118));
sachd.kevnm = char(h(119:134));
sachd.khole = char(h(135:142));
sachd.ko = char(h(143:150));
sachd.ka = char(h(151:158));
sachd.kt0 = char(h(159:166));
sachd.kt1 = char(h(167:174));
sachd.kt2 = char(h(175:182));
sachd.kt3 = char(h(183:190));
sachd.kt4 = char(h(191:198));
sachd.kt5 = char(h(199:206));
sachd.kt6 = char(h(207:214));
sachd.kt7 = char(h(215:222));
sachd.kt8 = char(h(223:230));
sachd.kt9 = char(h(231:238));
sachd.kf = char(h(239:246));
sachd.kuser0 = char(h(247:254));
sachd.kuser1 = char(h(255:262));
sachd.kuser2 = char(h(263:270));
sachd.kcmpnm = char(h(271:278));
sachd.knetwk = char(h(279:286));
sachd.kdatrd = char(h(287:294));
sachd.kinst = char(h(295:302));
%% read in amplitudes
if nargout == 2
sacdata = fread(fid,sachd.npts,'single');
end
fclose(fid);
end
\ No newline at end of file
Utilities/single_wig7.m 0 → 100644
浏览文件 @ 52d4ee57
function h = single_wig7(time,data,offset,upcolor,dncolor,isupfill,isdnfill )
%SINGLE_WIG2 Summary of this function goes here
% Detailed explanation goes here
dark = [0 0 0];
gray = [0.7 0.7 0.7];
if nargin <= 3
upcolor = dark;
dncolor = gray;
isupfill = 'yes';
isdnfill = 'yes';
elseif nargin == 4
dncolor = gray;
isupfill = 'yes';
isdnfill = 'yes';
end
data = reshape(data,[],1);
time = reshape(time,[],1);
[t0,d0] = curveintersect(time, data, time, zeros(size(time)));
datatmp = [zeros(size(d0));data];
[timetmp,indx] = sort([t0;time],'ascend');
datatmp = datatmp(indx);
timetmp = [timetmp(1);timetmp;timetmp(end)];
datatmp = [0;datatmp;0];
dataup = datatmp;
dataup(find(dataup<0)) = 0;
dataup(1) = -max(abs(dataup))*1e-4;
dataup(end) = -max(abs(dataup))*1e-4;
datadn = datatmp;
datadn(find(datadn>0)) = 0;
datadn(1) = max(abs(datadn))*1e-4;
datadn(end) = max(abs(datadn))*1e-4;
if strcmpi(isupfill,'yes') || strcmpi(isupfill,'y')
h(1) = patch(timetmp,dataup + offset,upcolor,'Facecolor',upcolor,'Edgecolor','none','LineWidth',0.1);
else
h(1) = plot(timetmp,dataup+offset,'color',upcolor,'Linewidth',0.1);
end
if strcmpi(isdnfill,'yes') || strcmpi(isdnfill,'y')
h(2) = patch(timetmp,datadn + offset,dncolor,'Facecolor',dncolor,'Edgecolor','none','LineWidth',0.1);
else
h(2) = plot(timetmp,datadn+offset,'color',dncolor,'Linewidth',0.1);
end
end
Utilities/time2depth.m 0 → 100644
浏览文件 @ 52d4ee57
%% This shell convert the waveform from time domain to depth domain
%% Written by Heyaoyao on 2025/7/4
function [H, amp] = time2depth(time, data, tshift, vp, k, rayp, minH, maxH)
vs = vp/k;
delta = time(2)-time(1);
time_cut = time(round(tshift/delta)+1:end);
data_cut = data(round(tshift/delta)+1:end);
nt = length(time_cut);
H_temp = zeros(1, nt);
amp_temp = zeros(1, nt);
for ii = 1:nt
diff_t = time_cut(ii);
H_temp(ii) = diff_t/(sqrt(vs^(-2)-rayp^2)-sqrt(vp^(-2)-rayp^2));
amp_temp(ii) = data_cut(ii);
end
H = linspace(minH, maxH, 2501);
amp = interp1(H_temp, amp_temp, H, 'linear');
% def time_to_depth_Ps(t0,time,data,rayp,arfa,beta,hmin=0,hmax=600):
% dt = time[1] - time[0]
% nt = len(time)
% it0 = int((t0-time[0]) /dt + 1 )
% index = np.arange(it0,nt)
% H = []
% amp = []
% for i,tnumber in enumerate(index):
% thist = (tnumber - it0) *dt
% H.append(thist/(np.sqrt(beta**-2 - rayp**2)-np.sqrt(arfa**-2 - rayp**2)))
% amp.append(data[tnumber])
% H,amp = np.array(H),np.array(amp)
% newH = np.linspace(0,600,6000)
% f = interp1d(H,amp,kind='linear',bounds_error=False,fill_value=0)
% newamp = f(newH)
% return newH,newamp
end
\ No newline at end of file
Utilities/y_Hkappa_modified.m 0 → 100644
浏览文件 @ 52d4ee57
function [ Hk ] = y_Hkappa_modified( RF, time, rayp, H, k, vp, weight )
% H-kappa search for P-receiver functions
% Chunquan Yu, @MIT-EAPS, 06/13
% RF: nt * ntrace matrix
% time: nt * 1 vector
% rayp: ntrace * 1 vector
% % H k range
% H = 20:0.1:80; % crustal thickness for searching
% k = 1.5:0.01:2.0; % kappa = Vp/Vs for searching
% vp = 6.3; % assumed crustal P-wave speed
% % weight for 3 phases: Ps, PpPs, PpSs+PsPs
% weight = [2/5 2/5 1/5];
% normalize data
for i = 1:size(RF,2)
RF(:,i) = RF(:,i)/(max(RF(:,i))-min(RF(:,i)));
end
delta = time(2)-time(1);
r1 = zeros(length(k),length(H));
r2 = zeros(length(k),length(H));
r3 = zeros(length(k),length(H));
for m = 1:length(H)
for n = 1:length(k)
% if H(m) < 40 || k(n) < 1.6 || k(n) > 1.9 || H(m) > 49
% continue
% end
vs = vp/k(n);
for i = 1:size(RF,2) % loop for all receiver functions
partP = sqrt(vp^(-2) - rayp(i)^2);
partS = sqrt(vs^(-2) - rayp(i)^2);
tPs = H(m)*(partS-partP);
tPpPs = H(m)*(partP+partS);
tPpSs = 2*H(m)*partS;
% if tPs < 5 || tPpPs < 18 || tPpPs > 18.5
% continue
% end
RFi = RF(:,i);
r1(n,m) = r1(n,m) + RFi(round((tPs-time(1))/delta) + 1);
r2(n,m) = r2(n,m) + RFi(round((tPpPs-time(1))/delta) + 1);
r3(n,m) = r3(n,m) + RFi(round((tPpSs-time(1))/delta) + 1);
end
end
end
% % Hk, r1, r2 should be positive, and r3 should be negative, otherwise
% % setting them to be zero
% r1(find(r1<0))=0;
% r2(find(r2<0))=0;
% r3(find(r3>0))=0;
Hk = weight(1)*r1 + weight(2)*r2 - weight(3)*r3;
% Hk(find(Hk<0)) = 0;
end
cal_RFs.m 0 → 100644
浏览文件 @ 52d4ee57
%% calculating the receiver function in time domain by Yaoyao He, 2024/7/6
clear; clc; close all;
% %% convert event set to station set
% evset_dir = "D:\SUSTech businesse\observational_seismology\XE_work\sac_eventset";
% staset_dir = "D:\SUSTech businesse\observational_seismology\XE_work\sac_stationset";
% if ~exist(staset_dir, 'dir')
% mkdir(staset_dir);
% end
% evlist = string(ls(fullfile(evset_dir,'2*')));
% for i = 1:length(evlist)
% % list_r_Decon_PF = textscan(fopen(fullfile(evset_dir, evlist(i), 'list_r_Decon_PF.txt'), 'r'),...
% % '%s %f %f %f %f %s %s %f %f %f %f %f %f %f %f %f %f\n', 'HeaderLines', 1);
% % decon_filename = char(list_r_Decon_PF{1});
% saclist = char(ls(fullfile(evset_dir,evlist(i),'*.SAC')));
% for j = 1:size(saclist, 1)
% parts = split(saclist(j,:), '.');
% sta_no = parts{2};
% if ~exist(fullfile(staset_dir,sta_no),'dir')
% mkdir(fullfile(staset_dir,sta_no));
% end
% copyfile(fullfile(evset_dir,evlist(i),saclist(j, :)), fullfile(staset_dir,sta_no));
% end
% end
%% show the raw data
% I/O dir
sacdir= "D:\YB\yunnan2\F0_6.0_waterlevel_0.01\3_selected_stationset_v2";
stlist = string(ls(fullfile(sacdir,'NK*')));
outputdir = fullfile(fileparts(sacdir), '3_selected_stationset_v2_multiband_RFs');
if ~exist(outputdir, 'dir')
mkdir(outputdir);
end
% basic parameters
new_delta = 0.1;
timewin = [-30 60];
time_cut = reshape(timewin(1):new_delta:timewin(2),[],1);
fl = 0.02; fh = 1; order = 2;
np = 1000;
mod = 'ak135'; thmax = 15;
em = refinemodel(set_vmodel_v2(mod), thmax);
R = 6371;
% decon parameters
deconwin = [-30 60];
time_decon = reshape(deconwin(1):new_delta:deconwin(2),[],1);
F0 = 6; TDEL = 30; itmax = 100; minerr = 10^(-2); waterlevel = 0.001 ;
sourcedir = "D:\YB\yunnan2\F0_6.0_waterlevel_0.01\2_z_seis_cutdata_180s_1st_selection";
for i = 1:length(stlist)
fid = fopen(fullfile(sacdir, stlist(i), 'list_r_Decon_PF.txt'), 'r');
if fid == -1
continue
end
list_r_Decon_PF = textscan(fid, '%s %f %f %f %f %s %s %f %f %f %f %f %f %f %f %f %f %f %f \n', ...
'Commentstyle', 'shell', 'HeaderLines', 1);
fclose(fid);
if isempty(list_r_Decon_PF{1})
continue
end
if ~exist(fullfile(outputdir, stlist(i)), 'dir')
mkdir(fullfile(outputdir, stlist(i)));
end
sacfiles = list_r_Decon_PF{1};
num_file = length(sacfiles);
for j = 1:num_file
parts = string(split(sacfiles(j), '.'));
file_R = strcat(parts(1), '.', parts(2), '.', parts(3), '.', parts(4), '.', 'R', '.', parts(6));
file_Z = strcat(parts(1), '.', parts(2), '.', parts(3), '.', parts(4), '.', 'Z', '.', parts(6));
[sachdr, data_raw_r] = irdsac(fullfile(sourcedir, strcat(parts(3), '.', parts(4)), file_R));
[sachdz, data_raw_z] = irdsac(fullfile(sourcedir, strcat(parts(3), '.', parts(4)), file_Z));
stla = sachdr.stla; stlo = sachdr.stlo;
evla = sachdr.evla; evlo = sachdr.evlo; evdp = sachdr.evdp;
% modified by Chunquan Yu (2019/11/09) to take into account
[dist,baz] = distance(stla,stlo, evla,evlo, [6378.137*180/pi/6371 0.0818191908426215]);
[rayp, taup, Xp]= y_get_p_P_Pdiff (evdp, np, em);
tt = taup + rayp.*Xp;
dd = Xp*180/pi;
theo_tt = interp1db(dist, dd,tt);
time_raw = [0:sachdr.npts-1]' * sachdr.delta + sachdr.b - sachdr.o - theo_tt;
data_cut_r = interp1(time_raw, data_raw_r, time_cut, 'linear', nan);
data_cut_z = interp1(time_raw, data_raw_z, time_cut, 'linear', nan);
% filtered data
data_cut_r_filter = filtering(data_cut_r, new_delta, fl, fh, order);
data_cut_z_filter = filtering(data_cut_z, new_delta, fl, fh, order);
%% Decon part
data_decon_r = interp1(time_cut, data_cut_r, time_decon, 'linear', 0);
data_decon_z = interp1(time_cut, data_cut_z, time_decon, 'linear', 0);
data_decon_r = detrend(data_decon_r) .* tukeywin(length(time_decon),new_delta);
data_decon_z = detrend(data_decon_z) .* tukeywin(length(time_decon),new_delta);
% Frequency domian
[eqr, rms_r, nwl_r] = decon_waterlevel_matlab_v3_trueamp( data_decon_r, data_decon_z, TDEL, new_delta, waterlevel, F0);
[aftn, rms_z, nwl_z] = decon_waterlevel_matlab_v3_trueamp( data_decon_z, data_decon_z, TDEL, new_delta, waterlevel, F0);
newsachdr = sachdr;
newsachdr.a = sachdr.b + theo_tt;
newsachdr.b = newsachdr.a + timewin(1);
newsachdr.e = newsachdr.a + timewin(2);
newsachdr.npts = length(time_decon);
newsachdr.delta = new_delta;
newsachdz = sachdz;
newsachdz.a = sachdz.b + theo_tt;
newsachdz.b = newsachdz.a + timewin(1);
newsachdz.e = newsachdz.a + timewin(2);
newsachdz.npts = length(time_decon);
newsachdz.delta = new_delta;
deconR_filename = strcat('gaussa_', num2str(F0), '_', file_R, '.PF.Decon');
deconZ_filename = strcat('gaussa_', num2str(F0), '_', file_Z, '.PF.Decon');
imksac(newsachdr, eqr, fullfile(outputdir, stlist(i), deconR_filename));
% imksac(newsachdz, aftn, fullfile(outputdir, stalist(i), deconZ_filename));
end
end
data_demo/NK.S0001/list_r_Decon 0 → 100644
浏览文件 @ 52d4ee57
NK.S0001.20210530Z133329.37T.R.sac.PF.Decon
NK.S0001.20210603Z100958.29T.R.sac.PF.Decon
NK.S0001.20210604Z134513.27T.R.sac.PF.Decon
NK.S0001.20210605Z225738.23T.R.sac.PF.Decon
NK.S0001.20210608Z050032.30T.R.sac.PF.Decon
NK.S0001.20210608Z105637.58T.R.sac.PF.Decon
NK.S0001.20210608Z152901.05T.R.sac.PF.Decon
NK.S0001.20210611Z022343.04T.R.sac.PF.Decon
NK.S0001.20210613Z092759.72T.R.sac.PF.Decon
NK.S0001.20210614Z024459.14T.R.sac.PF.Decon
NK.S0001.20210616Z044307.45T.R.sac.PF.Decon
NK.S0001.20210620Z110824.06T.R.sac.PF.Decon
NK.S0001.20210621Z221416.82T.R.sac.PF.Decon
NK.S0001.20210624Z022456.11T.R.sac.PF.Decon
NK.S0001.20210624Z144956.44T.R.sac.PF.Decon
NK.S0001.20210703Z170427.76T.R.sac.PF.Decon
NK.S0001.20210704Z161847.79T.R.sac.PF.Decon
NK.S0001.20210709Z133110.14T.R.sac.PF.Decon
NK.S0001.20210710Z004358.95T.R.sac.PF.Decon
NK.S0001.20210710Z005110.81T.R.sac.PF.Decon
NK.S0001.20210718Z000935.15T.R.sac.PF.Decon
NK.S0001.20210718Z143419.28T.R.sac.PF.Decon
NK.S0001.20210719Z152231.00T.R.sac.PF.Decon
NK.S0001.20210719Z203049.87T.R.sac.PF.Decon
NK.S0001.20210719Z230210.92T.R.sac.PF.Decon
NK.S0001.20210721Z142559.55T.R.sac.PF.Decon
NK.S0001.20210721Z174402.86T.R.sac.PF.Decon
NK.S0001.20210726Z035205.18T.R.sac.PF.Decon
NK.S0001.20210726Z120910.01T.R.sac.PF.Decon
NK.S0001.20210729Z061549.19T.R.sac.PF.Decon
NK.S0001.20210729Z203244.85T.R.sac.PF.Decon
NK.S0001.20210801Z043127.10T.R.sac.PF.Decon
NK.S0001.20210802Z050121.98T.R.sac.PF.Decon
NK.S0001.20210802Z150425.06T.R.sac.PF.Decon
NK.S0001.20210803Z203340.53T.R.sac.PF.Decon
NK.S0001.20210803Z204316.36T.R.sac.PF.Decon
NK.S0001.20210803Z221057.60T.R.sac.PF.Decon
NK.S0001.20210804Z025623.23T.R.sac.PF.Decon
NK.S0001.20210804Z044022.78T.R.sac.PF.Decon
NK.S0001.20210804Z055507.92T.R.sac.PF.Decon
NK.S0001.20210805Z050607.47T.R.sac.PF.Decon
NK.S0001.20210806Z110818.99T.R.sac.PF.Decon
NK.S0001.20210807Z145130.90T.R.sac.PF.Decon
NK.S0001.20210810Z212151.79T.R.sac.PF.Decon
NK.S0001.20210811Z150257.46T.R.sac.PF.Decon
NK.S0001.20210811Z174612.73T.R.sac.PF.Decon
NK.S0001.20210812Z155334.33T.R.sac.PF.Decon
NK.S0001.20210814Z012745.16T.R.sac.PF.Decon
NK.S0001.20210814Z115743.45T.R.sac.PF.Decon
NK.S0001.20210818Z101005.24T.R.sac.PF.Decon
data_demo/NK.S0001/list_r_Decon_PF.txt 0 → 100644
浏览文件 @ 52d4ee57
#filename theo_tt tshift obs_tt polarity stnm netwk rayp stla stlo stel evla evlo evdp dist az baz snr0 xcoeff0
NK.S0001.20210608Z105637.58T.R.sac.PF.Decon 582.757055 0.000000 582.757055 1 S0001 NK 0.062247 25.587839 99.890030 2.564900 -5.245700 151.316498 134.000000 58.605788 304.168399 114.080331 1.000000 1.000000
NK.S0001.20210719Z152231.00T.R.sac.PF.Decon 564.792477 0.000000 564.792477 1 S0001 NK 0.065775 25.587839 99.890030 2.564900 -3.360000 146.869995 1.000000 53.888271 305.158040 115.263242 1.000000 1.000000
NK.S0001.20210721Z142559.55T.R.sac.PF.Decon 563.195382 0.000000 563.195382 1 S0001 NK 0.065932 25.587839 99.890030 2.564900 -3.340000 146.591202 0.000000 53.649172 305.279949 115.442410 1.000000 1.000000
NK.S0001.20210812Z155334.33T.R.sac.PF.Decon 385.315972 0.000000 385.315972 1 S0001 NK 0.078748 25.587839 99.890030 2.564900 6.206400 127.156502 49.799999 32.441300 309.503417 121.790786 1.000000 1.000000
NK.S0001.20210811Z174612.73T.R.sac.PF.Decon 380.182576 0.000000 380.182576 1 S0001 NK 0.078918 25.587839 99.890030 2.564900 6.442000 126.612297 51.299999 31.874705 309.705282 122.100244 1.000000 1.000000
NK.S0001.20210802Z050121.98T.R.sac.PF.Decon 492.777671 0.000000 492.777671 1 S0001 NK 0.071785 25.587839 99.890030 2.564900 -4.503200 134.009903 10.000000 44.657918 313.857674 127.200690 1.000000 1.000000
NK.S0001.20210604Z134513.27T.R.sac.PF.Decon 419.582391 0.000000 419.582391 1 S0001 NK 0.077128 25.587839 99.890030 2.564900 0.340000 126.289703 8.600000 35.800042 316.624452 130.450297 1.000000 1.000000
NK.S0001.20210616Z044307.45T.R.sac.PF.Decon 462.245898 0.000000 462.245898 1 S0001 NK 0.074181 25.587839 99.890030 2.564900 -3.558000 129.523804 7.000000 40.840774 316.908082 130.927709 1.000000 1.000000
NK.S0001.20210611Z022343.04T.R.sac.PF.Decon 403.715027 0.000000 403.715027 1 S0001 NK 0.077579 25.587839 99.890030 2.564900 0.051900 124.297600 60.099998 34.672379 318.980568 133.345470 1.000000 1.000000
NK.S0001.20210608Z050032.30T.R.sac.PF.Decon 384.525026 0.000000 384.525026 1 S0001 NK 0.077403 25.587839 99.890030 2.564900 0.452600 123.647797 196.800003 33.945468 319.316662 133.756152 1.000000 1.000000
NK.S0001.20210709Z133110.14T.R.sac.PF.Decon 392.780661 0.000000 392.780661 1 S0001 NK 0.077563 25.587839 99.890030 2.564900 0.131200 123.636597 129.800003 34.181129 319.636171 134.142217 1.000000 1.000000
NK.S0001.20210726Z035205.18T.R.sac.PF.Decon 399.555214 0.000000 399.555214 1 S0001 NK 0.078137 25.587839 99.890030 2.564900 -0.713900 121.981903 26.900000 33.791129 322.337394 137.392146 1.000000 1.000000
NK.S0001.20210726Z120910.01T.R.sac.PF.Decon 399.280243 0.000000 399.280243 1 S0001 NK 0.078119 25.587839 99.890030 2.564900 -0.732600 121.993202 30.400000 33.812761 322.340845 137.396475 1.000000 1.000000
NK.S0001.20210703Z170427.76T.R.sac.PF.Decon 470.770438 0.000000 470.770438 1 S0001 NK 0.073504 25.587839 99.890030 2.564900 14.101900 56.839298 10.000000 41.932992 67.327380 262.660311 1.000000 1.000000
NK.S0001.20210614Z024459.14T.R.sac.PF.Decon 533.059179 0.000000 533.059179 1 S0001 NK 0.068415 25.587839 99.890030 2.564900 12.821300 48.600899 10.000000 49.826451 67.272047 265.337197 1.000000 1.000000
NK.S0001.20210801Z043127.10T.R.sac.PF.Decon 620.652688 0.000000 620.652688 1 S0001 NK 0.060376 25.587839 99.890030 2.564900 36.360001 27.049999 10.000000 62.057555 77.508092 299.282072 1.000000 1.000000
data_demo/NK.S0001/outpara.txt 0 → 100644
浏览文件 @ 52d4ee57
#Parameters used for processing
#filter
filter_type = Two-Pass
fl = 0.020000
fh = 1.000000
order = 2
data_demo/NK.S0002/list_r_Decon 0 → 100644
浏览文件 @ 52d4ee57
NK.S0002.20210528Z232104.11T.R.sac.PF.Decon
NK.S0002.20210529Z010241.08T.R.sac.PF.Decon
NK.S0002.20210530Z133329.37T.R.sac.PF.Decon
NK.S0002.20210603Z100958.29T.R.sac.PF.Decon
NK.S0002.20210604Z134513.27T.R.sac.PF.Decon
NK.S0002.20210605Z225738.23T.R.sac.PF.Decon
NK.S0002.20210608Z050032.30T.R.sac.PF.Decon
NK.S0002.20210608Z105637.58T.R.sac.PF.Decon
NK.S0002.20210608Z152901.05T.R.sac.PF.Decon
NK.S0002.20210611Z022343.04T.R.sac.PF.Decon
NK.S0002.20210613Z092759.72T.R.sac.PF.Decon
NK.S0002.20210614Z024459.14T.R.sac.PF.Decon
NK.S0002.20210616Z044307.45T.R.sac.PF.Decon
NK.S0002.20210619Z115905.54T.R.sac.PF.Decon
NK.S0002.20210620Z110824.06T.R.sac.PF.Decon
NK.S0002.20210621Z014411.86T.R.sac.PF.Decon
NK.S0002.20210624Z022456.11T.R.sac.PF.Decon
NK.S0002.20210624Z144956.44T.R.sac.PF.Decon
NK.S0002.20210625Z182838.30T.R.sac.PF.Decon
NK.S0002.20210628Z015640.29T.R.sac.PF.Decon
NK.S0002.20210701Z150634.87T.R.sac.PF.Decon
data_demo/NK.S0002/list_r_Decon_PF.txt 0 → 100644
浏览文件 @ 52d4ee57
#filename theo_tt tshift obs_tt polarity stnm netwk rayp stla stlo stel evla evlo evdp dist az baz snr0 xcoeff0
NK.S0002.20210604Z134513.27T.R.sac.PF.Decon 419.545912 0.000000 419.545912 1 S0002 NK 0.077130 25.595501 99.903458 3.121000 0.340000 126.289703 8.600000 35.795766 316.647847 130.475058 1.000000 1.000000
NK.S0002.20210616Z044307.45T.R.sac.PF.Decon 462.211621 0.000000 462.211621 1 S0002 NK 0.074184 25.595501 99.903458 3.121000 -3.558000 129.523804 7.000000 40.836611 316.929063 130.949367 1.000000 1.000000
data_demo/NK.S0002/outpara.txt 0 → 100644
浏览文件 @ 52d4ee57
#Parameters used for processing
#filter
filter_type = Two-Pass
fl = 0.020000
fh = 1.000000
order = 2
data_demo/NK.S0003/list_r_Decon 0 → 100644
浏览文件 @ 52d4ee57
差异被折叠。
data_demo/NK.S0003/list_r_Decon_PF.txt 0 → 100644
浏览文件 @ 52d4ee57
差异被折叠。
data_demo/NK.S0003/outpara.txt 0 → 100644
浏览文件 @ 52d4ee57
差异被折叠。
data_demo/NK.S0004/SAC_list_r_Decon_PF_STACK 0 → 100644
浏览文件 @ 52d4ee57
差异被折叠。
data_demo/NK.S0004/list_r_Decon 0 → 100644
浏览文件 @ 52d4ee57
差异被折叠。
data_demo/NK.S0004/list_r_Decon_PF.txt 0 → 100644
浏览文件 @ 52d4ee57
差异被折叠。
data_demo/NK.S0004/list_r_Decon_STACK 0 → 100644
浏览文件 @ 52d4ee57
差异被折叠。
data_demo/NK.S0004/list_r_Decon_STACK_PF.txt 0 → 100644
浏览文件 @ 52d4ee57
差异被折叠。
data_demo/NK.S0004/outpara.txt 0 → 100644
浏览文件 @ 52d4ee57
差异被折叠。
data_demo/NK.S0005/list_r_Decon 0 → 100644
浏览文件 @ 52d4ee57
差异被折叠。
data_demo/NK.S0005/list_r_Decon_PF.txt 0 → 100644
浏览文件 @ 52d4ee57
差异被折叠。
data_demo/NK.S0005/outpara.txt 0 → 100644
浏览文件 @ 52d4ee57
差异被折叠。
data_demo/NK.S0006/list_r_Decon 0 → 100644
浏览文件 @ 52d4ee57
差异被折叠。
data_demo/NK.S0006/list_r_Decon_PF.txt 0 → 100644
浏览文件 @ 52d4ee57
差异被折叠。
data_demo/NK.S0006/outpara.txt 0 → 100644
浏览文件 @ 52d4ee57
差异被折叠。
prof_analysis_of_RFs.m 0 → 100644
浏览文件 @ 52d4ee57
差异被折叠。
Markdown 格式
0% 或
您添加了 0 人 到此讨论。请谨慎行事。
请先完成此评论的编辑！
请 注册 或者 后发表评论