Commit 5f690223 authored by thomas.forbriger's avatar thomas.forbriger
Browse files

master|cross [MERGE]: fix issue #11 and use libdatrwxx

Merge branch 'cross'

- ts/wf/cross now provides proper scaling and a reasonable definition of
  sample times in both modes (fixes issue #11)
- ts/wf/cross now uses libdatrwxx and provides access to all supported input
  and output file formats
parents bfd8a465 db59f6a7
......@@ -195,7 +195,7 @@ evelox suseix: %x: %.o
# C++ programs
# ------------
sigval sigscale geophone gatherdiff cross autocorr: \
sigval sigscale geophone gatherdiff autocorr: \
%: %.o
$(LINK.cpp) $^ $(LOADLIBES) $(LDLIBS) -o $@ \
$(LIBDATRWXX) -ltfxx
......@@ -235,7 +235,7 @@ noisymize: \
$(LINK.cpp) $^ $(LOADLIBES) $(LDLIBS) -o $@ \
$(LIBTSIOXX)
fregra foutra deconv: \
cross fregra foutra deconv: \
%: %.o
$(LINK.cpp) $^ $(LOADLIBES) $(LDLIBS) -o $@ \
-lfourierxx $(LIBTSIOXX) -lfftw3
......
......@@ -4,11 +4,11 @@
* ----------------------------------------------------------------------------
*
* \author Thomas Forbriger
* \date 26/01/2005
* \date 16/11/2016
*
* cross correlation
*
* Copyright (c) 2005 by Thomas Forbriger (BFO Schiltach)
* Copyright (c) 2005, 2016 by Thomas Forbriger (BFO Schiltach)
*
* ----
* This program is free software; you can redistribute it and/or modify
......@@ -29,34 +29,44 @@
* REVISIONS and CHANGES
* - 26/01/2005 V1.0 Thomas Forbriger
* - 18/09/2015 V1.1 provide convolution as alternative operation
* - 16/11/2016 V1.2 provide new features:
* - full libdatrwxx support
* - trace selectors
*
* ============================================================================
*/
#define CROSS_VERSION \
"CROSS V1.1 cross correlation"
"CROSS V1.2 cross correlation and convolution"
#include <fstream>
#include <iostream>
#include <list>
#include <string>
#include <tfxx/commandline.h>
#include <tfxx/xcmdline.h>
#include <tfxx/error.h>
#include <sffxx.h>
#include <tsxx/tsxx.h>
#include <tfxx/rangelist.h>
#include <tfxx/rangestring.h>
#include <tsioxx/sfftimeseries.h>
#include <tsioxx/inputoperators.h>
#include <tsioxx/outputoperators.h>
#include <tsxx/correlate.h>
#include <tsxx/convolve.h>
#include <datrwxx/sff.h>
#include <datrwxx/writeany.h>
#include <datrwxx/readany.h>
#include <aff/seriesoperators.h>
using std::cout;
using std::cerr;
using std::endl;
typedef std::list<std::string> Tnames;
typedef ts::TDsfftimeseries Tts;
typedef aff::Series<double> Tseries;
typedef ts::sff::SFFTimeSeries<Tseries> Tts;
typedef Tts::Tdttimeseries Tdttimeseries;
struct Options {
bool verbose, overwrite, debug, convolve;
std::string inputformat, outputformat;
}; // struct Options
int main(int iargc, char* argv[])
......@@ -66,8 +76,10 @@ int main(int iargc, char* argv[])
char usage_text[]=
{
CROSS_VERSION "\n"
"usage: cross [-v] [-o] [-D] [-convolve] file [file [...]] outfile" "\n"
"usage: cross [-v] [-o] [-D] [-convolve] [--itype t] [--otype f]\n"
" file [f:F] [t:T] [file [f:F] [t:T] [...]] outfile\n"
" or: cross --help|-h" "\n"
" or: cross --xhelp" "\n"
};
// define full help text
......@@ -77,14 +89,51 @@ int main(int iargc, char* argv[])
"-D debug mode" "\n"
"-o overwrite existing output file" "\n"
"-convolve apply convolution rather than correlation\n"
"-itype t file format for input files\n"
"-otype t file format for output files\n"
"\n"
"file ... input file(s)" "\n"
" t:T select traces T, where T may be any range\n"
" specification like '3-4' or '5,6,7-12,20'\n"
" f:F specifies an input file format differing from\n"
" the format selected by \"--itype\"\n"
"outfile output file" "\n"
"\n"
"--xhelp print details on supported I/O formats\n"
"\n"
"The first trace read will be taken as a reference. All other\n"
"traces will be cross-correlated with the reference." "\n"
"traces will be cross-correlated with the reference or convolved\n"
"with the reference. Convolution and cross-correlation are understood\n"
"as an approximation to the corresponding integral form. Consequently\n"
"the sampling interval of time series must be identical and the\n"
"result of the discrete operation will be scaled with the sampling\n"
"interval." "\n"
"\n"
"Data time will be handled as follows:\n"
"\n"
"Convolution\n"
" The reference series is understood as the impulse response of a\n"
" filter which s to be applied to the other input series. If a\n"
" source date is given in the file of the reference trace, this will\n"
" be taken as the origin of time, allowing for acausal impulse\n"
" response. If it is missing, the first sample of the reference (filter\n"
" impulse response) is understood to be at lag=0s. The date of the\n"
" first sample of the output trace is set to account for possible\n"
" acausal parts of the filter response. Source definition is taken\n"
" only from the first of the input files and is passed unaltered to\n"
" the output.\n"
"\n"
"Cross-correlation\n"
" Each input trace is cross-correlated with the reference trace.\n"
" When calculating the lag time, absolute time of the first sample\n"
" of each of the traces is taken into account. Absolute time of the\n"
" output is set with respect to the source time. The source time is\n"
" taken from the file header of the reference trace. If a source time\n"
" is missing there, it will be set to a default value.\n"
"\n"
"If option -convolve is selected, all traces will be convolved with\n"
"the reference trace.\n"
"\n"
};
// define commandline options
......@@ -101,8 +150,27 @@ int main(int iargc, char* argv[])
{"D",arg_no,"-"},
// 4: convolve rather than correlate
{"convolve",arg_no,"-"},
// 5: detailed information on I/O formats
{"xhelp",arg_no,"-"},
// 6: convolve rather than correlate
{"itype",arg_yes,"sff"},
// 7: convolve rather than correlate
{"otype",arg_yes,"sff"},
{NULL}
};
/*----------------------------------------------------------------------*/
/* define file specific options
*/
static const char formatkey[]="f";
static const char tracekey[]="t";
// define commandline argument modifier keys
static const char* cmdlinekeys[]
={formatkey, tracekey, 0};
/*----------------------------------------------------------------------*/
// no arguments? print usage...
if (iargc<2)
......@@ -115,107 +183,224 @@ int main(int iargc, char* argv[])
Commandline cmdline(iargc, argv, options);
// help requested? print full help text...
if (cmdline.optset(0))
if (cmdline.optset(0) || cmdline.optset(5))
{
cerr << usage_text << endl;
cerr << help_text << endl;
cerr << endl;
datrw::supported_data_types(cerr);
cerr << endl;
if (cmdline.optset(5))
{
cerr << endl;
datrw::online_help(cerr);
}
exit(0);
}
/* dummy operation: print option settings
for (int iopt=0; iopt<2; iopt++)
{
cout << "option: '" << options[iopt].opt_string << "'" << endl;
if (cmdline.optset(iopt)) { cout << " option was set"; }
else { cout << "option was not set"; }
cout << endl;
cout << " argument (string): '" << cmdline.string_arg(iopt) << "'" << endl;
cout << " argument (int): '" << cmdline.int_arg(iopt) << "'" << endl;
cout << " argument (long): '" << cmdline.long_arg(iopt) << "'" << endl;
cout << " argument (float): '" << cmdline.float_arg(iopt) << "'" << endl;
cout << " argument (double): '" << cmdline.double_arg(iopt) << "'" << endl;
cout << " argument (bool): '";
if (cmdline.bool_arg(iopt))
{ cout << "true"; } else { cout << "false"; }
cout << "'" << endl;
}
while (cmdline.extra()) { cout << cmdline.next() << endl; }
// dummy operation: print rest of command line
while (cmdline.extra()) { cout << cmdline.next() << endl; }
*/
/* ---------------------------------------------------------------------- */
// evaluate command line parameters
Options opt;
opt.verbose=cmdline.optset(1);
opt.overwrite=cmdline.optset(2);
opt.debug=cmdline.optset(3);
opt.convolve=cmdline.optset(4);
opt.inputformat=cmdline.string_arg(6);
opt.outputformat=cmdline.string_arg(7);
if (opt.verbose) { cout << CROSS_VERSION << endl; }
TFXX_assert(cmdline.extra(), "missing input file!");
std::string outfile=cmdline.next();
TFXX_assert(cmdline.extra(), "missing output file!");
Tnames infiles;
while (cmdline.extra())
tfxx::cmdline::Tparsed arguments=parse_cmdline(cmdline, cmdlinekeys);
TFXX_assert(arguments.size()>1, "missing output file!");
// extract output file name and remove from list
std::string outfile=arguments.back().name;
arguments.pop_back();
if ((arguments.size()>1) && opt.verbose)
{
infiles.push_back(outfile);
outfile=cmdline.next();
cout
<< "NOTICE: file specific information (SRCE line and file FREE)" << endl
<< " will be taken from first file only!" << endl;
}
/* ---------------------------------------------------------------------- */
// open output file
// ----------------
if (opt.verbose) { cout << "open output file " << outfile << endl; }
std::ofstream ofs(outfile.c_str());
sff::SFFostream<Tts::Tseries> os(ofs);
// check if output file exists and open
if (!opt.overwrite) { datrw::abort_if_exists(outfile); }
std::ios_base::openmode oopenmode
=datrw::oanystream::openmode(opt.outputformat);
std::ofstream ofs(outfile.c_str(), oopenmode);
datrw::oanystream os(ofs, opt.outputformat, opt.debug);
sff::SRCE srce=sff::srce_reference();
os << srce;
Tts reference, data;
/*
* the first trace to be read has to be stored as the reference trace
*/
Tts reference;
bool referenceread=false;
Tnames::const_iterator I=infiles.begin();
while (I!=infiles.end())
// file header information for reference trace
ts::sff::FileHeader referencefileheader;
// cycle through list of file names found on command line
tfxx::cmdline::Tparsed::const_iterator infile=arguments.begin();
while (infile!=arguments.end())
{
if (opt.verbose) { cout << "open " << std::string(*I) << endl; }
std::ifstream ifs(I->c_str());
datrw::isffstream is(ifs, opt.debug);
if (opt.verbose) { cout << "open " << std::string(infile->name) << endl; }
std::string inputformat=opt.inputformat;
if (infile->haskey(formatkey))
{
inputformat=infile->value(formatkey);
}
std::ios_base::openmode iopenmode
=datrw::ianystream::openmode(inputformat);
std::ifstream ifs(infile->name.c_str(), iopenmode);
datrw::ianystream is(ifs, inputformat);
// handle file header
// ------------------
ts::sff::FileHeader fileheader;
is >> fileheader;
// cycle through traces of input file
// ----------------------------------
// setup trace selection
typedef tfxx::RangeList<int> Trangelist;
bool doselect=infile->haskey(tracekey);
Trangelist traceranges=
tfxx::string::rangelist<Trangelist::Tvalue>(infile->value(tracekey));
int iftrace=0;
while (is.good())
{
if (referenceread) {
if (opt.verbose) { cout << "read trace" << endl; }
Tts data;
Tts::Tseries samples;
is >> samples >> data.header;
data.series()=samples;
sff::INFO info;
is >> info;
if (opt.verbose) { cout << data.header.line() << endl; }
TFXX_assert((data.header.dt == reference.header.dt),
"sampling intervals do not match!");
Tts result;
if (opt.convolve)
{
result=ts::convolve(reference,data);
result.header.auxid="corr";
} else {
result=ts::correlate(reference,data);
result.header.auxid="conv";
}
// result *= result.header.dt;
result.header=data.header;
result.header.date=srce.date
+(result.first()*libtime::double2time(result.header.dt))
+(data.header.date-reference.header.date);
os << result.series() << result.header;
if (is.hasinfo()) { os << info; }
++iftrace;
if (doselect && (!traceranges.contains(iftrace)))
{
if (opt.verbose)
{ std::cout << " skip trace #" << iftrace << std::endl; }
is.skipseries();
}
else
{
if (opt.verbose) { cout << "read reference trace" << endl; }
is >> reference.series() >> reference.header;
if (opt.verbose) { cout << reference.header.line() << endl; }
referenceread=true;
}
if (opt.verbose) { std::cout << "process trace #" << iftrace << ": "; }
if (!referenceread)
{
// handle file header of reference
// -------------------------------
// set a SRCE definition if not present in input and mode is
// cross-correlation
sff::FREE freeblock(fileheader.free());
freeblock.append(CROSS_VERSION);
// SRCE date is a reference time for cross-correlograms
if (!opt.convolve)
{
if (!fileheader.hassrce())
{
sff::SRCE srce;
srce.date=libtime::TAbsoluteTime(2000);
fileheader.srce(srce);
freeblock.append("added SRCE definition to file header");
}
}
fileheader.free(freeblock);
os << fileheader;
referencefileheader=fileheader;
// read actual reference time series
// ---------------------------------
if (opt.verbose)
{ std::cout << "read reference trace" << std::endl; }
is >> reference;
referenceread=true;
if (opt.verbose) { cout << reference.header.wid2().line() << endl; }
// remember
}
else
{
if (opt.verbose) { cout << "read trace" << endl; }
Tts data;
is >> data;
if (opt.verbose) { cout << data.header.wid2().line() << endl; }
TFXX_assert((data.header.wid2().dt == reference.header.wid2().dt),
"sampling intervals do not match!");
Tts result(data);
sff::WID2 wid2line(data.header.wid2());
if (opt.convolve)
{
result=ts::convolve(reference,data);
wid2line.auxid="conv";
// set time of first sample
// ------------------------
if (referencefileheader.hasfree())
{
// offset of impulse response
libtime::TRelativeTime t0(reference.header.wid2().date
-referencefileheader.srce().date);
if (reference.header.wid2().date
<referencefileheader.srce().date)
{
wid2line.date -= t0;
}
else
{
wid2line.date += t0;
}
}
} else {
result=ts::correlate(reference,data);
wid2line.auxid="corr";
/*
* Time delay of current signal with respect to first sample of
* reference signal to be understood as a signed value:
*
* tor = signal.date-ref.date
*
* Duration of reference signal: Tr
* Duration of current signal: Ts
*
* Duration of cross-correlogram: Tr+Ts
*
* lag time of first sample:
*
* tl1=tor-(Tr+Ts)/2
*/
libtime::TRelativeTime
Tr(libtime::double2time(reference.header.wid2().dt
*(reference.header.wid2().nsamples-1)));
libtime::TRelativeTime
Ts(libtime::double2time(data.header.wid2().dt
*(data.header.wid2().nsamples-1)));
libtime::TRelativeTime tor(data.header.wid2().date
-reference.header.wid2().date);
libtime::TAbsoluteTime
torigin=referencefileheader.srce().date-((Tr+Ts)/2);
if (data.header.wid2().date
>=reference.header.wid2().date)
{
torigin += tor;
}
else
{
torigin -= tor;
}
wid2line.date=torigin;
}
result *= wid2line.dt;
result.header.wid2(wid2line);
/*
result.header.wid2().date=fileheader.srce.date
+(result.series.first()
*libtime::double2time(result.header.wid2().dt))
+(data.header.wid2().date-reference.header.wid2().date);
*/
if (opt.verbose) { cout << result.header.wid2().line() << endl; }
os << result.header << Tts::Tcoc(result);
} // else; if (!referenceread)
} // else; if (doselect && (!traceranges.contains(itrace)))
}
++I;
++infile;
}
}
......
......@@ -50,7 +50,7 @@ edit: flist; vim $<
clean: ;
-find . -name \*.bak | xargs --no-run-if-empty /bin/rm -v
-/bin/rm -vf flist *.xxx *.xxx.* *.sff *.TZ *.TR *.grn *.su *.bin
-/bin/rm -vf *.ps *.pdf *.asc *.ascii
-/bin/rm -vf *.ps *.pdf *.asc *.ascii
#----------------------------------------------------------------------
......@@ -134,6 +134,12 @@ convolution%.ps: Makefile.convolution
#----------------------------------------------------------------------
# test cross-correlation
cross%.ps: Makefile.crosscorrelation
$(MAKE) -f $< $@
#----------------------------------------------------------------------
%.psp: %.ps; gv $<; /bin/rm -fv $<
%.pdf: %.ps; ps2pdf $<
%.rw.ps: %.ps; gs -sDEVICE=pswrite -dNOPAUSE -dBATCH -sOutputFile=$@ $<
......
......@@ -33,8 +33,8 @@ convolution_%.sff: conv_%.sff conv_signal.sff
FILES=$(addsuffix .sff,$(addprefix convolution_delta,1 2 5 10 19))
DFILES=$(addsuffix .sff,$(addprefix conv_delta,1 2 5 10 19))
convolution_all_delta.ps: conv_signal.sff $(DFILES) $(FILES)
stuplox -s x -i -d $@/ps $^
stuplox -a -t -s x -i -d $@/ps $^
convolution_%.ps: conv_signal.sff conv_%.sff convolution_%.sff
stuplox -s x -i -Y 0.8 -d $@/ps $^
stuplox -a -t -c dta -s x -i -Y 0.8 -d $@/ps $^
# ----- END OF Makefile.convolution -----
# this is <Makefile.crosscorrelation>
# ----------------------------------------------------------------------------
#
# Copyright (c) 2016 by Thomas Forbriger (BFO Schiltach)
#
# test correlation operation
#
# REVISIONS and CHANGES
# 17/11/2016 V1.0 Thomas Forbriger
#
# ============================================================================
#
# provided targets:
#
# cross_files_box.ps cross_box.ps cross_boxrev.ps
# cross_files_sin3.ps cross_sin3.ps cross_sin3rev.ps
# cross_files_delta.ps cross_delta.ps cross_deltarev.ps
# cross_files_noise.ps cross_noise.ps cross_noiserev.ps
#
FILETYPE=ascii
DT=0.1
noise.raw.$(FILETYPE):
siggenx 14 $@ -v -o -ot $(FILETYPE) -a 1 -T 500. -d $(DT)
sin3.raw.$(FILETYPE):
siggenx 3 $@ -v -o -ot $(FILETYPE) -a 1 -T 50. -Ta 10. -Te 15. -d $(DT)
box.raw.$(FILETYPE):
siggenx 8 $@ -v -o -ot $(FILETYPE) -a 1 -T 50. -Ta 10. -n 50 -d $(DT)
delta.raw.$(FILETYPE):
siggenx 8 $@ -v -o -ot $(FILETYPE) -a 1 -T 50. -Ta 10. -n 1 -d $(DT)
STDDATE=2016/07/12 08:30:00
DATEDEL=80
DELDATE=$(shell date -d '$(STDDATE) $(DATEDEL) seconds' +'%Y/%m/%d %H:%M:%S')
%.std.$(FILETYPE): %.raw.$(FILETYPE)
sehefixx $< $@.xxx --verbose --overwrite \
--itype $(FILETYPE) --otype $(FILETYPE) -st "$(STDDATE)"
cooset -v -o --itype $(FILETYPE) --otype $(FILETYPE) -sot "$(STDDATE)" \
$@.xxx $@
%.alt.$(FILETYPE): %.raw.$(FILETYPE)
sehefixx $< $@.xxx --verbose --overwrite \
--itype $(FILETYPE) --otype $(FILETYPE) -st "$(DELDATE)"
cooset -v -o --itype $(FILETYPE) --otype $(FILETYPE) -sot "$(STDDATE)" \
$@.xxx $@
%.del.$(FILETYPE): %.$(FILETYPE)
printf "del 20.\nend\n" | tidofi -v -o -cs \
-type $(FILETYPE) -Type $(FILETYPE) $@ $<
%d.cross.$(FILETYPE): %.std.$(FILETYPE) %.std.del.$(FILETYPE)
cross -v -o --itype $(FILETYPE) --otype $(FILETYPE) \
$^ $@
%.cross.$(FILETYPE): %.std.$(FILETYPE) %.std.del.$(FILETYPE) \
%.alt.$(FILETYPE) %.alt.del.$(FILETYPE)
cross -v -o --itype $(FILETYPE) --otype $(FILETYPE) \
$^ $@
%rev.cross.$(FILETYPE): %.alt.del.$(FILETYPE) %.alt.$(FILETYPE) \
%.std.del.$(FILETYPE) %.std.$(FILETYPE)
cross -v -o --itype $(FILETYPE) --otype $(FILETYPE) \
$^ $@
cross_files_%.ps: %.std.$(FILETYPE) %.std.del.$(FILETYPE) \
%.alt.$(FILETYPE) %.alt.del.$(FILETYPE)
stuplox -g -c dta -ty $(FILETYPE) -a -t -g \
-s x -i -d $@/cps -Y 0.8 -N -C -l 1,2,4 $^
cross_%.ps: %.cross.$(FILETYPE)
stuplox -g -c dta -ty $(FILETYPE) -X 'lag / s' -st -g \
-s x -i -d $@/cps -Y 0.8 -N -C -l 1,2,4 $<
%.ps: %.$(FILETYPE)
stuplox -st -g -c dta -ty $(FILETYPE) -E -d $@/cps \
-Y 0.8 -N -C -l 1,2,4 $<
%.psp: %.ps; gv $<; /bin/rm -fv $<
# ----- END OF Makefile.crosscorrelation -----
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment