Downsample and deconvolve instrument response

In this example, requested by Celso Alvizuri, we:

1. Load data from the Alaska Earthquake Center (University of Alaska Fairbanks) waveform database

2. Clean the data to remove spikes, gaps and trend

3. Downsample to a common sampling frequency

4. Deconvolve the instrument response (pads the data, applies Tukey window, filters, FFT, deconvolve, IFFT, unpad)

5. Plot resulting waveforms

6. Save waveforms to SAC format

%% Specify Master stations db
dbMaster = '/aerun/sum/db/dbmaster/master_stations';

%% Load data
startTime = datenum(2009,2,15,19,33,20);
endTime = datenum(2009,2,15,19,40,0);
ds = datasource('uaf_continuous');
nslc(1) = scnlobject('COLA','BHZ'); % Has sampling rate of 20 Hz
nslc(2) = scnlobject('RSO','EHZ'); % Has sampling rate of 100 Hz
w=waveform(ds, nslc, startTime, endTime);

%% Clean data
w = medfilt1(w,3); % remove any spikes of sample length 1
w = fillgaps(w, 'interp');
w = detrend(w);

%% Downsample to a common sampling frequency
target_fsamp = 20; % samples per second
for c=1:numel(w)
    current_fsamp = round(get(w(c),'freq'));
    w(c) = resample(w(c), 'mean', current_fsamp/target_fsamp );
end

%% Deconvolve instrument response from AEC Master stations database
% response_apply was written by Mike
filterObj = filterobject('b',[0.5 5],2);
wFilt = filtfilt(filterObj,w);
wCorrected = response_apply(wFilt,filterObj,'antelope',dbMaster);

%% Plot result
plot_panels(wCorrected)

%% Save as SAC
for c=1:numel(w)
    sta=nslc(c).station;
    chan=nslc(c).channel;
    fname=sprintf('%s.%s.sac',sta,chan);
    savesac(w(c),'.',fname);
end

After this the files can be loaded, plotted and headers listed with sac:

> sac

SAC> r COLA.BHZ.sac

SAC> pl

SAC> lh

(Glenn Thompson, 2016/05/19)