Skip to content

Commit

Permalink
Bug fixed
Browse files Browse the repository at this point in the history
  • Loading branch information
Karina Loviknes committed Aug 18, 2020
1 parent f0e5ea6 commit 800f9ef
Show file tree
Hide file tree
Showing 7 changed files with 75 additions and 48 deletions.
14 changes: 7 additions & 7 deletions EGF_tc_toolbox/RUN.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,27 +2,27 @@
settingsfile='settingsfile.txt';

% ESTIMATE THE GREEN'S FUNCTION:
egfs = estimate_GF(settingsfile);
%egfs = estimate_GF(settingsfile);

% MEASURE THE TIME SHIFT
delays = measure_timeshift(settingsfile);

% PLOT DAILY CROSS CORRELATION AND TIME SHIFTS
h = plot_egf_td(settingsfile);
h = plot_egf_td(settingsfile,'all');
h = plot_egf_td(settingsfile,'Daily');
h = plot_egf_td(settingsfile,'Frequency');
%h = plot_egf_td(settingsfile,'Daily');
%h = plot_egf_td(settingsfile,'Frequency');

% INVERT TO FIND THE TIME DELY OF EACH STATION
[fdelay fdelayc] = invert_TD(settingsfile,'plot');
%% INVERT TO FIND THE TIME DELY OF EACH STATION
%[fdelay fdelayc] = invert_TD(settingsfile,'plot');

% Plot the results with distance, filters and ampliude spectrum
h = plot_distance(settingsfile);
h = apply_filterband(settingsfile,'stack');
h = apply_filterband(settingsfile,'daily');

%%
% Give the drift on a specified date:
dato = {'2017-07-01'}
dato = {'2015-10-01'}
[tdd tddc dayn]= date_select('settingsfile.txt',char(dato));

% Corrcet for measured time shift:
Expand Down
2 changes: 1 addition & 1 deletion EGF_tc_toolbox/make_reference.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@
% ref = the reference trace
% numref = the number of reference traces
%
% Written by Karina Løviknes
% Written by Karina Loviknes
%

number_of_days = length(daily(:,1)); % Number of days
Expand Down
23 changes: 13 additions & 10 deletions EGF_tc_toolbox/measure_timeshift.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -144,12 +144,13 @@
% Normalize
r_neg = r_neg1/max(abs(r_neg1));
r_pos = r_pos1/max(abs(r_pos1));

% Find the static timing error (the unsymmetry of the
% reference):
delay_stat(it) = cl(find(abs(reff)==max(abs(reff))));
% Measure over the entire waveform
refwhn = reff/max(abs(reff));

% Find the static timing error (the unsymmetry of the
% reference): NOT READY!
[td_ref,lag_ref] = cross_conv(r_pos,r_neg,Fq);
delay_stat(it) = 0; %lag_ref(find(abs(td_ref)==max(abs(td_ref))))
% delay_stat(it) = cl(find(abs(reff)==max(abs(reff))));

% Loop over the daily cross correlations
for d=1:num_corr
Expand All @@ -172,7 +173,6 @@
% Normalize
s_neg = s_neg1/max(abs(s_neg1));
s_pos = s_pos1/max(abs(s_pos1));
% Measure over the entire waveform
crconvfn = egff(k,:)/max(abs(egff(k,:)));

% Calculate timing difference:
Expand Down Expand Up @@ -218,13 +218,13 @@
delay_pos(k)<=(-(delay_neg(k)-2))
% Calculate the average delay:
delay_dyn0(k) = (delay_pos(k)-delay_neg(k))/2;
delay_dyn(k) = 0;
delay_dyn(k) = (delay_pos(k)-delay_neg(k))/2;

type=[type 's'];
else
% If not the delay is set to zero or the
%delay of the previous day:
delay_dyn0(k) = NaN;
delay_dyn0(k) = NaN;
if k>1
% The delay can not measured, and is
% therefore set as the same as the
Expand All @@ -233,7 +233,7 @@
else
% If it is the first day the delay is
% set to zero
delay_dyn(k) = NaN;
delay_dyn(k) = 0;
end
type = [type '0'];

Expand Down Expand Up @@ -327,12 +327,15 @@
end
timeshift = dd_fit_eval;
else
dd_fit = polyfit(dd1,delay_dyn,1);
dd_fit_eval = polyval(dd_fit,dd1);

timeshift = delay_dyn;
end

for kk = 1:k
% Correct for found timing errors:
t01 = -dd_fit_eval(kk)/Fq;
t01 = -timeshift(kk)/Fq;
omgc = exp([0:Lc-1]*Fq/Lc*-1i*2*pi*t01);
shift1 = fft(egf(kk,:)).*omgc;
shift2 = ifft(shift1,'symmetric');
Expand Down
7 changes: 3 additions & 4 deletions EGF_tc_toolbox/plot_distance.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
% Sub-function: read_settings.m, str2filename.m, read_info_pz.m and
% calc_dist.m
%
% Written by Karina Løviknes
% Written by Karina Loviknes
%

% Read all the default values from the settings file
Expand Down Expand Up @@ -92,7 +92,7 @@
[longB,latB] = read_info_pz(pz_file_B,'Coordinates');

% Calculate distance:
in_st_dist(sp) = calc_dist(latA,longA,latB,longB)*1000;
in_st_dist(sp) = calc_dist(latA,longA,latB,longB);
end
% Make sure the stations are cross correlted with the rigth number of stations:
if ii >= num_stat_cc
Expand All @@ -103,10 +103,9 @@
end
end

maxdist = max(in_st_dist)
maxdist = max(in_st_dist);
btwdist = maxdist/sp;
mindist = min(in_st_dist);
in_st_dist

h = figure;
for ii=1:sp
Expand Down
56 changes: 40 additions & 16 deletions EGF_tc_toolbox/plot_egf_td.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,12 +13,12 @@
%
% Sub-function: read_settings.m
%
% Written by Karina L??viknes
% Written by Karina Loviknes
%

% Default values from settings file
[network, stations, first_day, last_day, channels, location,...
num_stat_cc, Fq, xaxis, yaxis, titl, bpf, lag_red, datesm] =...
num_stat_cc, Fq, xaxis, yaxis, titl, bpf, lag_red, datesm, fit] =...
read_settings(settingsfile, 'PLOT');

validateattributes(stations,{'cell'},{'nonempty'});
Expand Down Expand Up @@ -106,14 +106,15 @@
file2 = load(filename2);
timedelay = file2.timedelay.timedelay;
timedelay0 = file2.timedelay.timedelay0;
linear_td = file2.timedelay.linear_td;
linear_td0 = file2.timedelay.linear_td;
ref = file2.timedelay.reference;
else
warning(['Cannot find a mat.file with a measured time ',...
'delay for stationpair ' pair '. Fileformat must ',...
'be: TD_' pair '_' dates2 '.mat' ])
end


% Reduce computional effort by only using +-lag_red time lag
zerolag = find(lag==0);
egf = EGF(:,zerolag-lag_red:zerolag+lag_red);
Expand Down Expand Up @@ -176,8 +177,13 @@
end

subplot(2,2,4)
plot(dd1,linear_td/Fq,dd1,timedelay0/Fq,'.')
legend('Linearly fitted timedelay','Measured timedelay')
if strcmp(fit,'linfit')
plot(dd1,linear_td/Fq,dd1,timedelay0/Fq,'.')
legend('Linearly fitted timedelay','Measured timedelay')
else
plot(dd1,timedelay/Fq,dd1,timedelay0/Fq)
legend('Continous timedelay','Measured timedelay')
end
axis([1 num_days yaxis])
xlabel('Days','FontSize', 16), ylabel('Time delay (s)',...
'FontSize', 16)
Expand Down Expand Up @@ -224,9 +230,14 @@
'timedelays'],'FontSize', 15)
end

subplot(2,2,4)
plot(dd1,linear_td/Fq,dd1,timedelay0/Fq,'.')
legend('Linearly fitted timedelay','Measured timedelay')
subplot(2,2,4)
if strcmp(fit,'linfit')
plot(dd1,linear_td/Fq,'b',dd1,timedelay0/Fq,'r.')
legend('Linearly fitted timedelay','Measured timedelay')
else
plot(dd1,timedelay/Fq,dd1,timedelay0/Fq)
legend('Continous timedelay','Measured timedelay')
end
axis([1 num_days yaxis])
xlabel('Days','FontSize', 16), ylabel('Time delay (s)',...
'FontSize', 16)
Expand Down Expand Up @@ -277,8 +288,13 @@
end

subplot(2,1,2)
plot(dd1,linear_td/Fq,dd1,timedelay0/Fq,'.')
legend('Linearly fitted timedelay','Measured timedelay')
if strcmp(fit,'linfit')
plot(dd1,linear_td/Fq,'b',dd1,timedelay0/Fq,'r.')
legend('Linearly fitted timedelay','Measured timedelay')
else
plot(dd1,timedelay/Fq,dd1,timedelay0/Fq)
legend('Continous timedelay','Measured timedelay')
end
axis([1 num_days yaxis])
xlabel('Days','FontSize', 16), ylabel('Time delay (s)',...
'FontSize', 16)
Expand Down Expand Up @@ -363,9 +379,13 @@
end

subplot(2,nsp,sp+nsp)
plot(dd1,linear_td/Fq,'b-',dd1,timedelay0/Fq,'r.')
legend('Linearly fitted timedelay',...
'Measured timedelay')
if strcmp(fit,'linfit')
plot(dd1,linear_td/Fq,'b',dd1,timedelay0/Fq,'r.')
legend('Linearly fitted timedelay','Measured timedelay')
else
plot(dd1,timedelay/Fq,dd1,timedelay0/Fq)
legend('Continous timedelay','Measured timedelay')
end
axis([1 num_days yaxis])
xlabel('Days','FontSize', 10), ylabel(...
'Time delay (s)','FontSize', 10)
Expand Down Expand Up @@ -396,9 +416,13 @@
end

subplot(2,nch,spp+nch)
plot(dd1,linear_td/Fq,'b-',dd1,timedelay0/Fq,'r.')
legend('Linearly fitted timedelay',...
'Measured timedelay')
if strcmp(fit,'linfit')
plot(dd1,linear_td/Fq,'b',dd1,timedelay0/Fq,'r.')
legend('Linearly fitted timedelay','Measured timedelay')
else
plot(dd1,timedelay/Fq,dd1,timedelay0/Fq)
legend('Continous timedelay','Measured timedelay')
end
axis([1 num_days yaxis])
xlabel('Days','FontSize', 10), ylabel(...
'Time delay (s)','FontSize', 10)
Expand Down
1 change: 1 addition & 0 deletions EGF_tc_toolbox/read_settings.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -267,6 +267,7 @@
varargout{4}=filterp;
varargout{5}=lag_red;
varargout{6}=datesm;
varargout{7}=fit;

elseif strcmp(varargin{1},'FILTER')
% Read the variables used to filter
Expand Down
20 changes: 10 additions & 10 deletions EGF_tc_toolbox/settingsfile.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,15 +18,15 @@ Fq = 10

% The filename format, must include station name, date and folder (if in
% another folder), if not specified default is [network.stationname.location.HHchannels.D.yyyy.ddd.SAC]
filename = stationname/network.stationname.location.HH.D.(DATE).mseed
filename = ../stationname/network.stationname.location.HH.D.(DATE).mseed

% Input file format, can be either sac (default) or miniseed
%fileformat = sac
fileformat = miniseed

% The date format used in the filename, default is 'yyyy-mm-dd'
% dateformat = yyyy.mm.dd
dateformat = yyyy.ddd
dateformat = yyyymmdd

%% EGF:
% Specify the values specifically for loading the files and for estimating
Expand Down Expand Up @@ -59,31 +59,31 @@ wl = 1
swl = 24

% Cross correlation overlap window length in %, the default is 100%
perco = 50
%perco = 50

%% TD:
% Specify the values for measuring time shifts
%datesm = 2015-09-01 2015-11-01

% Bandpass filter to apply on the daily noise correlations and reference trace
%before measuring time delays, if not specified no filter will be applied
%filterm = [0.01 1.25]
filterm = [0.1429 0.5]

% Number of iterations to run the measuring process, must be an integer,
% default is 3
iterations = 3

% Only use +-lag_red time lag to reduce computational effect, default is
% 2000 (in samples)
lag_red = 500
lag_red = 2000

% Specification about the time period the reference is stacked over. The
% options are the whole period (default), months (the first and last day of
% the correlation period must be specified), surrounding_days, firstdays
% (number of days must be specified), and increase (for using the signal of
% the first day as reference in the first iteration and then a specified
% numdays for the next iteration):
stackperiod = increasing 1 3
stackperiod = firstdays 30

% Specification about which part of the signal should be used to measure
% the timeshift, the options are all (default), separated, whole, positive
Expand Down Expand Up @@ -111,14 +111,14 @@ reference_clock_station = N2ST

%% PLOT:
% Specify the values for plotting:
xaxis = -4 4
yaxis = -4 4
xaxis = -150 150
yaxis = -150 150

% Specify wheter the title should be the name of the plot ('name',default) or a letter (alphabet):
titl = alphabet
%titl = alphabet

% Filter
%filterp = [10 25]
%filterp = [0.1429 0.5]

%% FILTER
% Specify values for applying different bandpass filters:
Expand Down

0 comments on commit 800f9ef

Please sign in to comment.