From 75e691de5f61323eb6965619fe494febd8b403cc Mon Sep 17 00:00:00 2001 From: Joanna Kocot Date: Fri, 5 Apr 2019 13:31:09 +0200 Subject: [PATCH] added TemplateMatchingDetection application --- src/StressMaps/LICENCE.txt | 2 +- src/TemplateMatchingDetection/LICENCE.txt | 15 + src/TemplateMatchingDetection/TM_EPOS.m | 83 +++++ .../compute_correlation.m | 68 +++++ src/TemplateMatchingDetection/find_events.m | 19 ++ src/TemplateMatchingDetection/jd2md.m | 38 +++ src/TemplateMatchingDetection/mk_template.m | 17 ++ src/TemplateMatchingDetection/read_xml.m | 141 +++++++++ src/TemplateMatchingDetection/readsac.m | 287 ++++++++++++++++++ 9 files changed, 669 insertions(+), 1 deletion(-) create mode 100644 src/TemplateMatchingDetection/LICENCE.txt create mode 100644 src/TemplateMatchingDetection/TM_EPOS.m create mode 100644 src/TemplateMatchingDetection/compute_correlation.m create mode 100644 src/TemplateMatchingDetection/find_events.m create mode 100644 src/TemplateMatchingDetection/jd2md.m create mode 100644 src/TemplateMatchingDetection/mk_template.m create mode 100644 src/TemplateMatchingDetection/read_xml.m create mode 100644 src/TemplateMatchingDetection/readsac.m diff --git a/src/StressMaps/LICENCE.txt b/src/StressMaps/LICENCE.txt index a35291a..e6eebd9 100644 --- a/src/StressMaps/LICENCE.txt +++ b/src/StressMaps/LICENCE.txt @@ -1,3 +1,3 @@ This application is a part of the IS-EPOS e-PLATFORM -Creative Commons Attribution-ShareAlike 4.0 International License: https://creativecommons.org/licenses/by-sa/4.0/ +This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License: https://creativecommons.org/licenses/by-sa/4.0/ diff --git a/src/TemplateMatchingDetection/LICENCE.txt b/src/TemplateMatchingDetection/LICENCE.txt new file mode 100644 index 0000000..f645b74 --- /dev/null +++ b/src/TemplateMatchingDetection/LICENCE.txt @@ -0,0 +1,15 @@ +This application is a part of the IS-EPOS e-PLATFORM + +Copyright 2019 Olivier Lengliné + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. \ No newline at end of file diff --git a/src/TemplateMatchingDetection/TM_EPOS.m b/src/TemplateMatchingDetection/TM_EPOS.m new file mode 100644 index 0000000..df8c553 --- /dev/null +++ b/src/TemplateMatchingDetection/TM_EPOS.m @@ -0,0 +1,83 @@ +function [T, yc] = TM_EPOS(pickinfo,event_wfm,continuousfile,nwin,npre,fmin,fmax, ... + min_dist,min_cor) + +% Check for integers +if(floor(nwin) ~= nwin); error('nwin is not an integer'); end +if(floor(npre) ~= npre); error('npre is not an integer'); end +if(floor(min_dist) ~= min_dist); error('min_dist is not an integer'); end +% Check frequencies +if(fmin>fmax); error('fmin > fmax'); end + + + +%% Parameters definition (example) +% % Filenames +% pickinfo = 'picks.xml'; % Picking files +% event_wfm= 'VN.TBVB..HHE.300000.SAC'; % Event waveform +% continuousfile= 'VN.TBVB..HHE.620000.SAC'; % Continuous waveform +% +% % Set Template window parameters +% nwin = 512; % Total number of samples in the template window +% npre = 50; % Number of points before the arrival time +% +% % Set filter parameters +% fmin = 1; % Minimum frequency +% fmax = 15; % Maximum frequency +% +% % Set parameters for declaring a detection +% min_dist = 200; % Minimum required distance (in samples) between 2 detections +% min_cor = 0.7; % Correlation threshold for declaring a detection (optional) + + +% Get the P-wave pick info +[tp,sta,ntwk,channel]=read_xml(pickinfo); + +% Read the continuous sac file +Fc=readsac(continuousfile); +if(Fc.tau==-1); error('The continuous sac file does not exist'); end +% And we will filter it. +fs = round(1./Fc.delta) ; fn = fs/2; +[b,a] = butter(4,[fmin fmax]./fn); +yc = Fc.trace - mean(Fc.trace); yc = filtfilt(b,a,yc); + + +% Read the event trace file +Fe=readsac(event_wfm); +if(Fe.tau==-1); error('The event sac file does not exist'); end + +% Check if the two waveforms correspond to the same station, channel, and +% have the same sampling frequency +if(Fe.delta ~= Fc.delta); error('Sampling frequencies are not the same'); end +if(Fe.kcmpnm ~= Fc.kcmpnm); error('Components are not the same'); end +if(Fe.kstnm ~= Fc.kstnm); error('Stations are not the same'); end + +% Check if picks correspond to the supplied waveform event file +if(strcmp(Fe.kstnm,sta)) + if(strcmp(Fe.kcmpnm,channel)) + if(strcmp(Fe.knetwk,ntwk)) + [mo,day]=jd2md(Fe.nzjday,Fe.nzyear); + t0 = datenum(Fe.nzyear,mo,day,Fe.nzhour,Fe.nzmin,Fe.sec); + t1 = t0 + (Fe.npts*Fe.delta)/(24*3600); + % Check if the picking time is within the event waveform + % file. It should be the case but just to be sure. + if( (tp > t0) && (tp < t1)) + % Extract the template waveform window + y=mk_template(Fe.trace,fmin,fmax,tp,t0,1./Fe.delta,nwin,npre); + % Compute the correlation + R=compute_correlation(yc,y); + % Find detections + T= find_events(R,min_dist,min_cor); + else + error('Picking time not within the event waveform duration') + end + else + error('Not the same Network') + end + else + error('Not the same Channel') + end +else + error('Not the same station') + +end + diff --git a/src/TemplateMatchingDetection/compute_correlation.m b/src/TemplateMatchingDetection/compute_correlation.m new file mode 100644 index 0000000..8ad1138 --- /dev/null +++ b/src/TemplateMatchingDetection/compute_correlation.m @@ -0,0 +1,68 @@ +function R=compute_correlation(v,H) + +% where H is the template waveform +% and v is the continuous signal +% we require that both signals are filtered similarly + +% Get fourier transform of the continuous signal +[e,nv] = size(v); +if(nv ==1); v = v'; nv = e; end % Make sure we have the good shape +f_v = fft(v); + + +% Transform template signal +[nh1,nh2] = size(H); % Get size +if(nh2 ==1); H = H'; nh2 = nh1; end % Make sure we have the good shape + +H = H - mean(H); % remove mean +H = H.*tukeywin(nh2,0.1)'; % Tapering +H = H./std(H); % Normalize by standard deviation + +npad = nv - nh2; +H = padarray(H,[0 npad],0,'post'); +f_H = fft(H); % Get the Fourier transform of H (flipped and same size +% as v. + +% Compute the standard deviation of continuous signal +s_v= std_v(v,nv,nh2); + +% Compute correlation +R = fliplr(real(ifft(f_H.*conj(f_v)))); + +% Get normalized values +R = R./s_v; + + + + +end + + + +function s_v= std_v(v,nv,nw) +% This function computes the stabdard deviation on the continuous signal +% on windows of size nw +% Compute std for the first window +s1 = sum(v(1:nw)); +s2 = v(1:nw)*v(1:nw)'; +s_v(1) = sqrt( (s2-s1*s1/nw)/(nw-1)); +% Initialize all values with this one +s_v(1:nv) = s_v(1); + +% Increment +for i = 2:nv-nw+1 + s1 = s1 - v(i-1) + v(i+nw-1); + s2 = s2 - v(i-1)*v(i-1) + v(i+nw-1)*v(i+nw-1); + s_v(i ) = sqrt( (s2-s1*s1/nw)/(nw-1)); +end + +% For the end of file +for i = nv-nw+2:nv + s1 = s1 - v(i-1); + s2 = s2 - v(i-1)*v(i-1); + s_v(i ) = sqrt( (s2-s1*s1/(nv-i))/(nv-i-1)); +end + +s_v = s_v * nw; + +end \ No newline at end of file diff --git a/src/TemplateMatchingDetection/find_events.m b/src/TemplateMatchingDetection/find_events.m new file mode 100644 index 0000000..a64eaa8 --- /dev/null +++ b/src/TemplateMatchingDetection/find_events.m @@ -0,0 +1,19 @@ +function T= find_events(R,min_dist,min_cor) + +% Take the enveloppe of the correlation signal +% (not very much used here but it s betetr in case we stack several +% correlation functions +R = abs(hilbert(R)); + +% Set the correrlation threshold based on the median absolute devaition if +% it has not been set +if nargin == 2 + min_cor = mean(R) + 8 * mad(R); +end + +% Build the time vector +t = 1:length(R); +% Find Peaks in the correlation function separated at least by min_dist and +% whose values are higher than min_cor +[~,T]=findpeaks(R,t,'MinpeakDistance',min_dist,'MinpeakHeight',min_cor); + diff --git a/src/TemplateMatchingDetection/jd2md.m b/src/TemplateMatchingDetection/jd2md.m new file mode 100644 index 0000000..c80ce37 --- /dev/null +++ b/src/TemplateMatchingDetection/jd2md.m @@ -0,0 +1,38 @@ +function [mo,day]=jd2md(jd,yr) +% function [mo,day]=jd2md(jd,) +% +% Function to transfer Julian day to date with the option to have leap years +% Input 'leap' or year number for if you want to use leap years +% Default is regular year +% +% Remark: a leap year is a year with 366 days: every 4th year is a leap year +% BUT every 100th is NOT _AND_ every 400th _IS_ again a leap year... +% +% György Hetényi +% 22 Nov 2004 +% ------------------------------------------------------------------ +for i=1:length(jd) + +if nargin<2 yr(i)=1; end +if ischar(yr(i))==1 & yr(i)=='leap' type='leap'; +elseif mod(yr(i),4)==0 & mod(yr(i),100)~=0 type='leap'; +elseif mod(yr(i),400)==0 type='leap'; +else type='regu'; +end + +clear mos +if type=='leap' +mos=[31,29,31,30,31,30,31,31,30,31,30,31]; +end +if type=='regu' +mos=[31,28,31,30,31,30,31,31,30,31,30,31]; +end + +k=1; +while jd(i)>sum(mos(1:k)) + k=k+1; +end +mo(i)=k; +day(i)=jd(i)-sum(mos(1:mo(i)-1)); + +end diff --git a/src/TemplateMatchingDetection/mk_template.m b/src/TemplateMatchingDetection/mk_template.m new file mode 100644 index 0000000..fc582d9 --- /dev/null +++ b/src/TemplateMatchingDetection/mk_template.m @@ -0,0 +1,17 @@ +function y=mk_template(y,fmin,fmax,tp,t0,fs,nwin,npre) + +% This function build the template waveform window + +% Filter the whole seismogram (to avoid edge problem) +fn = fs/2; +[b,a] = butter(4,[fmin fmax]./fn); +y = y - mean(y); +y = filtfilt(b,a,y); + +% Get the first point of the seismogram for windowing P-wave +I = round((tp-t0)*24*3600*fs); +I = I -npre; + +y = y(I:I+nwin-1); +y = y .* tukeywin(nwin,0.05); + diff --git a/src/TemplateMatchingDetection/read_xml.m b/src/TemplateMatchingDetection/read_xml.m new file mode 100644 index 0000000..087f00b --- /dev/null +++ b/src/TemplateMatchingDetection/read_xml.m @@ -0,0 +1,141 @@ +function [tp,sta,ntwk,channel]=read_xml(xmlfilename) + +% Load the xml File architecture +F = parseXML(xmlfilename); + +k1=xml_field(F,'eventParameters'); H = F.Children(k1); +k2=xml_field(H,'event'); H = H.Children(k2); +k3=xml_field(H,'pick'); +% Several picks are possible +for ik = 1:length(k3) + H_test = H.Children(k3(ik)); + % Test if we have the P wave + k=xml_field(H_test,'phaseHint'); + % If this is the P-wave read the informations + if(strcmp(H_test.Children(k).Children(1).Data,'P')) + k=xml_field(H_test,'waveformID'); A = H_test.Children(k); + [sta,ntwk,channel] = read_xml_attributes(A); + k=xml_field(H_test,'time'); G=H_test.Children(k); + k=xml_field(G,'value'); G = G.Children(k); + tp = datenum(G.Children(1).Data,'yyyy-mm-ddTHH:MM:SS.FFF'); + + end + +end + + + +end + +% Read waveform attributes +function [sta,ntwk,channel] = read_xml_attributes(A) +na = length(A.Attributes); +for ia = 1:na + if(strcmp(A.Attributes(ia).Name,'channelCode')) + channel = A.Attributes(ia).Value; + end + if(strcmp(A.Attributes(ia).Name,'networkCode')) + ntwk = A.Attributes(ia).Value; + end + if(strcmp(A.Attributes(ia).Name,'stationCode')) + sta = A.Attributes(ia).Value; + end +end + +end + +% Get the current xml field with the name keystring +function [ikey]=xml_field(F,keystring) + +nchilds = length(F.Children); +ikey = []; +for i = 1:nchilds + if(strcmp(F.Children(i).Name,keystring)) + ikey = [ikey i]; + end +end +nmatch = length(ikey); +if(nmatch == 0) + error('No field %s', keystring) +end + +end + +function theStruct = parseXML(filename) +% PARSEXML Convert XML file to a MATLAB structure. +% THis functions is from MATLAB xmlread doc +try + tree = xmlread(filename); +catch + error('Failed to read XML file %s.',filename); +end + +% Recurse over child nodes. This could run into problems +% with very deeply nested trees. +try + theStruct = parseChildNodes(tree); +catch + error('Unable to parse XML file %s.',filename); +end + +end + +% ----- Local function PARSECHILDNODES ----- +function children = parseChildNodes(theNode) +% Recurse over node children. +children = []; +if theNode.hasChildNodes + childNodes = theNode.getChildNodes; + numChildNodes = childNodes.getLength; + allocCell = cell(1, numChildNodes); + + children = struct( ... + 'Name', allocCell, 'Attributes', allocCell, ... + 'Data', allocCell, 'Children', allocCell); + + for count = 1:numChildNodes + theChild = childNodes.item(count-1); + children(count) = makeStructFromNode(theChild); + end +end + +end + +% ----- Local function MAKESTRUCTFROMNODE ----- +function nodeStruct = makeStructFromNode(theNode) +% Create structure of node info. + +nodeStruct = struct( ... + 'Name', char(theNode.getNodeName), ... + 'Attributes', parseAttributes(theNode), ... + 'Data', '', ... + 'Children', parseChildNodes(theNode)); + +if any(strcmp(methods(theNode), 'getData')) + nodeStruct.Data = char(theNode.getData); +else + nodeStruct.Data = ''; +end + +end + +% ----- Local function PARSEATTRIBUTES ----- +function attributes = parseAttributes(theNode) +% Create attributes structure. + +attributes = []; +if theNode.hasAttributes + theAttributes = theNode.getAttributes; + numAttributes = theAttributes.getLength; + allocCell = cell(1, numAttributes); + attributes = struct('Name', allocCell, 'Value', ... + allocCell); + + for count = 1:numAttributes + attrib = theAttributes.item(count-1); + attributes(count).Name = char(attrib.getName); + attributes(count).Value = char(attrib.getValue); + end +end + +end \ No newline at end of file diff --git a/src/TemplateMatchingDetection/readsac.m b/src/TemplateMatchingDetection/readsac.m new file mode 100644 index 0000000..7b7c0ed --- /dev/null +++ b/src/TemplateMatchingDetection/readsac.m @@ -0,0 +1,287 @@ +function F=readsac(files) +% F=readsac('files') +% +% Read a SAC-file or a set of SAC-files. +% +% Input : +% 'files' corresponds either to the complete name of one SAC-file, +% either to a set of SAC-files with the logical expression *. +% Output : +% F is a structure (length equal to the number of SAC-files read). +% It contains ALL the fields of the SAC's header and the trace in +% the field "trace". +% +% EXAMPLE : +% F=readsac('2004.02.23-17.31.00.ABH.SHZ.SAC') +% reads the file 2004.02.23-17.31.00.ABH.SHZ.SAC and saves it into +% the structure F (dimension: 1*1). +% F=readsac('Data/*BHZ*.SAC') +% reads all the files of the form *BHZ*.SAC in the directory "Data". +% The size of F equals the number of these files. +% +% From J. Vergne and G. Hetenyi + + +%------------------------------------------------------------------ +% By default, the signals are read in mode 'r' +modlect='r'; + +% Determine the type of "files" +rep_files=fileparts(files); +list_files=dir(files); + +if length(list_files)==0 +% disp('"readsac": File(s) do not exist'); + F.tau=-1; +else + +for ifile=1:length(list_files) + + nomfich=fullfile(rep_files,list_files(ifile).name); + clear dsac; + + % Read signals + % ------------ + % Following the signal type, reading procedure can be different: + % 'l' - IEEE floating point with little-endian byte ordering + % 'b' - IEEE floating point with big-endian byte ordering + + [C,MAXSIZE,ENDIAN]=computer; + if ENDIAN=='L' bool_l=0; + else bool_l=1; + end + + fidev=fopen(nomfich,modlect,'l'); +if fidev > 0 + + h1=fread(fidev,70,'float'); % -------------- + h2=fread(fidev,40,'long'); % reading header + h3=fread(fidev,192,'uchar'); % -------------- + + % testing byte-order, h2(7) must! be 6. + if h2(7)~=6 + bool_l=1; + fclose(fidev); + fidev=fopen(nomfich,modlect,'b'); + h1=fread(fidev,70,'float'); + h2=fread(fidev,40,'long'); + h3=fread(fidev,192,'uchar'); + end + + % reading trace + tamp=fread(fidev,inf,'float'); + + dsac.h1=h1; + dsac.h2=h2; + dsac.h3=h3; + + + % PART 1: reading float-type parameters + % ------------------------------------------ + % %- comment are from original version, + % % comments ares from SAC-homepage + + dsac.delta=h1(1); %- sampling time interval + dsac.depmin=h1(2); % minimum value of dependent variable + dsac.depmax=h1(3); % maximum value of dependent variable + dsac.scale=h1(4); % multiplying scale factor for dependent variable (not currently used) + dsac.odelta=h1(5); % observed increment if different from nominal value + + dsac.b=h1(6); %- begin time (d�calage du 1er �chantillon) (beginning value of independent variable) + dsac.e=h1(7); % ending value of independent variable + dsac.o=h1(8); %- event origin time (seconds relative to reference time) + dsac.a=h1(9); % first arrival time (seconds relative to reference time) + dsac.internal10=h1(10); % INTERNAL + + dsac.t0=h1(11); %- user defined time picks or markers + dsac.t1=h1(12); %- (seconds relative to reference time) + dsac.t2=h1(13); %- + dsac.t3=h1(14); %- + dsac.t4=h1(15); %- + dsac.t5=h1(16); % + dsac.t6=h1(17); % + dsac.t7=h1(18); % + dsac.t8=h1(19); % + dsac.t9=h1(20); % + + dsac.f=h1(21); % fini or end of event time (seconds relative to reference time) + dsac.resp0=h1(22); % instrument response parameters (not currently used) + dsac.resp1=h1(23); % + dsac.resp2=h1(24); % + dsac.resp3=h1(25); % + dsac.resp4=h1(26); % + dsac.resp5=h1(27); % + dsac.resp6=h1(28); % + dsac.resp7=h1(29); % + dsac.resp8=h1(30); % + dsac.resp9=h1(31); % + + + dsac.stla=h1(32); %- station latitude (degrees, north positive) + dsac.stlo=h1(33); %- station longitude (degrees, east positive) + dsac.stel=h1(34); %- station elevation (meters) + dsac.stdp=h1(35); % station depth below surface (meters)(not currently used) + + dsac.evla=h1(36); %- event latitude (degrees, north positive) + dsac.evlo=h1(37); %- event longitude (degrees, east positive) + dsac.evel=h1(38); % event elevation (meters)(not currently used) + dsac.evdp=h1(39); %- event depth below surface (meters) + dsac.mag=h1(40); % event magnitude + + + dsac.user0=h1(41); %- used defined variable storage area, floating! + dsac.user1=h1(42); %- + dsac.user2=h1(43); %- + dsac.user3=h1(44); % + dsac.user4=h1(45); % + dsac.user5=h1(46); % + dsac.user6=h1(47); % + dsac.user7=h1(48); % + dsac.user8=h1(49); % + dsac.user9=h1(50); % + + dsac.dist=h1(51); %- station to event distance (km) + dsac.az=h1(52); %- event to station azimuth (degrees) + dsac.baz=h1(53); %- station to event azimuth (degrees) + dsac.gcarc=h1(54); %- station to event great circle arc length (degrees) + dsac.internal55=h1(55); % INTERNAL + + dsac.internal56=h1(56); % INTERNAL + dsac.depmen=h1(57); % mean value of dependent variable + dsac.cmpaz=h1(58); %- component azimuth (degrees clockwise from north) + dsac.cmpinc=h1(59); %- component incidence angle (degrees from vertical) + + dsac.xminimum=h1(60); % minimum value of X (spectral files only) + dsac.xmaximum=h1(61); % maximum value of X (spectral files only) + dsac.yminimum=h1(62); % minimum value of Y (spectral files only) + dsac.ymaximum=h1(63); % maximum value of Y (spectral files only) + dsac.unused64=h1(64); % UNUSED + dsac.unused65=h1(65); % UNUSED + dsac.unused66=h1(66); % UNUSED + dsac.unused67=h1(67); % UNUSED + dsac.unused68=h1(68); % UNUSED + dsac.unused69=h1(69); % UNUSED + dsac.unused70=h1(70); % UNUSED + + + % PART 2: reading long-type parameters + % ------------------------------------ + + % GMT time corresponding to reference (0) time in file + dsac.nzyear=h2(1); %- year + dsac.nzjday=h2(2); %- julian day + dsac.nzhour=h2(3); %- hour + dsac.nzmin=h2(4); %- minute + dsac.nzsec=h2(5); % second + dsac.nzmsec=h2(6); % millisecond + dsac.sec=h2(5)+h2(6)/1000; %- full second (NOT an original SAC-attribute!) + + dsac.nvhdr=h2(7); % header version number: 6! + dsac.norid=h2(8); % origin ID (CSS 3.0) + dsac.nevid=h2(9); % event ID (CSS 3.0) + + dsac.npts=h2(10); %- number of points per data component + + dsac.internal81=h2(11);% INTERNAL + dsac.nwfid=h2(12); % waveform ID (CSS 3.0) + dsac.nxsize=h2(13); % spectral length (spectral files only) + dsac.nysize=h2(14); % spectral width (spectral files only) + dsac.unused85=h2(15); % UNUSED + + dsac.iftype=h2(16); % type of file (required)(see SAC-page) + dsac.idep=h2(17); % type of dependent variable (see SAC-page) + dsac.iztype=h2(18); % reference time equivalence (see SAC-page) + dsac.unused89=h2(19); % UNUSED + dsac.iinst=h2(20); % type of recording instrument (not currently used) + + dsac.istreg=h2(21); % station geogrpahic region (not currently used) + dsac.ievreg=h2(22); % event geographic location (not currently used) + dsac.ievtyp=h2(23); % type of event (see SAC-page) + dsac.iqual=h2(24); % quality of data (not currently used)(see SAC-page) + dsac.isynth=h2(25); % synthetic data flag (not currently used) + + dsac.imagtyp=h2(26); % magnitude type + dsac.imagsrc=h2(27); % source of magnitude information + dsac.unused98=h2(28); % UNUSED + dsac.unused99=h2(29); % UNUSED + dsac.unused100=h2(30); % UNUSED + + dsac.unused101=h2(31); % UNUSED + dsac.unused102=h2(32); % UNUSED + dsac.unused103=h2(33); % UNUSED + dsac.unused104=h2(34); % UNUSED + dsac.unused105=h2(35); % UNUSED + + dsac.leven=h2(36); % TRUE if data is evenly spaced + dsac.lspol=h2(37); % TRUE if station components have positive polarity (left-hand rule) + dsac.lovrok=h2(38); % TRUE if it is okay to overwrite this file on disk + dsac.lcalda=h2(39); % TRUE if DIST,AZ,BAZ and GCARC are to be calculated form station and event coordinates + dsac.unused110=h2(40); % UNUSED + + + % PART 3: reading uchar-type parameters + % ------------------------------------- + + imin1=min(find(h3(1:8)==0 | h3(1:8)==32)); + imin2=min(find(h3(9:24)==0 | h3(9:24)==32)); + + dsac.kstnm=rm_blanc(h3(1:1+imin1-1))'; %- station name + dsac.kevnm=rm_blanc(h3(9:9+imin2-1))'; %- event name + + dsac.khole=rm_blanc(h3(25:32))'; % hole identification if nuclear event + dsac.ko=rm_blanc(h3(33:40))'; % event origin time identification + dsac.ka=rm_blanc(h3(41:48))'; % first arrival time identification + + dsac.kt0=rm_blanc(h3(49:56))'; %- user defined time pick identifications, h1(11:20) + dsac.kt1=rm_blanc(h3(57:64))'; %- + dsac.kt2=rm_blanc(h3(65:72))'; %- + dsac.kt3=rm_blanc(h3(73:80))'; %- + dsac.kt4=rm_blanc(h3(81:88))'; %- + dsac.kt5=rm_blanc(h3(89:96))'; % + dsac.kt6=rm_blanc(h3(97:104))'; % + dsac.kt7=rm_blanc(h3(105:112))'; % + dsac.kt8=rm_blanc(h3(113:120))'; % + dsac.kt9=rm_blanc(h3(121:128))'; % + + dsac.kf=rm_blanc(h3(129:136))'; % fini identification + + dsac.kuser0=rm_blanc(h3(137:144))'; %- user defined variable storage area + dsac.kuser1=rm_blanc(h3(145:152))'; %- + dsac.kuser2=rm_blanc(h3(153:160))'; %- + + dsac.kcmpnm=rm_blanc(h3(161:168))'; %- component name + dsac.knetwk=rm_blanc(h3(169:176))'; % name of seismic network + dsac.kdatrd=rm_blanc(h3(177:184))'; % date data was read onto computer + + dsac.kinst=rm_blanc(h3(185:192))'; %- generic name of recording instrument + + + % reading trace + % ------------- + + dsac.trace=tamp; + + dsac.tau = 1; % Added to check if the file exist + fclose(fidev); + +else + %disp(['"readsac": file ' nomfich ' : reading error - file does not exist']) + dsac.tau=-1; +end + F(ifile)=dsac; +end + + +end + +% ------------------------------------------------------------- + +function ch1=rm_blanc(ch) + +% looks for whitespace in 'ch' and removes them + if ischar(ch) + ch1=ch(find(double(ch)~=32) & double(ch)~=0); + else + ch1=ch(find(ch~=32 & ch~=0)); + ch1=char(ch1); + end