%% plot 24 hour data time series and helicorder record
%%%%%%%%%%%%%%%%%%%%%
% Used in June 2009 for NM Tech Field Course



%% read in 24 hrs data
%%%%%%%%%%%%%%%%%%%%%

cd '/Users/geop572/Desktop/GEOP572_MATLAB/Kilauea2008/SEGY/' % move to data directory

start_time = [187 0 0 0]; % [jday hour min sec]
window = 86400; % 24 hours or 86400 s of data
das = '92C8'; % name of seismic station to read in. das = datalogger (92C8 and 9919 for R186.01)

[seisV,hdr] = segycut(start_time,window,das,'1'); % start_time,window,das,channel

ts = (1:length(seisV))/hdr.sps; % create time series vector
seis_sens = 800; % sensitivity of Guralp T40 is 800 V/m/s
atod = hdr.atod; % Anolog to digital. converts counts to volts. 
sens_gain = hdr.gain; % programmed gain setting
seis_calib = atod/seis_sens/sens_gain*1e6; % total conversion of counts to um/s

%% now plot the calibrated data trace for 24 hours
figure(1); clf
decimation = 10; %resamples data as not to include all data points. 
plot(ts(1:decimation:end)/60/60,seisV(1:decimation:end)*seis_calib,'k') % plot decimated data trace for 24 hours (k = blck)
ylabel('velocity (\mum/s)')
xlabel('time (hours)')
title(['das: ' das ' start time: ' num2str(start_time)])


%% now plot a traditional helicorder record of the same data
figure(2); clf
decimation = 2;
ts_jday = hdr.day + ts/60/60/24; % time axis for the ogram script must be decimal julian day
scaling = 8; % scaling is set to the window, larger numbers amplify the traces by more
ogram(ts_jday(1:decimation:end),seisV(1:decimation:end),scaling) % this is J. Johnson's script
ylabel('time(hours)')
xlabel('minutes')














%% 3 axis data, convert vertical data to disp. and acc., and do spectra 
%%%%%%%%%%%%%%%%%
% read in ten minutes and 3 channels of data for a single das
%%%%%%%%%%%%%%%%%

start_time = [187 0 0 0]; % [jday hour min sec] of start
window = 600; % read in ten minutes of data
das = '92C8'; % name of seismic station to read in

[seisV,hdr] = segycut(start_time,window,das,'1');
[seisN,hdr] = segycut(start_time,window,das,'2');
[seisE,hdr] = segycut(start_time,window,das,'3');

ts = (1:length(seisV))/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; % total conversion of counts to um/s

%% now plot the three components on a single figure
figure(3)
subplot(3,1,1)
plot(ts,seisV*seis_calib)
ylabel('vertical velocity (\mum/s)')
title(['das: ' das ' start time: ' num2str(start_time)])

subplot(3,1,2)
plot(ts,seisN*seis_calib)
ylabel('N/S velocity (\mum/s)')

subplot(3,1,3)
plot(ts,seisE*seis_calib)
ylabel('E/W velocity (\mum/s)')
xlabel('time (s)')

%% now plot vertical velocity, displacement, and acceleration
figure(4)
subplot(3,1,1)
velV = seisV-mean(seisV); % get rid of DC offset
plot(ts,velV*seis_calib)
ylabel('velocity (\mum/s)')
title(['das: ' das ' start time: ' num2str(start_time)])
sps = hdr.sps; % sample rate of trace data

subplot(3,1,2)
dispV = detrend(cumsum(velV)/sps); % convert to displacement and get rid of trend
plot(ts,dispV*seis_calib)
ylabel('displacement (\mum)')

subplot(3,1,3)

accV = gradient(velV)*sps; % convert to acceleration by differentiating
plot(ts,accV*seis_calib)
ylabel('acceleration (\mum/s^2)')
xlabel('time (s)')

%% now plot overyling frequency spectra for velocity, displacement, and acceleration
figure(5); clf
nfft = 512; % fast fournier transform; 512 = 2^9 why mult of 2? b/c fft is more efficient
% Why not 2^8? 256 = 2^8 and we are sampling at 100Hz this gives 2.56 s.
% 512 allows us to see lower frequencies
velV_spectra = abs(fft(velV,nfft));
dispV_spectra = abs(fft(dispV,nfft));
accV_spectra = abs(fft(accV,nfft));
fs = (1:length(velV_spectra))/length(velV_spectra)*sps; % fs = sampling freq

loglog(fs,velV_spectra,fs,dispV_spectra,fs,accV_spectra)
legend('velocity','displacement','acceleration',3)
xlim([0 50])
grid on








%% apply filters to data and plot their traces and spectra
%%%%%%%%%%%%%%%%%
% read in ten minutes of data, single component and filter data
%%%%%%%%%%%%%%%%%
start_time = [187 0 0 0]; % [jday hour min sec] of start
window = 600; % read in ten minutes of data
das = '92C8'; % name of seismic station to read in

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

ts = (1:length(seisV))/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; % total conversion of counts to um/s

%% now plot data filtered traces in several bands
figure(6); clf

sps = hdr.sps; % get sample rate from header
npoles = 2; % specify order of filter
low = 10; % corner frequency at 10 Hz
high = 20; % corner frequency at 20 Hz
%%%%% signal processing tool box with butter is necessary for following operations
% first do bandpass
[B1,A1] = butter(npoles,[low high]/(sps/2))
% now do highpass
[B2,A2] = butter(npoles,[high]/(sps/2),'high')
% now do lowpass
[B3,A3] = butter(npoles,[low]/(sps/2),'low')

ts = (1:length(seisV))/sps;

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

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

subplot(3,2,4)
highpassV = filter(B2,A2,seisV);
plot(ts,highpassV,'g')
ylabel('velocity (\mum/s)')
title(['highpass filtered above ' num2str(high) ' Hz'])

subplot(3,2,6)
lowpassV = filter(B3,A3,seisV);
plot(ts,lowpassV,'r')
ylabel('velocity (\mum/s)')
title(['bandpass filtered below ' num2str(low) ' Hz'])
xlabel('time (s)')


%% now plot overyling frequency spectra for filtered velocity
figure(7); clf
nfft = 512;
raw_spectra = abs(fft(seisV,nfft));
bandpass_spectra = abs(fft(bandpassV,nfft));
highpass_spectra = abs(fft(highpassV,nfft));
lowpass_spectra = abs(fft(lowpassV,nfft));
fs = (1:length(raw_spectra))/length(raw_spectra)*sps;

loglog(fs,raw_spectra,fs,bandpass_spectra,fs,highpass_spectra,fs,lowpass_spectra)
legend('raw velocity','bandpass','highpass','lowpass',3)
xlim([0 50])
grid on






%% read in multi-station data and display filtered traces
%%%%%%%%%%%%%%%%%
% read in ten minutes of data, single component and filter data
%%%%%%%%%%%%%%%%%
start_time = [187 0 0 0]; % [jday hour min sec] of start
window = 60; % read in one minute of data

das = '92C8'; % name of seismic station to read in
[seisVs(:,1),hdr] = segycut(start_time,window,das,'1');
das = '9919'; % name of seismic station to read in
[seisVs(:,2),hdr] = segycut(start_time,window,das,'1');


ts = (1:length(seisVs))/hdr.sps; % create time series vector


%% now plot raw data traces for 1 minute window
figure(8); clf

for k=1:2
    tmp_seistrace = seisVs(:,k)/(max(seisVs(:,k))-min(seisVs(:,k)));
    tmp_seistrace = tmp_seistrace - mean(tmp_seistrace);
    plot(ts,-k/1.5 + tmp_seistrace,'k')
    hold on
end
set(gca,'ytick',[]) % get rid of ordinate
title(['(unfiltered vertical seismogram) - start time: ' num2str(start_time)])
xlabel('time (s)')


axis tight


%% now plot low pass filtered data traces for 1 minute window
figure(9); clf

low = .125; % corner frequency at 8 s
high = .25; % corner frequency at 4 s
[B,A] = butter(npoles,[low high]/(sps/2))

for k=1:2
    tmp_seistrace = filter(B,A,seisVs(:,k))
    tmp_seistrace = tmp_seistrace/(max(tmp_seistrace)-min(tmp_seistrace));
    tmp_seistrace = tmp_seistrace - mean(tmp_seistrace);
    plot(ts,-k/1.5 + tmp_seistrace,'k')
    hold on
end
set(gca,'ytick',[]) % get rid of ordinate
title(['(lowpass filtered vertical seismograms) - start time: ' num2str(start_time)])
xlabel('time (s)')

axis tight


%%