%% drumplot cookbook
% In this example you will learn how to create a waveform object by loading
% a SAC file, extract 1 hour of data from that waveform object, and make
% a helical drum recorder plot using the drumplot class. 
% 
% Another feature of drumplot is that it allows detected events, from a
% Catalog object, to be superimposed on the helicorder plot. In the example
% below this Catalog is generated by running an STA/LTA detection algorithm
% on the extracted 1-hour waveform object.

%% Load 1 day of data for REF.EHZ on 2009.03.22

% here we can either use scnlobject or ChannelTag to describe the
% net-sta-loc-chan we want. waveform understands both.
ctag = ChannelTag('AV.REF.-.EHZ')

% starttime in MATLAB datenum format
snum = datenum(2009,3,22);

% endtime (1 day after start)
enum=snum+1;

% create the datasource object - in this case it is a SAC file
ds = datasource('sac', 'SACDATA/REF.EHZ.2009-03-22.sac');

% create the waveform object with this datasource, ChannelTag, starttime
% and endtime
w=waveform(ds, ctag, snum, enum);


%% fill gaps, detrend, band pass filter

% in case there are gaps in the time series (marked by NaN) we can interpolate a
% meaningful value
w = fillgaps(w, 'interp');

% detrend the time series - this removes linear drift 
w = detrend(w);

% create a Butterworth bandpass filter from 0.5 to 15 Hz, 2 poles
fobj = filterobject('b', [0.5 15], 2);

% apply the filter in both directions (acausal) - this is a zero phase
% filter which is helpful because it doesn't disperse different frequency 
% components. the caveat is that it can spread energy so arrivals may appear
% slighter earlier than they actually are
w = filtfilt(fobj, w);

%% plot the waveform object
figure
plot(w)

%% extract the first 1 hour of data and plot with 5 minutes per line
starttime = get(w,'start');
w2=extract(w, 'time', starttime, starttime + 1/24 )

% make a drumplot object. mpl means minutes per line and is set here to 5.
h2 = drumplot(w2, 'mpl', 5);

% plot the drumplot object - many events are visible, this is an earthquake
% swarm
plot(h2)

%% Run an STA/LTA detector on the waveform object to detect events
% Compute STA/LTA ratio. In this example, trigger "ON" when STA/LTA exceeds
% 3, trigger "OFF" when this drops back below 1.5. If the trigger was on
% for at least 2 seconds, declare it as an event. 
% All the events get saved into a Catalog object.

% set the STA/LTA detector
sta_seconds = 0.7; % STA time window 0.7 seconds
lta_seconds = 7.0; % LTA time window 7 seconds
thresh_on = 3; % Event triggers "ON" with STA/LTA ratio exceeds 3
thresh_off = 1.5; % Event triggers "OFF" when STA/LTA ratio drops below 1.5
minimum_event_duration_seconds = 2.0; % Trigger must be on at least 2 secs
pre_trigger_seconds = 0; % Do not pad before trigger
post_trigger_seconds = 0; % Do not pad after trigger
event_detection_params = [sta_seconds lta_seconds thresh_on thresh_off ...
    minimum_event_duration_seconds];
%%
% run the STA/LTA detector. lta_mode = 'frozen' means the LTA stops
% updating when trigger is "ON". 
[detectionObject,sta,lta,sta_to_lta] = Detection.sta_lta(w2, 'edp', event_detection_params, ...
    'lta_mode', 'frozen');

% not sure what this is for
set(gca, 'XLim', [44*60 48*60])

%% Plot detected events on top of the continuous drumplot
h3 = drumplot(w2, 'mpl', 5, 'detections', detectionObject);
plot(h3)

