% Simple demo code to illustrate process of reading and filtering a set of
% SAC files.
%
% Set your MATLAB working directory to a directory containing SAC files.
% This code will read in the SAC files, and tape, bandpass, and decimate
% the files.
%
% Created for the USArray Data Processing and Analysis Short Course
% August 15-19, 2011, Northwester University
% author: % Bob Woodward, IRIS
%
clear
%
% Use wildcards to build a list of files that we want to read
dirlist = dir('*LHZ*SAC');
fprintf('Number of SAC files selected: %d\n',length(dirlist));
%
% loop over the selected files and only read the headers
%
maxpts = 0;
tic
for i=1:length(dirlist)
    s = rsac_h(dirlist(i).name);
    maxpts = max(maxpts,s.npts);
    sachdr(:,i) = s(:);
end;
toc
fprintf('scanned file headers, maxpts: %d\n',maxpts);
%
% preallocate memory for speed
sacdat = zeros(maxpts,length(dirlist));
lat    = zeros(length(dirlist),1);
lon    = zeros(length(dirlist),1);
tbeg   = zeros(length(dirlist),1);
npts   = zeros(length(dirlist),1);
dt     = zeros(length(dirlist),1);
tdiff  = zeros(length(dirlist),1);
%
% read each timeseries and build a single matrix containing all the time series
%
tic
for i=1:length(dirlist)
    s = rsac(dirlist(i).name); % this reads the sac file into a structure
    % save some important metadata
    npts(i) = s.npts; % number of points in the timeseries
    tbeg(i) = datenum(s.nzyear, 0, s.nzjday, s.nzhour, s.nzmin, s.nzsec+s.nzmsec/1000.) + s.b; % start time of time series
    dt(i)   = s.delta; % sample interval, in seconds
    lat(i)  = s.stla;  % station latitude
    lon(i)  = s.stlo;  % station longitude
    % add this seismogram as next column of matrix
    sacdat(1:npts(i),i) = s.d(1:npts(i)); % put this seismogram into a column of the matrix
end
toc
fprintf('read data\n');
%
% pre-processing loop
%
for i=1:length(lat)
    sacdat(1:npts(i),i) = sacdat(1:npts(i),i) - mean(sacdat(1:npts(i),i)); % demean the data
    t = tukeywin(npts(i),0.05); % compute a 5% taper
    sacdat(1:npts(i),i) = t .* sacdat(1:npts(i),i); % apply the taper
    %
    % plot some intermediate results
    %
    clf
    hold on
    subplot(2,1,1);
    plot (sacdat(1:npts(i),i));  % plot raw trace
    
    nq = 0.5 / dt(i); % Nyquist frequency
    % apply high pass filter
    co = 0.004 / nq; % normalized freq, in Hz
    [b a] = butter(4,co,'high');
    sacdat(1:npts(i),i) = filter(b,a,sacdat(1:npts(i),i));
    % apply low pass filter
    co = 0.02 / nq; % normalized freq, in Hz
    [b a] = butter(4,co,'low');
    sacdat(1:npts(i),i) = filter(b,a,sacdat(1:npts(i),i));
    %
    % plot filtered trace
    
    subplot(2,1,2);
    plot (sacdat(1:npts(i),i));
    
    fprintf('%s\n',sachdr(i).kstnm);
    pause
end;
save mydata sacdat lat lon tbeg npts dt