%%
clear all;
start_time = [171 21 0 0]; % [jday hour min sec] of start
window = 3600*3; % 2 hours at 3600
das = 'A881'; % name of seismic station to read in

[seisTBAG,hdr] = segycut(start_time,window,das,'1');

%tvec=1:(1/hdr.sps):(853999/hdr.sps);

tvec=linspace(0,3,hdr.nsamp);
%% convert data to real physical stuff

%ts = (1:length(seisTBAG))/hdr.sps; % create time series vector
seis_sens = 800; % sensitivity is 800 V/m/s
atod = hdr.atod; % digitization V/count
sens_gain = hdr.gain; % programmed gain setting
seis_calib = atod/seis_sens/sens_gain*1e6;
%% apply butter filter

sps = hdr.sps; % get sample rate from header
npoles = 2; % specify order of filter
low = 10; % corner frequency at 10 Hz
high = 49; % corner frequency at 20 Hz
%%%%% signal processing tool box with butter is necessary for following operations

[B1,A1] = butter(npoles,[low high]/(sps/2));%bandpass filter
% % now do highpass
% [B2,A2] = butter(npoles,[high]/(sps/2),'high');
% % now do lowpass
% [B3,A3] = butter(npoles,[low]/(sps/2),'low');

%%
figure(3)
subplot(2,1,1) % plot unfiltered velocity trace
plot(tvec,seisTBAG,'k')
title(['das: ' das ' start time: ' num2str(start_time)])
ylabel('velocity (\mum/s)')

subplot(2,1,2)
bandpassV = filter(B1,A1,seisTBAG);
plot(tvec,bandpassV,'b')
ylabel('velocity (\mum/s)')
title(['bandpass filtered ' num2str(low) ' to ' num2str(high) ' Hz'])

% subplot(4,1,3)
% highpassV = filter(B2,A2,seisTBAG);
% plot(tvec,highpassV,'g')
% ylabel('velocity (\mum/s)')
% title(['highpass filtered above ' num2str(high) ' Hz'])
% 
% subplot(4,1,4)
% lowpassV = filter(B3,A3,seisTBAG);
% plot(tvec,lowpassV,'r')
% ylabel('velocity (\mum/s)')
% title(['bandpass filtered below ' num2str(low) ' Hz'])
xlabel('time (s)')

%% make spectrogram

data=bandpassV;
window=[];%this should be a factor of 2 number
overlap=0;%this should be a factor of 2 number
nfft=1024;
sps=hdr.sps;
figure(4)
spectrogram(data,window,overlap,nfft,sps,'yaxis');
axis xy; axis tight; colormap(jet); view(0,90)