% Spectrograms
% Script written by Carolyn Parcheta and Brando
% June 23, 2009 NMTech Course

%% Objective: pre-code
% 1. Load Data
% 2. Convert from segy to matlab
% 3. Convert to digital and then counts to volts, etc
% 4. detrend (removes mean)
% 5. Butter filter bandpass 0.01-25 Hz
% 6. spectrogram(input, window) = will reshape (partition data in one hour
% long windows) and fft (fast fournier transform)
% 7. have some fun after!

%% READ IN DATA
cd '/Users/geop572/Desktop/SEGY/' % move to data directory

start_time = [171 21 0 0]; % [jday hour min sec]
window = 60*60*3; % seconds of data
das = '9026'; % name of seismic station datalogger to read in. das = datalogger = 9026
Station_Name = 'TLAV';

[wiggleV,hdr] = segycut(start_time,window,das,'1'); % start_time,window,das,channel
[wiggleN,hdr] = segycut(start_time,window,das,'2');
[wiggleE,hdr] = segycut(start_time,window,das,'3');

sps = hdr.sps;
ts = (1:length(wiggleV))/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

%% Detrend/ Remove Mean

wiggleV_detrend = detrend(wiggleV); 
wiggleN_detrend = detrend(wiggleN);
wiggleE_detrend = detrend(wiggleE);

% detrend forces signal to oscillate over zero. (e.g., infrasound sensor
% oscillates over 10 millivolts and not zero

%% Bandpass Filter 0.01 - 49.99 Hz

npoles = 2;
low = 49.99; % corner frequency at 49.99 Hz. You cannot use 50 Hz b/c sps is 100Hz and you cannot get info above 49.99 Hz. (100Hz/2 = threshold)
high = 0.01; % corner frequency at 1 Hz
% signal processing tool box with butter is necessary for following operations
% first do bandpass

[B1,A1] = butter(npoles,[high low]/(sps/2)); % The order of high and low MATTERS. lowest number must come first.
% now do highpass

bandpass_wiggleV = filter(B1,A1,wiggleV);
bandpass_wiggleN = filter(B1,A1,wiggleN);
bandpass_wiggleE = filter(B1,A1,wiggleE);


%% Spectrograms
nfft = 512; 

figure(1); clf

subplot(3,2,1)
spectrogram(bandpass_wiggleV, 10000 , [] ,nfft,sps,'yaxis');
Ylabel('Hz')
Xlabel('time(s)')
colorbar
title(['TLAV Vertical Component' '  | start time: ' num2str(start_time)])

subplot(3,2,2)
plot(ts,bandpass_wiggleV*seis_calib)
ylabel('vertical velocity (\mum/s)')
Xlabel('time(s)')
title(['TLAV Vertical Component' ' | start time: ' num2str(start_time)])

subplot(3,2,3)
spectrogram(bandpass_wiggleN, 10000 , [] ,nfft,sps,'yaxis');
Ylabel('Hz')
Xlabel('time(s)')
colorbar
title(['TLAV North Component' '  | start time: ' num2str(start_time)])

subplot(3,2,4)
plot(ts,bandpass_wiggleN*seis_calib)
ylabel('N/S velocity (\mum/s)')
Xlabel('time(s)')
title(['TLAV North Component' ' | start time: ' num2str(start_time)])

subplot(3,2,5)
spectrogram(bandpass_wiggleE, 10000 , [] ,nfft,sps,'yaxis');
Ylabel('Hz')
Xlabel('time(s)')
colorbar
title(['TLAV East Component' '  | start time: ' num2str(start_time)])

subplot(3,2,6)
plot(ts,bandpass_wiggleE*seis_calib)
ylabel('E/W velocity (\mum/s)')
xlabel('time (s)')
title(['TLAV East Component' ' | start time: ' num2str(start_time)])

% Couldn't get this to work:
% [Sv,Fv,Tv,Pv] = spectrogram(bandpass_wiggleV, 10000, 200, 512, 100);
% pcolor(Tv,Fv,Pv);
% shading interp
% So we went with this:
% spectrogram(bandpass_wiggleV, 10000 , [] ,nfft,sps,'yaxis');