croposp.cc 12.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
/*! \file croposp.cc
 * \brief Cross power spectral density
 * 
 * ----------------------------------------------------------------------------
 * 
 * \author Thomas Forbriger
 * \date 25/12/2018
 * 
 * Cross power spectral density
 * 
 * Copyright (c) 2018 by Thomas Forbriger (BFO Schiltach) 
12
13
14
15
16
17
 *
 * ----
 * 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. 
18
 * 
19
20
21
22
23
24
25
26
27
 * 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/>.
 * ----
 *
28
29
30
31
32
33
34
35
36
 * REVISIONS and CHANGES 
 *  - 25/12/2018   V1.0   Thomas Forbriger
 * 
 * ============================================================================
 */
#define CROPOSP_VERSION \
  "CROPOSP   V1.0   Cross power spectral density"

#include <iostream>
37
38
39
#include <vector>
#include <string>
#include <sstream>
40
#include <tfxx/commandline.h>
41
#include <tfxx/error.h>
42
#include <tfxx/misc.h>
43
#include <tfxx/stringfunc.h>
44
#include <tsioxx/cmdlinefiles.h>
45
#include <tsxx/tscollection.h>
46
47
48
49
50

using std::cout;
using std::cerr;
using std::endl;

51
struct Options {
52
  bool verbose, trim, debug;
53
  std::string inputformat, labelpattern;
54
  double datetolerance;
55
  bool logscale;
56
57
58
  unsigned int n_per_decade;
  bool compute_psd, compute_npsd, compute_coherency, compute_transfer;
  std::string outfile_psd, outfile_npsd, outfile_coherency, outfile_transfer;
59
60
}; // struct Options

61
62
63
64
65
66
67
// 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
68
typedef ts::TimeSeriesCollection<double> TCollection;
69
70
// a vector to gold meta data for each time series
typedef std::vector<Meta> TMetaVector;
71
72

/* ====================================================================== */
73
74
75
76
77
78
79
80
81
82
83
84
  
// define usage information
char reference_sleeman_et_al_2006[]=
{
  "Reinoud Sleeman, Arie van Wettum, and Jeannot Trampert, 2006.\n"
  "  Three-Channel Correlation Analysis: A New Technique to Measure\n"
  "  Instrumental Noise of Digitizers and Seismic Sensors.\n"
  "  Bull. Seism. Soc. Am., Vol. 96, No. 1, pp. 258–271.\n"
  "  doi: 10.1785/0120050032\n"
  "  http://www.geo.uu.nl/~seismain/pdf/bssa06-inst.pdf (accessed 2019-01-25)\n" };

/* ====================================================================== */
85

86
87
88
89
90
91
92
93
94
95
96
97
98
99
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)

/* ---------------------------------------------------------------------- */

100
void report_collection(const TCollection& collection,
101
                       const TMetaVector& metadata,
102
                       const bool& debug=false)
103
{
104
105
  TFXX_assert(collection.size() == metadata.size(),
              "data inconsitency; report this as a program bug");
106
  TCollection::const_iterator i_series=collection.begin();
107
108
109
  TMetaVector::const_iterator i_meta=metadata.begin();
  while (i_series != collection.end() &&
         i_meta != metadata.end())
110
111
112
113
114
115
116
117
118
  {
    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;
119
    cout << "  " << i_meta->label << endl;
120
121
122
123
124
125
    TFXX_debug(debug, "report_collection",
               TFXX_value(i_series->f())
               << " " <<
               TFXX_value(i_series->l())
               << " " <<
               TFXX_value(i_series->size()));
126
    ++i_series;
127
    ++i_meta;
128
129
130
131
132
  }
} // void report_collection(const Tcollection& collection)

/* ====================================================================== */

133
134
135
136
137
138
139
int main(int iargc, char* argv[])
{

  // define usage information
  char usage_text[]=
  {
    CROPOSP_VERSION "\n"
140
141
    "usage: croposp [-verbose] [-itype f] [-trim] [-datetolerance t]" "\n"
    "               [-label p]" "\n"
142
    "               [-psd f] [-npsd f] [-transfer f] [-coherency f]" "\n"
143
    "       file [t:sel] [f:format] [n:label] [file [t:s] [f:f] [n:l]] [...]\n"
144
145
146
147
148
149
    "   or: croposp --help|-h" "\n"
  };

  // define full help text
  char help_text[]=
  {
150
    "-verbose           be verbose\n"
151
    "-DEBUG             produce debugging output\n"
152
153
154
155
156
    "-itype f           set default input file format (this value can\n"
    "                   be overwritten by the format key at a file name)\n"
    "-trim              trim time series to identical begin and length\n"
    "-datetolerance t   set tolerance as a fraction of sampling interval\n"
    "                   when comparing time series header date values\n"
157
158
159
160
161
162
163
164
    "-label p           set pattern for time series label\n"
    "                   %C: channel (WID2 header field)\n"
    "                   %S: station (WID2 header field)\n"
    "                   %I: instype (WID2 header field)\n"
    "                   %A: auxid (WID2 header field)\n"
    "                   %N: label set for file name\n"
    "                   %F: file name\n"
    "                   %NT: number of trace in file\n"
165
166
    "-log n             map averages to logarithmic scale with \"n\"\n"
    "                   samples pre decade\n"
167
    "\n"
168
169
170
171
172
173
174
175
176
177
178
    "output options:\n"
    "-psd f        compute power spectral density and write to file \"f\"\n"
    "-npsd f       compute power spectral density of incoherent component\n"
    "              and write to file \"f\"\n"
    "              Sleeman et al (2006, eq. 14)\n"
    "-transfer f   compute amplitude transfer function for pairs of\n"
    "              signals and write to file \"f\"\n"
    "              Sleeman et al (2006, eq. 13)\n"
    "-coherency f  compute coherecy for pairs of signals\n"
    "              and write to file \"f\"\n"
    "\n"
179
    "keys to be appended to file names if required:\n"
180
181
    "t:sel        trace selection\n"
    "f:format     file format\n"
182
    "n:label      label to be used in plot legend\n"
183
184
185
186
187
188
189
190
191
  };

  // define commandline options
  using namespace tfxx::cmdline;
  static Declare options[]= 
  {
    // 0: print help
    {"help",arg_no,"-"},
    // 1: verbose mode
192
193
194
    {"verbose",arg_no,"-"},
    // 2: default input file format
    {"itype",arg_yes,"sff"},
195
196
197
198
    // 3: trim time series
    {"trim",arg_no,"-"},
    // 4: default input file format
    {"datetolerance",arg_yes,"0."},
199
200
    // 5: debug mode
    {"DEBUG",arg_no,"-"},
201
202
    // 6: pattern for trace label
    {"pattern",arg_yes,"%S:%C:%A:%I"},
203
204
205
206
207
208
209
210
    // 7: compute PSD
    {"psd",arg_yes,"-"},
    // 8: compute PSD of incoherent component
    {"npsd",arg_yes,"-"},
    // 9: compute transfer function
    {"transfer",arg_yes,"-"},
    // 10: compute coherency
    {"coherency",arg_yes,"-"},
211
212
    // 11: compute coherency
    {"log",arg_yes,"1"},
213
214
215
216
217
    {NULL}
  };

  // define command line keys for input files
  static const char tracekey[]="t";
218
  static const char formatkey[]="f";
219
  static const char labelkey[]="n";
220
221

  // define commandline argument modifier keys
222
  static const char* cmdlinekeys[]={tracekey, formatkey, labelkey, 0};
223

224
225
226
227
228
  /* ====================================================================== */
  /* analyze command line
   * --------------------
   */

229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
  // no arguments? print usage...
  if (iargc<2) 
  {
    cerr << usage_text << endl;
    exit(0);
  }

  // collect options from commandline
  Commandline cmdline(iargc, argv, options);

  // help requested? print full help text...
  if (cmdline.optset(0))
  {
    cerr << usage_text << endl;
    cerr << help_text << endl;
244
245
246
247
    cerr << endl;
    cerr << "References" << endl
      <<    "----------" << endl;
    cerr << reference_sleeman_et_al_2006 << endl;
248
249
    exit(0);
  }
250
251
252
253

  Options opt;
  opt.verbose=cmdline.optset(1);
  opt.inputformat=cmdline.string_arg(2);
254
255
  opt.trim=cmdline.optset(3);
  opt.datetolerance=cmdline.double_arg(4);
256
  opt.debug=cmdline.optset(5);
257
  opt.labelpattern=cmdline.string_arg(6);
258

259
260
261
262
263
264
265
266
267
  opt.compute_psd=cmdline.optset(7);
  opt.outfile_psd=cmdline.string_arg(7);
  opt.compute_npsd=cmdline.optset(8);
  opt.outfile_npsd=cmdline.string_arg(8);
  opt.compute_transfer=cmdline.optset(9);
  opt.outfile_transfer=cmdline.string_arg(9);
  opt.compute_coherency=cmdline.optset(10);
  opt.outfile_coherency=cmdline.string_arg(10);

268
269
270
271
272
273
  opt.logscale=cmdline.optset(11);
  opt.n_per_decade=cmdline.int_arg(11);

  TFXX_assert(opt.n_per_decade>0,
              "number of samples per decade must be finite and positive");

274
275
  TFXX_assert(((opt.datetolerance >= 0.) && (opt.datetolerance <=1.)),
              "datetolerance is out of accepted range");
276
277
278
279
280
  
  // extract commandline arguments
  TFXX_assert(cmdline.extra(), "missing input file");
  tfxx::cmdline::Tparsed arguments=parse_cmdline(cmdline, cmdlinekeys);

281
282
283
284
285
  /* ====================================================================== */
  /* read data, prepare labels and trim time series if requested
   * -----------------------------------------------------------
   */

286
287
288
289
290
  // read data files
  if (opt.verbose) 
  {
    cout << "Read time series data" << endl;
  }
291
292
293
294
295
296
297
298
299
300
301
302
303
304
  ts::sff::TDFileList input_file_list;
  tfxx::cmdline::Tparsed::const_iterator file=arguments.begin();
  while (file != arguments.end())
  {
    std::string format=opt.inputformat;
    if (file->haskey(formatkey))
    {
      format=file->value(formatkey);
    }
    input_file_list.push_back(ts::sff::readDSFF(*file, opt.verbose, 
                                                tracekey, format));
    ++file;
  }

305
  // collect traces
306
  TCollection collection_of_series;
307
  TMetaVector vector_of_metadata;
308
309
310
311
312
313
314
315
  typedef ts::sff::DFile::Tfile Tfile;
  typedef Tfile::Ttracevector Ttracevector;
  ts::sff::TDFileList::const_iterator i_file=input_file_list.begin();
  while (i_file != input_file_list.end())
  {
    Ttracevector::const_iterator i_trace=i_file->data.begin();
    while (i_trace != i_file->data.end())
    {
316
      collection_of_series.push_back(*i_trace);
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
      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);
      std::ostringstream oss;
      oss << i_trace->traceindex();
      metadata.label=patsubst(metadata.label, "%NT", oss.str());
      if (metadata.filearguments.haskey(labelkey))
      {
        metadata.label=patsubst(metadata.label, "%N", 
                                metadata.filearguments.value(labelkey));
      }
      vector_of_metadata.push_back(metadata);
339
340
341
342
343
344
345
346
347
      ++i_trace;
    }
    ++i_file;
  }

  // report traces
  if (opt.verbose)
  {
    cout << "Time series to be analyzed:" << endl;
348
349
    report_collection(collection_of_series,
                      vector_of_metadata);
350
351
  }

352
353
354
  // trim time series
  if (opt.trim)
  {
355
356
    TFXX_assert(collection_of_series.overlap(),
                "time series do not overlap");
357
358
    collection_of_series.trim_to_date();
    collection_of_series.trim_to_nsamples();
359
360
361
362
363

    // report traces
    if (opt.verbose)
    {
      cout << "Time series after trimming:" << endl;
364
365
      report_collection(collection_of_series,
                        vector_of_metadata, opt.debug);
366
    }
367
368
369
370
371
372
373
374
  }

  // check consistency
  ::sff::WID2compare comparer(::sff::Fdt | ::sff::Fnsamples | ::sff::Fdate);
  comparer.setdatetolerance(opt.datetolerance);
  TFXX_assert(collection_of_series.are_consistent(comparer),
              "time series are inconsistent");

375
376
  /* ====================================================================== */

377
378
379
380
381
  TFXX_assert(!opt.compute_psd, "-psd is not yet implemented");
  TFXX_assert(!opt.compute_npsd, "-npsd is not yet implemented");
  TFXX_assert(!opt.compute_transfer, "-transfer is not yet implemented");
  TFXX_assert(!opt.compute_coherency, "-coherency is not yet implemented");

382
383
384
} // main()

/* ----- END OF croposp.cc ----- */