-
Notifications
You must be signed in to change notification settings - Fork 5
/
tremometer_control.m
161 lines (131 loc) · 4.83 KB
/
tremometer_control.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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Tremometer 1.0 %
% %
% Matlab code to automatically detect and characterize harmonic tremor %
% in continuous seismic data using a pitch-detection algorithm %
% %
% Diana C. Roman %
% Last modified July 31, 2017 %
% %
% Please cite: Roman, D.C. (2017), Automated detection and characterization of harmonic tremor %
% in continuous seismic data. Geophys. Res. Lett.,44,doi: 10.1002/2017GL073715. %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set the following variables %
% (or leave as default values) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x=load('example.asc'); %Assign one day of (instrument-corrected and filtered) data to x
yr=2015; %Set start date (four digit year) of input data
mo=03; %Set start date (two digit month) of input data
dy=22; %Set start date (two digit day) of input data
fs=40; %set sampling frequency (in Hz)
minfreq=0.5; %set the minimum frequency of the fundamental (in Hz). Default is 0.5 Hz
minHSI=30; %set the minimum Harmonic Strength Index for fundamental and two overtones. Default is 30.
calmodeflag=0; %Toggle periodogram plots on (1) or off (0). Default is on (1).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Don't change anything below here %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nHarm=4;
s=1440;
for i=1:s
startpt(i)=1+((i-1)*fs*60);
endpt(i)=i*fs*60;
end
final=zeros(14,s);
for j=1:s
[harmpow,harmfreq,interharmpow,harmstrength]=tremometer(x(startpt(j):endpt(j)),fs,nHarm);
final(:,j)=vertcat(harmfreq,harmpow,interharmpow,harmstrength);
end
final=final';
% Check minimum frequency criterion
for j=1:s
if final(j,1) > minfreq
final(j,15)=1;
else final(j,15)=0;
end
end
% Check minimum HSI criterion for each harmonic
for j=1:s
if final(j,12) > minHSI
final(j,16)=1;
else final(j,16)=0;
end
end
for j=1:s
if final(j,13) > minHSI
final(j,17)=1;
else final(j,17)=0;
end
end
for j=1:s
if final(j,14) > minHSI
final(j,18)=1;
else final(j,18)=0;
end
end
% Check if all criteria are met and flag tremor episodes, then add them to a list of detections with times
for j=1:s
if final(j,15)+final(j,16)+final(j,17)+final(j,18)==4
final(j,19)=1;
else final(j,19)=0;
end
end
%Make a table of harmonic tremor detections
ind1 = final(:,19) == 1;
date=datetime(yr,mo,dy,0,0,0);
t1=date;
t2=date+minutes(1439);
t=t1:minutes(1):t2;
t=t';
detection_times=t(ind1,:);
fund_freq=final(ind1,1);
harm_freq_1=final(ind1,2);
harm_freq_2=final(ind1,3);
harm_freq_3=final(ind1,4);
harm_pow_1=final(ind1,5);
harm_pow_2=final(ind1,6);
harm_pow_3=final(ind1,7);
harm_pow_4=final(ind1,8);
interharm_pow_1=final(ind1,9);
interharm_pow_2=final(ind1,10);
interharm_pow_3=final(ind1,11);
HSI_1=final(ind1,12);
HSI_2=final(ind1,13);
HSI_3=final(ind1,14);
clc
detections=table(detection_times,fund_freq,harm_freq_1,harm_freq_2,harm_freq_3,HSI_1,HSI_2,HSI_3)
% Make some plots
plot(t,final(:,12))
datetick('x')
ylabel('Harmonic Strength Index')
xlabel('Time of Day')
%hline=refline([0 minHSI]);
%hline.Color = 'r';
yline(minHSI);
%end
figure
plot(detection_times,harm_freq_1, 'o')
ylabel('Frequency (Hz)')
xlabel('Time of Day')
%end
%Calibration mode - make periodogram plot for all detections
if calmodeflag==1;
for j=1:s
if final(j,19) == 1
figure
snip=x(startpt(j):endpt(j));
snip=snip-mean(snip);
n=length(snip);
w=kaiser(n,38);
[Pxx,F]=periodogram(snip,w,n,fs);
plot(F,Pxx)
hold('on')
plot(final(j,1:3),final(j,5:7),'*')
ylabel('Power')
xlabel('Frequency (Hz)')
end
end
end
%tidy up
clear calmodeflag date detection_times dy endpt F fs fund_freq harm_freq_1 harm_freq_2 harm_freq_3 harm_pow_1 harm_pow_2 harm_pow_3 harm_pow_4 harmfreq harmpow harmstrength hline HSI_1 HSI_2 HSI_3 i ind1 interharmpow interharm_pow_1 interharm_pow_2 interharm_pow_3 j minfreq minHSI mo n nHarm prcdata Pxx s samplrate snip startpt t t1 t2 temp_detections w x yr