croposp.cc 14.9 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
#include <psdxx/psd.h>
47
48
49
50
51

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

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

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

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

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

/* ====================================================================== */
95

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

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

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

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

143
144
145
146
147
148
149
int main(int iargc, char* argv[])
{

  // define usage information
  char usage_text[]=
  {
    CROPOSP_VERSION "\n"
150
151
    "usage: croposp [-verbose] [-itype f] [-trim] [-datetolerance t]" "\n"
    "               [-label p]" "\n"
152
    "               [-psd f] [-npsd f] [-transfer f] [-coherency f]" "\n"
153
    "       file [t:sel] [f:format] [n:label] [file [t:s] [f:f] [n:l]] [...]\n"
154
155
156
157
158
159
    "   or: croposp --help|-h" "\n"
  };

  // define full help text
  char help_text[]=
  {
160
    "-verbose           be verbose\n"
161
    "-DEBUG             produce debugging output\n"
162
163
164
165
166
    "-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"
167
168
169
170
171
172
173
174
    "-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"
175
176
    "-log n             map averages to logarithmic scale with \"n\"\n"
    "                   samples pre decade\n"
177
    "\n"
178
179
180
181
182
183
184
185
186
187
188
    "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"
189
    "keys to be appended to file names if required:\n"
190
191
    "t:sel        trace selection\n"
    "f:format     file format\n"
192
    "n:label      label to be used in plot legend\n"
193
194
195
196
197
198
199
200
201
  };

  // define commandline options
  using namespace tfxx::cmdline;
  static Declare options[]= 
  {
    // 0: print help
    {"help",arg_no,"-"},
    // 1: verbose mode
202
203
204
    {"verbose",arg_no,"-"},
    // 2: default input file format
    {"itype",arg_yes,"sff"},
205
206
207
208
    // 3: trim time series
    {"trim",arg_no,"-"},
    // 4: default input file format
    {"datetolerance",arg_yes,"0."},
209
210
    // 5: debug mode
    {"DEBUG",arg_no,"-"},
211
212
    // 6: pattern for trace label
    {"pattern",arg_yes,"%S:%C:%A:%I"},
213
214
215
216
217
218
219
220
    // 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,"-"},
221
222
    // 11: compute coherency
    {"log",arg_yes,"1"},
223
224
225
226
227
    {NULL}
  };

  // define command line keys for input files
  static const char tracekey[]="t";
228
  static const char formatkey[]="f";
229
  static const char labelkey[]="n";
230
231

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

234
235
236
237
238
  /* ====================================================================== */
  /* analyze command line
   * --------------------
   */

239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
  // 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;
254
255
256
257
    cerr << endl;
    cerr << "References" << endl
      <<    "----------" << endl;
    cerr << reference_sleeman_et_al_2006 << endl;
258
259
    exit(0);
  }
260
261
262
263

  Options opt;
  opt.verbose=cmdline.optset(1);
  opt.inputformat=cmdline.string_arg(2);
264
265
  opt.trim=cmdline.optset(3);
  opt.datetolerance=cmdline.double_arg(4);
266
  opt.debug=cmdline.optset(5);
267
  opt.labelpattern=cmdline.string_arg(6);
268

269
270
271
272
273
274
275
276
277
  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);

278
279
280
281
282
283
  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");

284
285
  TFXX_assert(((opt.datetolerance >= 0.) && (opt.datetolerance <=1.)),
              "datetolerance is out of accepted range");
286
287
288
289
290
  
  // extract commandline arguments
  TFXX_assert(cmdline.extra(), "missing input file");
  tfxx::cmdline::Tparsed arguments=parse_cmdline(cmdline, cmdlinekeys);

291
292
293
294
295
  /* ====================================================================== */
  /* read data, prepare labels and trim time series if requested
   * -----------------------------------------------------------
   */

296
297
298
299
300
  // read data files
  if (opt.verbose) 
  {
    cout << "Read time series data" << endl;
  }
301
302
303
304
305
306
307
308
309
310
311
312
313
314
  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;
  }

315
  // collect traces
316
  TCollection collection_of_series;
317
  TMetaVector vector_of_metadata;
318
319
320
321
322
323
324
325
  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())
    {
326
      collection_of_series.push_back(*i_trace);
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
      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);
349
350
351
352
353
354
355
356
357
      ++i_trace;
    }
    ++i_file;
  }

  // report traces
  if (opt.verbose)
  {
    cout << "Time series to be analyzed:" << endl;
358
359
    report_collection(collection_of_series,
                      vector_of_metadata);
360
361
  }

362
363
364
  // trim time series
  if (opt.trim)
  {
365
366
    TFXX_assert(collection_of_series.overlap(),
                "time series do not overlap");
367
368
    collection_of_series.trim_to_date();
    collection_of_series.trim_to_nsamples();
369
370
371
372
373

    // report traces
    if (opt.verbose)
    {
      cout << "Time series after trimming:" << endl;
374
375
      report_collection(collection_of_series,
                        vector_of_metadata, opt.debug);
376
    }
377
378
379
380
381
382
383
384
  }

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

385
386
  /* ====================================================================== */

387
388
389
390
391
  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");

392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
  // 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;
    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;
      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())
  } // if (opt.compute_psd)

450
451
452
} // main()

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