croposp.cc 17.5 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
#include <fstream>
38
39
40
#include <vector>
#include <string>
#include <sstream>
41
#include <tfxx/commandline.h>
42
#include <tfxx/error.h>
43
#include <tfxx/misc.h>
44
#include <tfxx/stringfunc.h>
45
#include <tsioxx/cmdlinefiles.h>
46
#include <tsxx/tscollection.h>
47
#include <psdxx/psd.h>
48
49
50
51
52

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

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

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

74
75
76
77
78
79
80
81
82
// 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;

83
/* ====================================================================== */
84
85
86
87
88
89
90
91
92
93
94
95
  
// 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" };

/* ====================================================================== */
96

97
98
99
100
101
102
103
104
105
106
107
108
109
110
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)

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

111
void report_collection(const TCollection& collection,
112
                       const TMetaVector& metadata,
113
                       const bool& debug=false)
114
{
115
116
  TFXX_assert(collection.size() == metadata.size(),
              "data inconsitency; report this as a program bug");
117
  TCollection::const_iterator i_series=collection.begin();
118
119
120
  TMetaVector::const_iterator i_meta=metadata.begin();
  while (i_series != collection.end() &&
         i_meta != metadata.end())
121
122
123
124
125
126
127
128
129
  {
    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;
130
    cout << "  " << i_meta->label << endl;
131
132
133
134
135
136
    TFXX_debug(debug, "report_collection",
               TFXX_value(i_series->f())
               << " " <<
               TFXX_value(i_series->l())
               << " " <<
               TFXX_value(i_series->size()));
137
    ++i_series;
138
    ++i_meta;
139
140
141
  }
} // void report_collection(const Tcollection& collection)

142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
/* ---------------------------------------------------------------------- */
/* 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)

196
197
/* ====================================================================== */

198
199
200
201
202
203
204
int main(int iargc, char* argv[])
{

  // define usage information
  char usage_text[]=
  {
    CROPOSP_VERSION "\n"
205
    "usage: croposp [-verbose] [-itype f] [-trim] [-datetolerance t]" "\n"
206
    "               [-label p] [-overwrite] [-DEBUG]" "\n"
207
    "               [-psd f] [-npsd f] [-transfer f] [-coherency f]" "\n"
208
    "       file [t:sel] [f:format] [n:label] [file [t:s] [f:f] [n:l]] [...]\n"
209
210
211
212
213
214
    "   or: croposp --help|-h" "\n"
  };

  // define full help text
  char help_text[]=
  {
215
    "-verbose           be verbose\n"
216
    "-DEBUG             produce debugging output\n"
217
218
219
220
221
    "-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"
222
223
224
225
226
227
228
229
    "-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"
230
231
    "-log n             map averages to logarithmic scale with \"n\"\n"
    "                   samples pre decade\n"
232
    "-overwrite         overwrite existing output files\n"
233
    "\n"
234
235
236
237
238
239
240
241
242
243
244
    "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"
245
    "keys to be appended to file names if required:\n"
246
247
    "t:sel        trace selection\n"
    "f:format     file format\n"
248
    "n:label      label to be used in plot legend\n"
249
250
251
252
253
254
255
256
257
  };

  // define commandline options
  using namespace tfxx::cmdline;
  static Declare options[]= 
  {
    // 0: print help
    {"help",arg_no,"-"},
    // 1: verbose mode
258
259
260
    {"verbose",arg_no,"-"},
    // 2: default input file format
    {"itype",arg_yes,"sff"},
261
262
263
264
    // 3: trim time series
    {"trim",arg_no,"-"},
    // 4: default input file format
    {"datetolerance",arg_yes,"0."},
265
266
    // 5: debug mode
    {"DEBUG",arg_no,"-"},
267
268
    // 6: pattern for trace label
    {"pattern",arg_yes,"%S:%C:%A:%I"},
269
270
271
272
273
274
275
276
    // 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,"-"},
277
    // 11: map to logarithmic frequency sampling
278
    {"log",arg_yes,"1"},
279
280
    // 12: overwrite existing output files
    {"overwrite",arg_no,"-"},
281
282
283
284
285
    {NULL}
  };

  // define command line keys for input files
  static const char tracekey[]="t";
286
  static const char formatkey[]="f";
287
  static const char labelkey[]="n";
288
289

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

292
293
294
295
296
  /* ====================================================================== */
  /* analyze command line
   * --------------------
   */

297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
  // 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;
312
313
314
315
    cerr << endl;
    cerr << "References" << endl
      <<    "----------" << endl;
    cerr << reference_sleeman_et_al_2006 << endl;
316
317
    exit(0);
  }
318
319
320
321

  Options opt;
  opt.verbose=cmdline.optset(1);
  opt.inputformat=cmdline.string_arg(2);
322
323
  opt.trim=cmdline.optset(3);
  opt.datetolerance=cmdline.double_arg(4);
324
  opt.debug=cmdline.optset(5);
325
  opt.labelpattern=cmdline.string_arg(6);
326

327
328
329
330
331
332
333
334
335
  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);

336
337
  opt.logscale=cmdline.optset(11);
  opt.n_per_decade=cmdline.int_arg(11);
338
  opt.overwrite=cmdline.optset(12);
339
340
341
342

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

343
344
  TFXX_assert(((opt.datetolerance >= 0.) && (opt.datetolerance <=1.)),
              "datetolerance is out of accepted range");
345
346
347
348
349
  
  // extract commandline arguments
  TFXX_assert(cmdline.extra(), "missing input file");
  tfxx::cmdline::Tparsed arguments=parse_cmdline(cmdline, cmdlinekeys);

350
351
352
353
354
  /* ====================================================================== */
  /* read data, prepare labels and trim time series if requested
   * -----------------------------------------------------------
   */

355
356
357
358
359
  // read data files
  if (opt.verbose) 
  {
    cout << "Read time series data" << endl;
  }
360
361
362
363
364
365
366
367
368
369
370
371
372
373
  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;
  }

374
  // collect traces
375
  TCollection collection_of_series;
376
  TMetaVector vector_of_metadata;
377
378
379
380
381
382
383
384
  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())
    {
385
      collection_of_series.push_back(*i_trace);
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
      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);
408
409
410
411
412
413
414
415
416
      ++i_trace;
    }
    ++i_file;
  }

  // report traces
  if (opt.verbose)
  {
    cout << "Time series to be analyzed:" << endl;
417
418
    report_collection(collection_of_series,
                      vector_of_metadata);
419
420
  }

421
422
423
  // trim time series
  if (opt.trim)
  {
424
425
    TFXX_assert(collection_of_series.overlap(),
                "time series do not overlap");
426
427
    collection_of_series.trim_to_date();
    collection_of_series.trim_to_nsamples();
428
429
430
431
432

    // report traces
    if (opt.verbose)
    {
      cout << "Time series after trimming:" << endl;
433
434
      report_collection(collection_of_series,
                        vector_of_metadata, opt.debug);
435
    }
436
437
438
439
440
441
442
443
  }

  // 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");

444
445
  /* ====================================================================== */

446
447
448
449
  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");

450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
  // provide container for frequency values (will be filled upon computation
  // of PSD values)
  psd::TDseries frequencies;

  /* ---------------------------------------------------------------------- */
  /* compute power spectral density
   * ------------------------------
   */

  // provide container for PSD values
  TNamedSeriesVector PSD_vector;

  if (opt.compute_psd || opt.compute_npsd || opt.compute_transfer ||
      opt.compute_coherency)
  {
    psd::DPSDComputer psd_computer;
466
467
    psd_computer.set_debug(opt.debug);
    psd_computer.set_verbose(opt.verbose);
468
469
470
471
472
473
474
475
    TCollection::const_iterator i_series=collection_of_series.begin();
    TMetaVector::const_iterator i_meta=vector_of_metadata.begin();
    while (i_series != collection_of_series.end() &&
           i_meta != vector_of_metadata.end())
    {
      psd::TDISeries interval_series;
      interval_series.data=*i_series;
      interval_series.interval=i_series->header.dt;
476
477
478
479
      TFXX_debug(opt.debug, "croposp (compute PSD)",
                 TFXX_value(i_series->header.dt)
                 << " "
                 << TFXX_value(interval_series.interval));
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
      psd::TDISeries psd=psd_computer.psd(interval_series);

      if (i_series==collection_of_series.begin())
      {
        if (opt.logscale)
        {
          frequencies=psd::log_frequency(psd.interval, psd.data.size(),
                                         opt.n_per_decade);
        }
        else
        {
          frequencies=psd::lin_frequency(psd.interval, psd.data.size());
        } // if (opt.logscale) ... else
      } // if (i_series==collection_of_series.begin())

      NamedSeries named_psd;
      if (opt.logscale)
      {
        named_psd.series=psd::log_sampling(psd, frequencies);
      } // if (opt.logscale)
      else
      {
        named_psd.series=psd.data;
      } // if (opt.logscale) ... else

      named_psd.label=i_meta->label;
      PSD_vector.push_back(named_psd);

      ++i_series;
      ++i_meta;
    } // while (i_series != collection_of_series.end() &&
      //        i_meta != vector_of_metadata.end())
512
513
514
515
516
517
518
519
520
521
522
    if (opt.compute_psd)
    {
      write_named_series(opt.outfile_psd,
                         "total power spectral density of signals",
                         frequencies,
                         PSD_vector,
                         opt.verbose,
                         opt.overwrite);
    } // if (opt.compute_psd)
  } // if (opt.compute_psd || opt.compute_npsd || opt.compute_transfer ||
    //  //      opt.compute_coherency)
523

524
525
526
} // main()

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