Downsample and deconvolve instrument response - geoscience-community-codes/GISMO GitHub Wiki

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