sigfit.cc 11 KB
Newer Older
thomas.forbriger's avatar
thomas.forbriger committed
1
2
3
4
5
/*! \file sigfit.cc
 * \brief fit signal by trial-signals
 * 
 * ----------------------------------------------------------------------------
 * 
thomas.forbriger's avatar
thomas.forbriger committed
6
 * $Id: sigfit.cc,v 1.9 2004-02-15 20:43:09 tforb Exp $
thomas.forbriger's avatar
thomas.forbriger committed
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
 * \author Thomas Forbriger
 * \date 28/01/2004
 * 
 * fit signal by trial-signals
 * 
 * Copyright (c) 2004 by Thomas Forbriger (BFO Schiltach) 
 * 
 * REVISIONS and CHANGES 
 *  - 28/01/2004   V1.0   Thomas Forbriger
 * 
 * ============================================================================
 */
#define SIGFIT_VERSION \
  "SIGFIT   V1.0   fit signal by trial-signals"
#define SIGFIT_CVSID \
thomas.forbriger's avatar
thomas.forbriger committed
22
  "$Id: sigfit.cc,v 1.9 2004-02-15 20:43:09 tforb Exp $"
thomas.forbriger's avatar
thomas.forbriger committed
23
24
25
26
27
28

#include <fstream>
#include <iostream>
#include <list>
#include <vector>
#include <tfxx/commandline.h>
thomas.forbriger's avatar
thomas.forbriger committed
29
#include <tfxx/error.h>
thomas.forbriger's avatar
thomas.forbriger committed
30
#include <sffxx.h>
31
32
#include <sffostream.h>
#include <tsxx/tsxx.h>
thomas.forbriger's avatar
thomas.forbriger committed
33
34
35
#include <tsxx/innerproduct.h>
#include <aff/array.h>
#include <aff/dump.h>
36
#include <aff/seriesoperators.h>
37
#include <linearxx/lapackxx.h>
thomas.forbriger's avatar
thomas.forbriger committed
38
39

typedef std::list<std::string> Tnamelist;
40
typedef ts::TDsfftimeseries Tbundle;
thomas.forbriger's avatar
thomas.forbriger committed
41
42
typedef Tbundle::Tseries Tseries;
typedef std::vector<Tbundle> Tbundlevec;
thomas.forbriger's avatar
thomas.forbriger committed
43
typedef aff::Array<Tbundle::Tseries::Tvalue> Tmatrix;
thomas.forbriger's avatar
thomas.forbriger committed
44
45
46
47
48
49
50
51

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

// structure to store commandline options
struct Options {
  bool verbose;
thomas.forbriger's avatar
thomas.forbriger committed
52
  double Tdate;
thomas.forbriger's avatar
thomas.forbriger committed
53
  bool truncate;
thomas.forbriger's avatar
thomas.forbriger committed
54
55
  std::string residualname;
  bool writeresidual;
thomas.forbriger's avatar
thomas.forbriger committed
56
57
}; // struct Options

thomas.forbriger's avatar
thomas.forbriger committed
58
59
60
61
62
63
64
65
// formatted number output
std::string formatfloat(const double &v)
{
  char seq[20];
  std::sprintf(seq, "%8.5f ", v);
  return(std::string(seq));
}

thomas.forbriger's avatar
thomas.forbriger committed
66
67
68
69
70
71
72
int main(int iargc, char* argv[])
{

  // define usage information
  char usage_text[]=
  {
    SIGFIT_VERSION "\n"
thomas.forbriger's avatar
thomas.forbriger committed
73
74
    "usage: sigfit [-v] [-Tdate v] [-truncate]" "\n"
    "              signal trial [trial ...]" "\n"
thomas.forbriger's avatar
thomas.forbriger committed
75
76
77
78
79
80
81
82
83
84
85
    "   or: sigfit --help|-h" "\n"
  };

  // define full help text
  char help_text[]=
  {
    SIGFIT_CVSID "\n"
    "\n"
    "-v           be verbose" "\n"
    "-Sramp       add a ramp to the set of trial signals" "\n"
    "-Sconst      add a constant to the set of trial signals" "\n"
thomas.forbriger's avatar
thomas.forbriger committed
86
87
    "-Tdate v     tolerance for comparison of date of first sample" "\n"
    "             v give the tolerance in units of the sampling interval" "\n"
thomas.forbriger's avatar
thomas.forbriger committed
88
    "-truncate    truncate all series to the same number of samples" "\n"
thomas.forbriger's avatar
thomas.forbriger committed
89
    "-residual f  write waveforms containing residual to file \"f\"" "\n"   
thomas.forbriger's avatar
thomas.forbriger committed
90
91
92
93
94
95
96
97
98
99
100
101
102
103
    "\n"
    "signal       signal to by explained by a linear combination" "\n"
    "             of trial signals" "\n"
    "trial        any set of trial signals"
  };

  // define commandline options
  using namespace tfxx::cmdline;
  static Declare options[]= 
  {
    // 0: print help
    {"help",arg_no,"-"},
    // 1: verbose mode
    {"v",arg_no,"-"},
thomas.forbriger's avatar
thomas.forbriger committed
104
    // 2: date tolerance
thomas.forbriger's avatar
thomas.forbriger committed
105
    {"Tdate",arg_yes,"0."},
thomas.forbriger's avatar
thomas.forbriger committed
106
107
    // 3: truncate to a common number of samples
    {"truncate",arg_no,"-"},
thomas.forbriger's avatar
thomas.forbriger committed
108
109
    // 4: output SFF file
    {"residual",arg_yes,"-"},
thomas.forbriger's avatar
thomas.forbriger committed
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
    {NULL}
  };

  // 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;
    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; }
*/

  // read options
  Options opt;
  opt.verbose=cmdline.optset(1);
thomas.forbriger's avatar
thomas.forbriger committed
158
  opt.Tdate=cmdline.double_arg(2);
thomas.forbriger's avatar
thomas.forbriger committed
159
  opt.truncate=cmdline.optset(3);
thomas.forbriger's avatar
thomas.forbriger committed
160
161
  opt.writeresidual=cmdline.optset(4);
  opt.residualname=cmdline.string_arg(4);
thomas.forbriger's avatar
thomas.forbriger committed
162

thomas.forbriger's avatar
thomas.forbriger committed
163
164
  /*----------------------------------------------------------------------*/

thomas.forbriger's avatar
thomas.forbriger committed
165
  // read file names
thomas.forbriger's avatar
thomas.forbriger committed
166
167
168
169
  TFXX_assert(cmdline.extra(),"ERROR: missing signal file name!");
  std::string signalname=cmdline.next();

  TFXX_assert(cmdline.extra(),"ERROR: missing trial signal file name!");
thomas.forbriger's avatar
thomas.forbriger committed
170
171
172
  Tnamelist namelist;
  while (cmdline.extra()) { namelist.push_back(cmdline.next()); }

thomas.forbriger's avatar
thomas.forbriger committed
173
174
175
176
177
178
179
180
181
182
183
184
  // read signal file
  Tbundle signal;
  {
    if (opt.verbose)
    { cout << "read signal file \"" << signalname << "\"" << endl; }
    std::ifstream is(signalname.c_str());
    sff::FileHeader fileheader(is);
    if (opt.verbose)
    { cout << "  read trace" << endl; }
    sff::InputWaveform<Tseries> inputwaveform(is);
    sff::TraceHeader traceheader=inputwaveform.header();
    signal.header=traceheader.wid2();
185
    signal=inputwaveform.series();
thomas.forbriger's avatar
thomas.forbriger committed
186
187
188
189
190
191
192
193
194
195
196
197
198
199
    bool last=traceheader.last();
    std::string wid2line=signal.header.line();
    if (opt.verbose)
    { 
      cout << "  " << wid2line.substr(0,69) << endl; 
      if (!last)
      {
        cout << "there are more traces in that files" << endl 
          << "ignoring them..." << endl;
      }
    }
  }
  
  // read trial files
thomas.forbriger's avatar
thomas.forbriger committed
200
  if (opt.verbose)
thomas.forbriger's avatar
thomas.forbriger committed
201
  { cout << namelist.size() << " trial data files to read" << endl; }
thomas.forbriger's avatar
thomas.forbriger committed
202
203
204
205
  Tbundlevec bundlevec;
  for (Tnamelist::const_iterator i=namelist.begin(); i!=namelist.end(); i++)
  {
    if (opt.verbose)
thomas.forbriger's avatar
thomas.forbriger committed
206
    { cout << "read trial data file \"" << *i << "\"" << endl; }
thomas.forbriger's avatar
thomas.forbriger committed
207
208
209
210
211
212
    std::ifstream is(i->c_str());
    sff::FileHeader fileheader(is);
    bool last=false;
    while(!last)
    {
      if (opt.verbose)
thomas.forbriger's avatar
thomas.forbriger committed
213
      { cout << "  read trial signal #" << bundlevec.size()+1 << endl; }
thomas.forbriger's avatar
thomas.forbriger committed
214
215
216
217
      sff::InputWaveform<Tseries> inputwaveform(is);
      sff::TraceHeader traceheader=inputwaveform.header();
      Tbundle bundle;
      bundle.header=traceheader.wid2();
218
      bundle=inputwaveform.series();
thomas.forbriger's avatar
thomas.forbriger committed
219
220
221
222
223
224
225
      bundlevec.push_back(bundle);
      last=traceheader.last();
      std::string wid2line=bundle.header.line();
      if (opt.verbose)
      { cout << "  " << wid2line.substr(0,69) << endl; }
    }
  }
thomas.forbriger's avatar
thomas.forbriger committed
226
227
228

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

thomas.forbriger's avatar
thomas.forbriger committed
229
230
231
232
  // truncate length if requested
  if (opt.truncate)
  {
    if (opt.verbose) { cout << "truncate signals if necessary..." << endl; }
233
    long int n=signal.last();
thomas.forbriger's avatar
thomas.forbriger committed
234
235
    for (Tbundlevec::const_iterator i=bundlevec.begin(); 
         i!=bundlevec.end(); i++)
236
237
    { n=n<i->last() ? n : i->last(); }
    if (signal.last()>n)
thomas.forbriger's avatar
thomas.forbriger committed
238
    {
239
240
      signal.setlastindex(n);
      signal.header.nsamples=signal.size();
thomas.forbriger's avatar
thomas.forbriger committed
241
242
243
244
245
246
      if (opt.verbose) { cout << "truncate signal to "
        << signal.header.nsamples << " samples" << endl; }
    }
    for (Tbundlevec::iterator i=bundlevec.begin(); 
         i!=bundlevec.end(); i++)
    { 
247
      if (i->last()>n)
thomas.forbriger's avatar
thomas.forbriger committed
248
      {
249
250
        i->setlastindex(n);
        i->header.nsamples=signal.size();
thomas.forbriger's avatar
thomas.forbriger committed
251
252
253
254
255
256
257
258
259
        if (opt.verbose) { cout << "truncate trial signal" << endl
          << i->header.line() << endl; }
      }
    }
  }

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

  // check header consistency
thomas.forbriger's avatar
thomas.forbriger committed
260
261
  sff::WID2compare compare(sff::Fnsamples | sff::Fdt | sff::Fdate);
  compare.setdatetolerance(opt.Tdate);
thomas.forbriger's avatar
thomas.forbriger committed
262
263
264
265
266
267
268
269
270
271
272
273
274
  if (opt.verbose) { cout << "checking consistency..." << endl; }
  for (Tbundlevec::const_iterator i=bundlevec.begin(); i!=bundlevec.end(); i++)
  {
    if (!compare (i->header,signal.header))
    {
      cerr << "ERROR: header signature mismatch:" << endl;
      cerr << "signal:" << endl;
      cerr << signal.header.line();
      cerr << "trial signal:" << endl;
      cerr << i->header.line();
      TFXX_abort("baling out...");
    }
  }
thomas.forbriger's avatar
thomas.forbriger committed
275
276
277
278
279
280
281
282
283
284
285
286

  /*----------------------------------------------------------------------*/
  
  // set up system of linear equations
  if (opt.verbose) { cout << "set up system of linear equations..." << endl; }
  int N=bundlevec.size();
  if (opt.verbose) { cout << "system is of size " << N << "x" << N << endl; }
  Tmatrix Matrix(N,N), rhs(N), coeff(N);
  for (int i=1; i<=N; ++i)
  {
    for (int k=i; k<=N; ++k)
    {
287
      Matrix(i,k)=ts::innerproduct(bundlevec[i-1], bundlevec[k-1]);
thomas.forbriger's avatar
thomas.forbriger committed
288
289
      Matrix(k,i)=Matrix(i,k);
    }
290
    rhs(i)=ts::innerproduct(bundlevec[i-1], signal);
thomas.forbriger's avatar
thomas.forbriger committed
291
  }
292

293
  // solve
294
  coeff=linear::lapack::dposv(Matrix,rhs);
295
296
297

  // set up synthetics
  Tbundle synthetics;
298
299
  synthetics=Tseries(signal.shape());
  synthetics.series()=0.;
300
  for (int i=1; i<=N; ++i)
301
  { synthetics += coeff(i) * bundlevec[i-1]; }
302
303
304
305
306
  synthetics.header=signal.header;
  synthetics.header.channel="synt";

  // set up residual
  Tbundle residual;
307
308
  residual=Tseries(signal.shape());
  residual=signal-synthetics;
309
310
  residual.header=signal.header;
  residual.header.channel="diff";
thomas.forbriger's avatar
thomas.forbriger committed
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327

  // write waveforms
  if (opt.writeresidual)
  {
    if (opt.verbose) 
    {
      cout << "write residual waveform to " 
      << opt.residualname << endl;
    }
    std::ofstream ofs(opt.residualname.c_str());
    sff::SFFostream<Tseries> os(ofs);
    os << signal;
    os << synthetics;
    os << residual;
    for (Tbundlevec::const_iterator i=bundlevec.begin(); 
         i!=bundlevec.end(); i++)
    { os << *i; }
thomas.forbriger's avatar
thomas.forbriger committed
328
329
330
331
332
333
334
335
336
    for (int i=1; i<=N; ++i)
    {
      Tbundle psyn=bundlevec[i-1];
      psyn=bundlevec[i-1].series()*coeff(i);
      char chan[6];
      std::sprintf(chan,"s%1.1d", i);
      psyn.header.channel=std::string(chan);
      os << psyn;
    }
thomas.forbriger's avatar
thomas.forbriger committed
337
  }
338
339
  
  // calculate rms values
340
341
  double signalrms=ts::rms(signal);
  double residualrms=ts::rms(residual);
thomas.forbriger's avatar
thomas.forbriger committed
342
  double explained=1.-(residualrms/signalrms);
343
  
thomas.forbriger's avatar
thomas.forbriger committed
344
345
346
347
  std::string timestring(signal.header.date.timestring());
  cout << "         time: " << timestring.substr(4,21) << endl;
  cout << "      channel: " << signal.header.channel << endl;
  cout << "      station: " << signal.header.station << endl;
thomas.forbriger's avatar
thomas.forbriger committed
348
  cout << "   instrument: " << signal.header.instype << endl;
thomas.forbriger's avatar
thomas.forbriger committed
349
350
351
352
  cout << "    signalrms: " << signalrms << endl;
  cout << "  residualrms: " << residualrms << endl;
  cout << "explained rms: " << explained << endl;
  cout << " coefficients: ";
353
354
355
  for (int i=1; i<=N; ++i)
  { cout << coeff(i) << "  "; }
  cout << endl;
thomas.forbriger's avatar
thomas.forbriger committed
356
357
358
359
360
361
  cout << " coefficients x rms: ";
  for (int i=1; i<=N; ++i)
  { 
    cout << coeff(i)*ts::rms(bundlevec[i-1]) << "  "; 
  }
  cout << endl;
362

thomas.forbriger's avatar
thomas.forbriger committed
363
364
365
  // reportline
  cout << signal.header.station << " "
    << signal.header.channel << " "
thomas.forbriger's avatar
thomas.forbriger committed
366
    << signal.header.instype << " "
thomas.forbriger's avatar
thomas.forbriger committed
367
368
369
370
371
372
373
374
    << timestring << " ";
  cout << formatfloat(explained);
  cout << formatfloat(signalrms);
  cout << formatfloat(residualrms);
  for (int i=1; i<=N; ++i)
  { 
    cout << formatfloat(coeff(i));
  }
thomas.forbriger's avatar
thomas.forbriger committed
375
376
377
378
  for (int i=1; i<=N; ++i)
  { 
    cout << formatfloat(coeff(i)*ts::rms(bundlevec[i-1])); 
  }
thomas.forbriger's avatar
thomas.forbriger committed
379
  cout << endl;
thomas.forbriger's avatar
thomas.forbriger committed
380
381
382
}

/* ----- END OF sigfit.cc ----- */