%% plot 24 hour data time series and helicorder record %%%%%%%%%%%%%%%%%%%%% % read in 24 hrs data %%%%%%%%%%%%%%%%%%%%% clear all, close all %cd '/Users/Jeff/Documents/DATA/KILAUEA2008/SEGY' % move to data directory cd '/Users/andrea/Desktop/SEGY_tungu09' % day=166:172; %start_time = [168 00 00 00]; %start_time = [171 22 00 00]; % [jday hour min sec] of start % start_time=[168 04 19 53]; %1-3Hz Explosion % start_time=[168 08 49 05]; %1-3Hz Explosion % start_time=[168 19 34 48]; %1-3Hz Explosion % start_time=[169 13 52 12]; %2-30Hz Rockfall %start_time=[170 18 13 37]; %1-11Hz Lahar % start_time=[171 10 30 24]; %1-4Hz Explosion %start_time=[171 22 00 00]; %1-25Hz Lahar % start_time=[172 08 28 31]; %1.5-5Hz LP % start_time=[171 19 04 46]; %1-25Hz ?? % start_time=[167 23 01 35]; % Background Noise % start_time=[163 04 30 42]; %4-12Hz VT start_time=[167 18 28 07]; %1.5-5Hz LP window = 3600*5/60/60; % 24 hours or 86400 s of data dasTLAV = '9026'; % name of seismic station to read in dasTLEF = 'A871'; dasTMAR = 'A270'; dasTBAG = 'A881'; dasTGRS = 'A860'; dasTBUG = 'AA32'; dasTREX = '9FCC'; % [seisV(:,1),hdr] = segycut([166 0 0 0],window,das,'1'); % [seisV(:,2),hdr] = segycut([167 0 0 0],window,das,'1'); % [seisV(:,3),hdr] = segycut([168 0 0 0],window,das,'1'); % [seisV(:,4),hdr] = segycut([169 0 0 0],window,das,'1'); % [seisV(:,5),hdr] = segycut([170 0 0 0],window,das,'1'); % [seisV(:,6),hdr] = segycut([171 0 0 0],window,das,'1'); % [seisV(:,7),hdr] = segycut([172 0 0 0],window,das,'1'); [seisV_TLAV,hdr_TLAV] = segycut(start_time,window,dasTLAV,'1'); [seisV_TLEF,hdr_TLEF] = segycut(start_time,window,dasTLEF,'1'); [seisV_TMAR,hdr_TMAR] = segycut(start_time,window,dasTMAR,'1'); [seisV_TBAG,hdr_TBAG] = segycut(start_time,window,dasTBAG,'1'); [seisV_TGRS,hdr_TGRS] = segycut(start_time,window,dasTGRS,'1'); [seisV_TBUG,hdr_TBUG] = segycut(start_time,window,dasTBUG,'1'); [seisV_TREX,hdr_TREX] = segycut(start_time,window,dasTREX,'1'); ts = (1:length(seisV_TBUG))/hdr_TBUG.sps; % create time series vector in seconds %% % Plotting seismic events all together clf subplot(7,1,1) plot(ts,seisV_TLAV,'k') %xlabel('time, hours') ylabel('counts') title('TLAV') % xlim([0,24]) % ylim([0,3000]) subplot(7,1,2) plot(ts,seisV_TLEF,'k') %xlabel('time, hours') ylabel('counts') title('TLEF') subplot(7,1,3) plot(ts,seisV_TMAR,'k') %xlabel('time, hours') ylabel('counts') title('MAR') subplot(7,1,4) plot(ts,seisV_TBAG,'k') %xlabel('time, hours') ylabel('counts') title('TBAG') subplot(7,1,5) plot(ts,seisV_TGRS,'k') %xlabel('time, hours') ylabel('counts') title('TGRS') % subplot(7,1,6) plot(ts,seisV_TBUG,'k') %xlabel('time, hours') ylabel('counts') title('TBUG') % % subplot(2,1,2) % plot(ts,seisV_TREX,'k') % %xlabel('time, hours') % ylabel('counts') % title('TREX') %% Signal Filtering sps = hdr_TLAV.sps; % get sample rate from header npoles = 2; % specify order of filter low = 1/4; % corner frequency at 1 Hz high = 4; % corner frequency at 3 Hz [B,A] = butter(npoles,[low high]/(sps/2)); tmp_seistraceTLAV = filter(B,A,seisV_TLAV); tmp_seistraceTLAV = tmp_seistraceTLAV/(max(tmp_seistraceTLAV)-min(tmp_seistraceTLAV)); tmp_seistraceTLAV = tmp_seistraceTLAV - mean(tmp_seistraceTLAV); tmp_seistraceTLEF = filter(B,A,seisV_TLEF); tmp_seistraceTLEF = tmp_seistraceTLEF/(max(tmp_seistraceTLEF)-min(tmp_seistraceTLEF)); tmp_seistraceTLEF = tmp_seistraceTLEF - mean(tmp_seistraceTLEF); tmp_seistraceTMAR = filter(B,A,seisV_TMAR); tmp_seistraceTMAR = tmp_seistraceTMAR/(max(tmp_seistraceTMAR)-min(tmp_seistraceTMAR)); tmp_seistraceTMAR = tmp_seistraceTMAR - mean(tmp_seistraceTMAR); tmp_seistraceTBAG = filter(B,A,seisV_TBAG); tmp_seistraceTBAG = tmp_seistraceTBAG/(max(tmp_seistraceTBAG)-min(tmp_seistraceTBAG)); tmp_seistraceTBAG = tmp_seistraceTBAG - mean(tmp_seistraceTBAG); tmp_seistraceTGRS = filter(B,A,seisV_TGRS); tmp_seistraceTGRS = tmp_seistraceTGRS/(max(tmp_seistraceTGRS)-min(tmp_seistraceTGRS)); tmp_seistraceTGRS = tmp_seistraceTGRS - mean(tmp_seistraceTGRS); tmp_seistraceTBUG = filter(B,A,seisV_TBUG); tmp_seistraceTBUG = tmp_seistraceTBUG/(max(tmp_seistraceTBUG)-min(tmp_seistraceTBUG)); tmp_seistraceTBUG = tmp_seistraceTBUG - mean(tmp_seistraceTBUG); tmp_seistraceTREX = filter(B,A,seisV_TREX); tmp_seistraceTREX = tmp_seistraceTREX/(max(tmp_seistraceTREX)-min(tmp_seistraceTREX)); tmp_seistraceTREX = tmp_seistraceTREX - mean(tmp_seistraceTREX); %% ts = (1:length(seisV))/hdr166.sps; % create time series vector in seconds ts172=(1:length(seisV172))/hdr166.sps; seis_sens = 800; % sensitivity is 800 V/m/s atod = hdr166.atod; % digitization V/count >> TRANSFORMA DE UN VALOR DIGITAL A UN ANALOGICO sens_gain = hdr166.gain; % programmed gain setting seis_calib = atod/seis_sens/sens_gain*1e6; % total conversion of counts to um/s >> CONVERSION DE CTAS A VALOR DE VELOCIDAD, MICROMETROS POR SEGUNDO % %% now plot the calibrated data trace for 24 hours % figure(1); clf % decimation = 10; % plot(ts(1:decimation:end)/60/60,seisV(1:decimation:end)*seis_calib,'k') % plot decimated data trace for 24 hours % ylabel('velocity (\mum/s)') % xlabel('time (hours)') % title(['das: ' das ' start time: ' num2str(start_time)]) %% Plot RSAM % Remove trace mean %norm166=seisV166-mean(seisV166); % figure(1) % now plot high pass filtered data traces for 1 minute window sps = hdr_TLEF.sps; % get sample rate from header npoles = 2; % specify order of filter low = 10; % corner frequency at 10 Hz high = (1/30); % corner frequency at 20 Hz %%%%% signal processing tool box with butter is necessary for following %%%%% operations % now do highpass % Sum individual time windows twin=10*60*100; tover=5*60*100; % Day 166 tmp_seistrace166 = filter(B,A,seisV166); tmp_seistrace166 = tmp_seistrace166/(max(tmp_seistrace166)-min(tmp_seistrace166)); tmp_seistrace166 = tmp_seistrace166 - mean(tmp_seistrace166); % Compute abs(waveform) tmp_seistrace166=abs(tmp_seistrace166); startsample166=1:tover:length(tmp_seistrace166)-twin; sumtrace166=zeros(length(startsample166)); for h=1:length(startsample166) sumtrace166(h)=(sum(tmp_seistrace166(startsample166(h):startsample166(h)+twin))); end % Day 167 tmp_seistrace167 = filter(B,A,seisV167); tmp_seistrace167 = tmp_seistrace167/(max(tmp_seistrace167)-min(tmp_seistrace167)); tmp_seistrace167 = tmp_seistrace167 - mean(tmp_seistrace167); tmp_seistrace167=abs(tmp_seistrace167); startsample167=1:tover:length(tmp_seistrace167)-twin; sumtrace167=zeros(length(startsample167)); for h=1:length(startsample167) sumtrace167(h)=(sum(tmp_seistrace167(startsample167(h):startsample167(h)+twin))); end % Day 168 tmp_seistrace168 = filter(B,A,seisV168); tmp_seistrace168 = tmp_seistrace168/(max(tmp_seistrace168)-min(tmp_seistrace168)); tmp_seistrace168 = tmp_seistrace168 - mean(tmp_seistrace168); tmp_seistrace168=abs(tmp_seistrace168); startsample168=1:tover:length(tmp_seistrace168)-twin; sumtrace168=zeros(length(startsample168)); for h=1:length(startsample168) sumtrace168(h)=(sum(tmp_seistrace168(startsample168(h):startsample168(h)+twin))); end % Day 169 tmp_seistrace169 = filter(B,A,seisV169); tmp_seistrace169 = tmp_seistrace169/(max(tmp_seistrace169)-min(tmp_seistrace169)); tmp_seistrace169 = tmp_seistrace169 - mean(tmp_seistrace169); tmp_seistrace169=abs(tmp_seistrace169); startsample169=1:tover:length(tmp_seistrace169)-twin; sumtrace169=zeros(length(startsample169)); for h=1:length(startsample169) sumtrace169(h)=(sum(tmp_seistrace169(startsample169(h):startsample169(h)+twin))); end % Day 170 tmp_seistrace170 = filter(B,A,seisV170); tmp_seistrace170 = tmp_seistrace170/(max(tmp_seistrace170)-min(tmp_seistrace170)); tmp_seistrace170 = tmp_seistrace170 - mean(tmp_seistrace170); tmp_seistrace170=abs(tmp_seistrace170); startsample170=1:tover:length(tmp_seistrace170)-twin; sumtrace170=zeros(length(startsample170)); for h=1:length(startsample170) sumtrace170(h)=(sum(tmp_seistrace170(startsample170(h):startsample170(h)+twin))); end % Day 171 tmp_seistrace171 = filter(B,A,seisV171); tmp_seistrace171 = tmp_seistrace171/(max(tmp_seistrace171)-min(tmp_seistrace171)); tmp_seistrace171 = tmp_seistrace171 - mean(tmp_seistrace171); tmp_seistrace171=abs(tmp_seistrace171); startsample171=1:tover:length(tmp_seistrace171)-twin; sumtrace171=zeros(length(startsample171)); for h=1:length(startsample171) sumtrace171(h)=(sum(tmp_seistrace171(startsample171(h):startsample171(h)+twin))); end % Day 172 tmp_seistrace172 = filter(B,A,seisV172); tmp_seistrace172 = tmp_seistrace172/(max(tmp_seistrace172)-min(tmp_seistrace172)); tmp_seistrace172 = tmp_seistrace172 - mean(tmp_seistrace172); tmp_seistrace172=abs(tmp_seistrace172); startsample172=1:tover:length(tmp_seistrace172)-twin; sumtrace172=zeros(length(startsample172)); for h=1:length(startsample172) sumtrace172(h)=(sum(tmp_seistrace172(startsample172(h):startsample172(h)+twin))); end % % Compute abs(waveform) % tmp_seistrace=abs(tmp_seistrace); % % % Sum individual time windows % twin=10*60*100; % tover=5*60*100; % % startsample=1:tover:length(tmp_seistrace)-twin; % sumtrace=zeros(length(startsample)); % % for h=1:length(startsample) % sumtrace(h)=(sum(tmp_seistrace(startsample(h):startsample(h)+twin))); % end %% Plot RSAM subplot(7,1,1) plot((startsample166-1)/100/60/60,sumtrace166,'k') %xlabel('time, hours') ylabel('counts') title('TLEF: Day 166') xlim([0,24]) % ylim([0,3000]) subplot(7,1,2) plot((startsample167-1)/100/60/60,sumtrace167,'k') %xlabel('time, hours') ylabel('counts') title('TLEF: Day 167') xlim([0,24]) % ylim([0,3000]) subplot(7,1,3) plot((startsample168-1)/100/60/60,sumtrace168,'k') %xlabel('time, hours') ylabel('counts') title('TLEF: Day 168') xlim([0,24]) % ylim([0,3000]) subplot(7,1,4) plot((startsample169-1)/100/60/60,sumtrace169,'k') %xlabel('time, hours') ylabel('counts') title('TLEF: Day 169') xlim([0,24]) % ylim([0,3000]) subplot(7,1,5) plot((startsample170-1)/100/60/60,sumtrace170,'k') %xlabel('time, hours') ylabel('counts') title('TLEF: Day 170') xlim([0,24]) % ylim([0,3000]) subplot(7,1,6) plot((startsample171-1)/100/60/60,sumtrace171,'k') %xlabel('time, hours') ylabel('counts') title('TLEF: Day 171') xlim([0,24]) % ylim([0,3000]) subplot(7,1,7) plot((startsample172-1)/100/60/60,sumtrace172,'k') xlim([0,24]) xlabel('time, hours') ylabel('counts') title('TLEF: Day 172') % ylim([0,3000]) %% now plot a traditional helicorder record of the same data figure(2); clf decimation = 2; scaling = 8; % scaling is set to the window, larger numbers amplify the traces by more ts_jday = hdr166.day + ts/60/60/24; % time axis for the ogram script must be decimal julian day ogram(ts_jday(1:decimation:end),seisV166(1:decimation:end),scaling) % this is J. Johnson's script figure(3); clf ts_jday = hdr167.day + ts/60/60/24; % time axis for the ogram script must be decimal julian day ogram(ts_jday(1:decimation:end),seisV167(1:decimation:end),scaling) % this is J. Johnson's script figure(4);clf ts_jday = hdr168.day + ts/60/60/24; % time axis for the ogram script must be decimal julian day ogram(ts_jday(1:decimation:end),seisV168(1:decimation:end),scaling) % this is J. Johnson's script figure(5);clf ts_jday = hdr169.day + ts/60/60/24; % time axis for the ogram script must be decimal julian day ogram(ts_jday(1:decimation:end),seisV169(1:decimation:end),scaling) % this is J. Johnson's script figure(6);clf ts_jday = hdr170.day + ts/60/60/24; % time axis for the ogram script must be decimal julian day ogram(ts_jday(1:decimation:end),seisV170(1:decimation:end),scaling) % this is J. Johnson's script figure(7);clf ts_jday = hdr171.day + ts/60/60/24; % time axis for the ogram script must be decimal julian day ogram(ts_jday(1:decimation:end),seisV171(1:decimation:end),scaling) % this is J. Johnson's script figure(8);clf ts_jday = hdr172.day + ts172/60/60/24; % time axis for the ogram script must be decimal julian day ogram(ts_jday(1:decimation:end),seisV172(1:decimation:end),scaling) % this is J. Johnson's script % ylabel('velocity (\mum/s)') % xlabel('time (hours)') % title(['das: ' das ' start time: ' num2str(start_time)]) %%