% Simple code to demonstrate using instrument response computations
% within MATLAB.
%

clear
%
% read a SAC file
s = rsac('2009.215.17.53.10.0000.TA.B17A..LHZ.R.SAC');
%
% deconvolve instrument response between 0.004 Hz and 0.050 Hz
p = transfer(s.d, s.npts, 1/s.delta, 'SAC_PZs_TA_B17A_LHZ__2007.319.00.00.00.0000_99999.9999.24.60.60.99999', .004, .050);

%
% alternatively, filter the raw data
t = tukeywin(s.npts,0.05); % compute a 5% taper
r = t .* s.d; % apply the tape
nq = 0.5 / s.delta; % Nyquist frequency
% apply high pass filter
co = 0.004 / nq; % normalized freq, in Hz
[b a] = butter(4,co,'high');
r = filter(b,a,r);
% apply low pass filter
co = 0.05 / nq; % normalized freq, in Hz
[b a] = butter(4,co,'low');
r = filter(b,a,r);
%
% normalize the traces, to simplify overplotting
s.d = s.d / max(abs(s.d));
p = p / max(abs(p));
r = r / max(abs(r));
%
% plot a comparison of deconvolved data and filtered data
clf
hold on
% plot(s.d) % raw data
plot(r) % filtered data
plot(p,'r')     % deconvolved, filtered data


%
% Some test code for reading in a pole-zero file and evaluating a respone
% 
% [zz pp constant] = read_sac_pole_zero('SAC_PZs_TA_B17A_LHZ__2007.319.00.00.00.0000_99999.9999.24.60.60.99999');
% 
% compute response between 0 to 200 Hz
% for j=1:s.npts+1
%     freq(j) = (j-1) * 200 / s.npts;
%     b(j)=generate_response(zz,pp,constant,freq(j));
% end
