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. Plots resulting waveforms

6. Save waveforms to SAC format

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

%% 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