Commit e35d1807 authored by thomas.forbriger's avatar thomas.forbriger
Browse files

Merge branch 'croposp' of gitserver into croposp

parents b7b40282 0cc6979f
......@@ -23,6 +23,7 @@ ts/wf
ts/hd
ts/refract
ts/cal
ts/croposp
ts/noise
ts/fidase
ts/lisousi
......
this is <COPYING>
============================================================================
croposp: cross power spectral density
============================================================================
Copyright (C) 2019 by Thomas Forbriger
----
These programs are free software; you can redistribute them and/or modify
them under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
These programs are distributed in the hope that they will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
----
----- END OF COPYING -----
# this is <Makefile.croposp>
# this is <Makefile>
# ----------------------------------------------------------------------------
#
# Copyright (c) 2019 by Thomas Forbriger (BFO Schiltach)
#
# provide test cases for proposp
# compile croposp
#
# ----
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# ----
#
# REVISIONS and CHANGES
# 01/01/2019 V1.0 Thomas Forbriger
# 06/02/2019 V1.0 Thomas Forbriger
#
# ============================================================================
#
all: cpsd1.bin cpsd2.bin cpsd3.bin
PROGRAMS=croposp
SCRIPTS=croposplot.py
INSTALLEDPROGRAMS=$(addprefix $(LOCBINDIR)/,$(PROGRAMS))
INSTALLEDSCRIPTS=$(addprefix $(LOCSCRIPTDIR)/,$(SCRIPTS))
.PHONY: all
all: install doc
.PHONY: list-programs
list-programs:
echo $(PROGRAMS) | tr ' ' '\n' | sort
.PHONY: doc
doc: doxydoc
.PHONY: install
install: $(INSTALLEDPROGRAMS) $(INSTALLEDSCRIPTS)
$(LOCBINDIR)/%: %
mkdir -pv $(LOCBINDIR)
/bin/mv -fv $< $(LOCBINDIR)
$(LOCSCRIPTDIR)/%: %
mkdir -pv $(LOCSCRIPTDIR)
/bin/cp -fv $< $(LOCSCRIPTDIR)
CHECKVAR=$(if $($(1)),,$(error ERROR: missing variable $(1)))
CHECKVARS=$(foreach var,$(1),$(call CHECKVAR,$(var)))
$(call CHECKVARS,LOCINCLUDEDIR LOCLIBDIR LOCBINDIR)
$(call CHECKVARS,TF_LINK_FORTRAN TF_LINK_PGPLOT)
# use STATIC=-static to produce statically linked binaries
STATIC=
FLAGS += $(MYFLAGS)
CPPFLAGS=$(addprefix -I,$(LOCINCLUDEDIR) $(subst :, ,$(SERVERINCLUDEDIR))) \
$(FLAGS)
CXXFLAGS+=-Wall $(FLAGS)
LDFLAGS=$(addprefix -L,$(LOCLIBDIR) $(subst :, ,$(SERVERLIBDIR))) $(STATIC)
#----------------------------------------------------------------------
# standard edit targets
.PHONY: clean
clean: ;
-find . -name \*.bak | xargs --no-run-if-empty /bin/rm -v
-find . -name \*.d | xargs --no-run-if-empty /bin/rm -v
-/bin/rm -vf flist *.o *.xxx *.ps *~ $(PROGRAMS)
TESTCASEMAKE=$(filter-out %.bak,$(wildcard testcases/Makefile*))
# the frame of doxygen documentation is placed in text files
DOXYTXT=$(wildcard doxygen*.txt */doxygen*txt)
EDITFILES=Makefile README $(wildcard *.cfg) COPYING $(DOXYTXT) \
$(filter-out %.bak,$(wildcard README* Makefile*))
EDITSRC=$(wildcard *.cc *.h *.py)
flist: $(EDITFILES) $(TF_EDIT) $(EDITSRC)
echo $(filter $(EDITFILES),$^) | tr ' ' '\n' | sort > $@
echo "---- source code ----" >> $@
echo $(filter $(EDITSRC),$^) | tr ' ' '\n' | sort >> $@
echo "----" >> $@
echo $(filter-out $(EDITFILES) $(EDITSRC),$^) \
| tr ' ' '\n' | sort >> $@
.PHONY: edit
edit: flist; vim $<
#----------------------------------------------------------------------
# pattern rules
#----------------------------------------------------------------------
# description and online texts
%.cc %.h: %_help.txt
echo "// DO NOT EDIT: this file is automatically derived from $<" \
> $(patsubst %.txt,%.h,$<)
echo "extern char $(patsubst %.txt,%,$<)[];" >> $(patsubst %.txt,%.h,$<)
echo "// DO NOT EDIT: this file is automatically derived from $<" \
> $(patsubst %.txt,%.cc,$<)
echo "#include \"$(patsubst %.txt,%.h,$<)\"" >> $(patsubst %.txt,%.cc,$<)
echo "char $(patsubst %.txt,%,$<)[]=" >> $(patsubst %.txt,%.cc,$<)
echo "{" >> $(patsubst %.txt,%.cc,$<)
cat $< | egrep -v '^#' | sed -e 's/"/\\"/g' \
| sed -e 's/$$/\\n"/' | sed -e 's/^/ "/'\
>> $(patsubst %.txt,%.cc,$<)
echo "};" >> $(patsubst %.txt,%.cc,$<)
#----------------------------------------------------------------------
# dependencies
%.d: %.cc
$(SHELL) -ec '$(CXX) -M $(CPPFLAGS) $< \
| sed '\''s,\($(notdir $*)\)\.o[ :]*,$(dir $@)\1.o $@ : ,g'\'' \
> $@; \
[ -s $@ ] || rm -f $@'
include $(patsubst %.txt,%.d,$(wildcard *_help.txt))
include $(patsubst %.cc,%.d,$(wildcard *.cc))
#----------------------------------------------------------------------
# typical combinations of libraries
LIBDATRWXX=-ldatrwxx -lsffxx -ltime++ -lgsexx -laff
# libgsl is required for random sequence generator in libtsxx
LIBTSIOXX=-ltsioxx -ltsxx -ltfxx $(LIBDATRWXX) -lgsl -lgslcblas
#----------------------------------------------------------------------
# binary executable targets
# -------------------------
OBJFILES=pairs_adapter.o \
pairs.o \
patsubst.o \
report_collection.o \
triples.o \
write_name_series.o
croposp: \
%: %.o $(OBJFILES)
$(LINK.cpp) $^ $(LOADLIBES) $(LDLIBS) -o $@ \
-lpsdxx -lfourierxx $(LIBTSIOXX) -lfftw3
#======================================================================
# test case
# ---------
alldata: cpsd1.bin cpsd2.bin cpsd3.bin
CPSD_NSAMPLESRAW=100000
#CPSD_NSAMPLESRAW=1000
......@@ -27,7 +166,7 @@ CPSD_AMP3=3.e-3
CPSD_AMPBG=4.e-5
# location of plot tool
CROPOSPLOT=../../../python/visu/croposplot.py
CROPOSPLOT=./croposplot.py
cpsd1raw.bin:
siggenx 12 $@ $(CPSD_SIGGENOPT) -a $(CPSD_AMP1)
......@@ -103,5 +242,52 @@ psd.pdf: psd.xxx
$(CROPOSPLOT) -v --grid -o $@ $<
%.pdp: %.pdf; evince $<; /bin/rm -fv $<
#
#======================================================================
# documentation part
# ------------------
#
# targets commonly used:
# ----------------------
#
# make doxyclean removes all documentation
# make doxydoc creates doxygen documentation in the DOXYWWWPATH
# make doxyview creates doxygen documentation and launches netscape to
# browse in the documentation
# make doxyconf edit the doxygen configuration file
#
# If you launch "make doxydoc" the documentation will be written to
# DOXYWWWPATH (see below). This is meant to export the documentation through
# your homepage. The doxyfull directory is just a symbolic link to this
# directory.
#
$(call CHECKVARS,TF_WWWBASEDIR TF_BROWSER)
DOXYWWWPATH=$(TF_WWWBASEDIR)/croposp
.PHONY: doxyclean doxyview doxydoc doxyconf
doxyclean: ;/bin/rm -rfv $(DOXYWWWPATH)
DOXYSRC=$(DOXYTXT) $(wildcard *.h *.cc *.f)
# create doxygen intermediate configuration
PWD=$(shell env pwd)
doxydoc.xxx: doxydoc.cfg
sed 's,<OUTPUTDIRECTORY>,$(DOXYWWWPATH),g;s,<STRIPFROMPATH>,$(PWD),g' \
$< > $@
# create commented version of doxygen configuration
doxycomm.xxx: doxydoc.cfg
/bin/cp -vf $< $@; doxygen -u $@
$(DOXYWWWPATH)/html/index.html: doxydoc.xxx $(DOXYSRC)
mkdir -vp $(DOXYWWWPATH)
doxygen $<
doxydoc: $(DOXYWWWPATH)/html/index.html
doxyview: $(DOXYWWWPATH)/html/index.html
$(TF_BROWSER) file:$< &
# ----- END OF Makefile.croposp -----
# ----- END OF Makefile -----
this is <README>
============================================================================
croposp: cross power spectral density
============================================================================
croposp provides:
- computation of power spectral density from time series data
- computation of coherency from pairs of time series data
- computation of power spectral density of incoherent signal component,
amplitude and phase transfer functions from triples of time series data
(Sleeman, 2006)
croposplot.py produces diagrams from output files produced by croposp.
References
----------
Reinoud Sleeman, Arie van Wettum, and Jeannot Trampert, 2006.
Three-Channel Correlation Analysis: A New Technique to Measure
Instrumental Noise of Digitizers and Seismic Sensors.
Bull. Seism. Soc. Am., Vol. 96, No. 1, pp. 258–271.
doi: 10.1785/0120050032
http://www.geo.uu.nl/~seismain/pdf/bssa06-inst.pdf (accessed 2019-01-25)
============================================================================
Installation
------------
This software is part of the project Seitosh. See README.1st in the root
directory of the collection or https://git.scc.kit.edu/Seitosh/Seitosh for
general installation instructions.
The Makefile supports program compilation and linking. Some binary libraries
are required.
The command
make all
will compile and install the binary executables as well as the doxygen
documentation.
Environment variables control where the results are stored and where
libraries and library header files are expected:
LOCLIBDIR defines location of binary libraries
LOCINCLUDEDIR defines location of C/C++ header files (prototypes)
LOCBINDIR defines location of binary executables
TF_WWWBASEDIR defines location of doxygen output
Dependencies:
Compilers required to build the programs:
C++ compiler
C/C++ preprocessor
doxygen (required to process source code documentation)
extern libraries/packages needed to compile the code
libfftw3: Fast Fourier Transformation
The header files for these libraries are required as well. Under OpenSuSE
you have to install the respective devel packages.
Seitosh libraries required to compile the code:
libaff
libdatrwxx
libfourier (libfourier, libfourierxx)
libgsexx
libsffxx
libtfxx
libtime (libtime, libtime++)
libtsxx
libpsdxx
============================================================================
The home of this software is
https://git.scc.kit.edu/Seitosh/Seitosh/tree/master/src/ts/croposp
Please send bug reports and suggestions to
Thomas.Forbriger@kit.edu
----- END OF README -----
......@@ -30,21 +30,17 @@
*
* ============================================================================
*/
/*
* version string is set in croposp.h
*
#define CROPOSP_VERSION \
"CROPOSP V1.0 Cross power spectral density"
*/
#include "croposp.h"
#include <iostream>
#include <fstream>
#include <vector>
#include <string>
#include <sstream>
#include <tfxx/commandline.h>
#include <tfxx/error.h>
#include <tfxx/misc.h>
#include <tfxx/stringfunc.h>
#include <tsioxx/cmdlinefiles.h>
#include <tsxx/tscollection.h>
#include <psdxx/psd.h>
#include <tfxx/misc.h>
using std::cout;
using std::cerr;
......@@ -66,48 +62,6 @@ struct Options {
std::string prefix_phase, prefix_coherency;
}; // struct Options
// a struct to hold meta data
struct Meta {
tfxx::cmdline::Filename filearguments;
std::string label;
}; // struct Meta
// type of time series collection to hold input data
typedef ts::TimeSeriesCollection<double> TCollection;
// a vector to gold meta data for each time series
typedef std::vector<Meta> TMetaVector;
// a struct to hold computation results for exactly one series
struct NamedSeries {
std::string label;
psd::TDseries series;
}; // struct NamedSeries
// a vector type to hold results
typedef std::vector<NamedSeries> TNamedSeriesVector;
/* ====================================================================== */
/*! provide index values to pairs
*
* cycle through a complete set
*/
class Pairs {
public:
Pairs(const unsigned int& n);
}; // class Pairs
/* ---------------------------------------------------------------------- */
/*! provide index values to triples
*
* cycle through a complete set
*/
class Triples {
public:
Triples(const unsigned int& n);
}; // class Triples
/* ====================================================================== */
// define usage information
......@@ -122,107 +76,6 @@ char reference_sleeman_et_al_2006[]=
/* ====================================================================== */
std::string patsubst(const std::string& templatestring,
const std::string& pattern,
const std::string& content)
{
std::string retval;
retval=tfxx::string::patsubst(templatestring, pattern,
tfxx::string::trimws(content));
return(retval);
} // std::string patsubst(const std::string& templatestring,
// const std::string& pattern,
// const std::string& content)
/* ---------------------------------------------------------------------- */
void report_collection(const TCollection& collection,
const TMetaVector& metadata,
const bool& debug=false)
{
TFXX_assert(collection.size() == metadata.size(),
"data inconsitency; report this as a program bug");
TCollection::const_iterator i_series=collection.begin();
TMetaVector::const_iterator i_meta=metadata.begin();
while (i_series != collection.end() &&
i_meta != metadata.end())
{
cout << " "
<< i_series->header.station << " "
<< i_series->header.channel << " "
<< i_series->header.auxid << " "
<< "dt=" << i_series->header.dt << "s "
<< "begin: " << i_series->header.date.timestring() << " "
<< "n: " << i_series->header.nsamples
<< endl;
cout << " " << i_meta->label << endl;
TFXX_debug(debug, "report_collection",
TFXX_value(i_series->f())
<< " " <<
TFXX_value(i_series->l())
<< " " <<
TFXX_value(i_series->size()));
++i_series;
++i_meta;
}
} // void report_collection(const Tcollection& collection)
/* ---------------------------------------------------------------------- */
/* output a vector of named series to file
* ---------------------------------------
*/
void write_named_series(const std::string& filename,
const std::string& comment,
psd::TDseries f,
const TNamedSeriesVector& nsv,
const bool& verbose,
const bool& overwrite=true)
{
if (verbose)
{
cout << endl
<< "output to file " << filename << ":" << endl
<< comment << endl;
}
if (!overwrite)
{
std::ifstream file(filename.c_str(), std::ios_base::in);
TFXX_assert((!file.good()),"ERROR: output file exists!");
}
TFXX_assert(nsv.size()>0, "data container is empty");
std::ofstream os(filename.c_str());
os << "# " << CROPOSP_VERSION << endl;
os << "# " << comment << endl;
for (unsigned int i=0; i<nsv.size(); ++i)
{
unsigned int fi=i+1;
os << "# #" << fi << ": " << nsv[i].label << endl;
TFXX_assert(nsv[i].series.size()==f.size(),
"series passed to output function are inconsistent")
}
// set first index to 0 - just in case
f.shift(-f.first());
unsigned int nsamples=nsv[0].series.size();
for (unsigned int i=0; i<nsamples; ++i)
{
os << f(i);
for (unsigned int j=0; j<nsv.size(); ++j)
{
const psd::TDseries::Tcoc& s=nsv[j].series;
os << " " << s(i+s.first());
}
os << endl;
}
} // void write_named_series(const std::string& filename,
// const std::string& comment,
// psd::TDseries f,
// const TNamedSeriesVector& nsv,
// const bool& verbose,
// const bool& overwrite=true)
/* ====================================================================== */
int main(int iargc, char* argv[])
{
......@@ -459,8 +312,8 @@ int main(int iargc, char* argv[])
}
// collect traces
TCollection collection_of_series;
TMetaVector vector_of_metadata;
croposp::TCollection collection_of_series;
croposp::TMetaVector vector_of_metadata;
typedef ts::sff::DFile::Tfile Tfile;
typedef Tfile::Ttracevector Ttracevector;
ts::sff::TDFileList::const_iterator i_file=input_file_list.begin();
......@@ -470,26 +323,26 @@ int main(int iargc, char* argv[])
while (i_trace != i_file->data.end())
{
collection_of_series.push_back(*i_trace);
Meta metadata;
croposp::Meta metadata;
metadata.filearguments=i_file->arguments;
metadata.label=opt.labelpattern;
metadata.label=patsubst(metadata.label, "%C",
i_trace->header.wid2().channel);
metadata.label=patsubst(metadata.label, "%S",
i_trace->header.wid2().station);
metadata.label=patsubst(metadata.label, "%A",
i_trace->header.wid2().auxid);
metadata.label=patsubst(metadata.label, "%I",
i_trace->header.wid2().instype);
metadata.label=patsubst(metadata.label, "%F",
metadata.filearguments.name);
metadata.label=croposp::patsubst(metadata.label, "%C",
i_trace->header.wid2().channel);
metadata.label=croposp::patsubst(metadata.label, "%S",
i_trace->header.wid2().station);
metadata.label=croposp::patsubst(metadata.label, "%A",
i_trace->header.wid2().auxid);
metadata.label=croposp::patsubst(metadata.label, "%I",
i_trace->header.wid2().instype);
metadata.label=croposp::patsubst(metadata.label, "%F",
metadata.filearguments.name);
std::ostringstream oss;
oss << i_trace->traceindex();
metadata.label=patsubst(metadata.label, "%NT", oss.str());
metadata.label=croposp::patsubst(metadata.label, "%NT", oss.str());
if (metadata.filearguments.haskey(labelkey))
{
metadata.label=patsubst(metadata.label, "%N",
metadata.filearguments.value(labelkey));
metadata.label=croposp::patsubst(metadata.label, "%N",
metadata.filearguments.value(labelkey));
}
vector_of_metadata.push_back(metadata);
++i_trace;
......@@ -545,7 +398,7 @@ int main(int iargc, char* argv[])
*/
// provide container for PSD values
TNamedSeriesVector PSD_vector;
croposp::TNamedSeriesVector PSD_vector;
if (opt.compute_psd || opt.compute_npsd ||
opt.compute_coherency)
......@@ -558,8 +411,8 @@ int main(int iargc, char* argv[])
psd_computer.set_padfactor(opt.padfactor);
psd_computer.set_overlap(opt.overlap);
TCollection::const_iterator i_series=collection_of_series.begin();
TMetaVector::const_iterator i_meta=vector_of_metadata.begin();
croposp::TCollection::const_iterator i_series=collection_of_series.begin();
croposp::TMetaVector::const_iterator i_meta=vector_of_metadata.begin();
while (i_series != collection_of_series.end() &&
i_meta != vector_of_metadata.end())
{
......@@ -585,7 +438,7 @@ int main(int iargc, char* argv[])
} // if (opt.logscale) ... else
} // if (i_series==collection_of_series.begin())
NamedSeries named_psd;
croposp::NamedSeries named_psd;
if (opt.logscale)
{
named_psd.series=psd::log_sampling(psd, frequencies);
......@@ -614,6 +467,11 @@ int main(int iargc, char* argv[])
} // if (opt.compute_psd || opt.compute_npsd || opt.compute_transfer ||
// // opt.compute_coherency)
/* ---------------------------------------------------------------------- */
// provide container for CPSD values
croposp::TNamedCPSDVector CPSD_vector;
} // main()
/* ----- END OF croposp.cc ----- */
/*! \file croposp.h
* \brief declaration of classes, types and functions used in croposp
*
* ----------------------------------------------------------------------------
*
* \author Thomas Forbriger
* \date 06/02/2019
*
* declaration of classes, types and functions used in croposp (prototypes)
*
* Copyright (c) 2019 by Thomas Forbriger (BFO Schiltach)
*
* ----
* This program is free