-
Notifications
You must be signed in to change notification settings - Fork 17
/
a5_plot_real_imag.m
170 lines (146 loc) · 6.42 KB
/
a5_plot_real_imag.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
% Plot the real and imaginary parts of the cross spectra to evaluate where signal
% is best and degree of bias from inhomogeneous noise sources.
%
% https://github.com/jbrussell
clear;
setup_parameters;
%======================= PARAMETERS =======================%
comp = 'ZZ'; %'ZZ'; %'RR'; %'TT';
period_lims = [2 100];
windir = 'window3hr';
pts_smooth = 20; % just for plotting purposes
issemilogx = 0;
%==========================================================%
stalist = parameters.stalist;
nsta = parameters.nsta;
winlength = parameters.winlength;
figpath = parameters.figpath;
%fig_winlength_path = [figpath,'window',num2str(winlength),'hr/fullStack/'];
% custom directory names
fig_winlength_path = [figpath,windir,'/fullStack/'];
%------------ PATH INFORMATION -------------%
ccf_path = parameters.ccfpath;
%ccf_winlength_path = [ccf_path,'window',num2str(winlength),'hr/'];
ccf_winlength_path = [ccf_path,windir,'/'];
ccf_singlestack_path = [ccf_winlength_path,'single/'];
ccf_daystack_path = [ccf_winlength_path,'dayStack/'];
ccf_monthstack_path = [ccf_winlength_path,'monthStack/'];
ccf_fullstack_path = [ccf_winlength_path,'fullStack/'];
ccf_stack_path = ccf_fullstack_path;
figpath = [fig_winlength_path,'bessel/'];
% create figure directory
if ~exist(fig_winlength_path)
mkdir(fig_winlength_path)
end
if ~exist(figpath)
mkdir(figpath)
end
ccf_path = [ccf_stack_path,'ccf',comp,'/',];
npairall = 0;
sum_all_psd = 0;
sum_all_psd_log = 0;
%------------ LOAD DATA AND CALCULATE AND PLOT PSD -------------%
for ista1=1:nsta % loop over all stations
sta1=char(stalist(ista1,:));
sta1dir=[ccf_path,sta1]; % dir to have all cross terms about this central station
nstapair = 0;
for ista2 = 1: nsta % loop over station pairs
sta2 = char(stalist(ista2,:));
% if same station, skip
if(strcmp(sta1,sta2))
continue
end
filename = sprintf('%s/%s_%s_f.mat',sta1dir,sta1,sta2);
if ~exist(filename,'file') % check that ccf file exists
disp(['not exist ',filename])
continue;
end
nstapair = nstapair + 1;
%----------- LOAD DATA -------------%
data = load(filename);
dt = data.stapairsinfo.dt;
ccf = data.coh_sum./data.coh_num;
ccfreal = smooth(real(ccf),pts_smooth);
ccfimag = smooth(imag(ccf),pts_smooth);
ccfamp = smooth(abs(ccf),pts_smooth);
ccf_win = data.coh_sum_win./data.coh_num;
ccfreal_win = smooth(real(ccf_win),1);
ccfimag_win = smooth(imag(ccf_win),1);
ccfamp_win = smooth(abs(ccf_win),1);
%------------------------%
% Distance between sta1 and sta2
sta1sta2_dist(nstapair) = deg2km(distance(data.stapairsinfo.lats(1),data.stapairsinfo.lons(1),data.stapairsinfo.lats(2),data.stapairsinfo.lons(2)));
% Check if reverse station pair has already been plotted
stapairinv = [sta2,'_',sta1];
if exist('existpair','var')
if find(strncmp(stapairinv,existpair,length(stapairinv)))
continue
end
end
% Update some other useful variables
dumsta2{nstapair} = sta2;
npairall = npairall + 1; % number of total station pairs
sta1sta2_dist_all(npairall) = sta1sta2_dist(nstapair); % vector containing distance between each station pair
existpair(npairall) = {[sta1,'_',sta2]};
%----------- PLOT REAL AND IMAGINARY -------------%
f201 = figure(201); clf; set(gcf, 'Color', 'w'); box on;
set(gcf,'position',[308 115 600 590]);
subplot(3,1,1)
T = length(ccfreal);
faxis = [0:(T-mod(T-1,2))/2 , -(T-mod(T,2))/2:-1]/dt/T;
ind = find(faxis>0);
if issemilogx
h(1) = semilogx((faxis(ind)),ccfreal(ind),'-','color',[0.5, 0.5, 0.5],'linewidth',2); hold on;
h(2) = semilogx((faxis(ind)),ccfreal_win(ind),'-k','linewidth',2);
else
h(1) = plot((faxis(ind)),ccfreal(ind),'-','color',[0.5, 0.5, 0.5],'linewidth',2); hold on;
h(2) = plot((faxis(ind)),ccfreal_win(ind),'-k','linewidth',2);
end
title(sprintf('REAL %s %s x-spectra %s ,distance: %f km',sta1,sta2,comp(1),sta1sta2_dist(nstapair)));
legend(h,{'full data';'windowed'});
xlim([1/period_lims(2) 1/period_lims(1)]);
%xlim([0.04 0.16])
xlabel('Frequency');
ylim([-0.03 0.03]);
%ylim([-0.01 0.01]);
set(gca,'linewidth',1.5);
subplot(3,1,2)
if issemilogx
h(1) = semilogx((faxis(ind)),ccfimag(ind),'-','color',[0.5, 0.5, 0.5],'linewidth',2); hold on;
h(2) = semilogx((faxis(ind)),ccfimag_win(ind),'-k','linewidth',2);
else
h(1) = plot((faxis(ind)),ccfimag(ind),'-','color',[0.5, 0.5, 0.5],'linewidth',2); hold on;
h(2) = plot((faxis(ind)),ccfimag_win(ind),'-k','linewidth',2);
end
title(sprintf('IMAGINARY %s %s x-spectra %s ,distance: %f km',sta1,sta2,comp(1),sta1sta2_dist(nstapair)));
xlim([1/period_lims(2) 1/period_lims(1)]);
%xlim([0.04 0.16])
xlabel('Frequency');
ylim([-0.03 0.03]);
%ylim([-0.01 0.01]);
set(gca,'linewidth',1.5);
subplot(3,1,3)
if issemilogx
h(1) = semilogx((faxis(ind)),ccfamp(ind),'-','color',[0.5, 0.5, 0.5],'linewidth',2); hold on;
h(2) = semilogx((faxis(ind)),ccfamp_win(ind),'-k','linewidth',2);
else
h(1) = plot((faxis(ind)),ccfamp(ind),'-','color',[0.5, 0.5, 0.5],'linewidth',2); hold on;
h(2) = plot((faxis(ind)),ccfamp_win(ind),'-k','linewidth',2);
end
title(sprintf('AMP %s %s x-spectra %s ,distance: %f km',sta1,sta2,comp(1),sta1sta2_dist(nstapair)));
xlim([1/period_lims(2) 1/period_lims(1)]);
%xlim([0.04 0.16])
xlabel('Frequency');
%ylim([-0.03 0.03]);
set(gca,'linewidth',1.5);
drawnow;
% pause;
%print(f201,'-dpdf',[figpath,'bessel_',comp,'_',sta1,'_',sta2,'.pdf']); % Save figure
if issemilogx
save2pdf([figpath,'log_bessel_',comp,'_',sta1,'_',sta2,'_',num2str(pts_smooth),'pts_comparewin.pdf'],f201,1000);
else
save2pdf([figpath,'bessel_',comp,'_',sta1,'_',sta2,'_',num2str(pts_smooth),'pts_comparewin.pdf'],f201,1000);
end
end % ista2
end % ista1
%pause;