sigfit.cc 17.9 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.19 2004-09-24 15:36:33 tforb Exp $
thomas.forbriger's avatar
thomas.forbriger committed
7
8
9
10
11
12
13
14
15
 * \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
thomas.forbriger's avatar
thomas.forbriger committed
16
17
 *  - 24/02/2004   V1.1   introduced search range stabilization
 *                        more output values
18
 *                 V1.2   correction: searchfac is a factor!
19
 *  - 16/03/2004   V1.3   add trend and offset fit
thomas.forbriger's avatar
thomas.forbriger committed
20
21
22
23
24
 *  - 24/09/2004   V1.4   
 *                        - allow to skip samples if this is helpful to remove
 *                          transient filter loading response
 *                        - remove correction from synthetics and signal
 *                          and make this process transparent
thomas.forbriger's avatar
thomas.forbriger committed
25
 *                        - provide exponential decay
thomas.forbriger's avatar
thomas.forbriger committed
26
27
28
29
 * 
 * ============================================================================
 */
#define SIGFIT_VERSION \
thomas.forbriger's avatar
thomas.forbriger committed
30
  "SIGFIT   V1.4   fit signal by trial-signals"
thomas.forbriger's avatar
thomas.forbriger committed
31
#define SIGFIT_CVSID \
thomas.forbriger's avatar
thomas.forbriger committed
32
  "$Id: sigfit.cc,v 1.19 2004-09-24 15:36:33 tforb Exp $"
thomas.forbriger's avatar
thomas.forbriger committed
33
34
35
36
37

#include <fstream>
#include <iostream>
#include <list>
#include <vector>
thomas.forbriger's avatar
thomas.forbriger committed
38
#include <cmath>
thomas.forbriger's avatar
thomas.forbriger committed
39
#include <tfxx/commandline.h>
thomas.forbriger's avatar
thomas.forbriger committed
40
#include <tfxx/error.h>
thomas.forbriger's avatar
thomas.forbriger committed
41
#include <sffxx.h>
42
43
#include <sffostream.h>
#include <tsxx/tsxx.h>
thomas.forbriger's avatar
thomas.forbriger committed
44
45
46
#include <tsxx/innerproduct.h>
#include <aff/array.h>
#include <aff/dump.h>
47
#include <aff/seriesoperators.h>
48
#include <linearxx/lapackxx.h>
thomas.forbriger's avatar
thomas.forbriger committed
49
50

typedef std::list<std::string> Tnamelist;
51
typedef ts::TDsfftimeseries Tbundle;
thomas.forbriger's avatar
thomas.forbriger committed
52
53
typedef Tbundle::Tseries Tseries;
typedef std::vector<Tbundle> Tbundlevec;
thomas.forbriger's avatar
thomas.forbriger committed
54
typedef aff::Array<Tbundle::Tseries::Tvalue> Tmatrix;
thomas.forbriger's avatar
thomas.forbriger committed
55
56
57
58
59
60
61
62

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
63
  double Tdate;
thomas.forbriger's avatar
thomas.forbriger committed
64
  bool truncate;
thomas.forbriger's avatar
thomas.forbriger committed
65
66
  std::string residualname;
  bool writeresidual;
thomas.forbriger's avatar
thomas.forbriger committed
67
68
  double searchrange;
  bool usesearchrange;
69
  bool equalsearch;
thomas.forbriger's avatar
thomas.forbriger committed
70
71
  bool fittrend, fitoffset, doskip, fitexp;
  double amptrend, ampoffset, skip, tcexp;
thomas.forbriger's avatar
thomas.forbriger committed
72
73
}; // struct Options

thomas.forbriger's avatar
unclear    
thomas.forbriger committed
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
/*======================================================================*/
/* functions */

// remove average
void removeavg(Tseries& s)
{
  double avg=0.;
  for (int i=s.first(); i<=s.last(); ++i) { avg += s(i); }
  avg /= s.size();
  for (int i=s.first(); i<=s.last(); ++i) { s(i) -= avg; }
}

// remove trend
void removetre(Tseries& s)
{
  double sum=0.;
  double jser=0.;
  double jqser=0.;
  for (int i=s.first(); i<=s.last(); ++i) 
  { 
    sum += s(i); 
    jser += i;
    jqser += i*i;
  }
thomas.forbriger's avatar
thomas.forbriger committed
98
  //double a,b;
thomas.forbriger's avatar
unclear    
thomas.forbriger committed
99
100
101
102
  sum /= s.size();
  for (int i=s.first(); i<=s.last(); ++i) { s(i) -= sum; }
}

thomas.forbriger's avatar
thomas.forbriger committed
103
104
105
106
107
108
109
110
// 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
unclear    
thomas.forbriger committed
111
112
113
/*======================================================================*/
/* main progran */

thomas.forbriger's avatar
thomas.forbriger committed
114
115
116
117
118
119
120
int main(int iargc, char* argv[])
{

  // define usage information
  char usage_text[]=
  {
    SIGFIT_VERSION "\n"
thomas.forbriger's avatar
thomas.forbriger committed
121
    "usage: sigfit [-v] [-Tdate v] [-truncate]" "\n"
thomas.forbriger's avatar
thomas.forbriger committed
122
123
    "              [-Sramp[=v]] [-Sconst[=v]] [-Sexp[=v]]" "\n"
    "              [-residual f]" "\n"
thomas.forbriger's avatar
thomas.forbriger committed
124
    "              [-searchrange[=r]] [-equalsearch] [-skip n]" "\n"
thomas.forbriger's avatar
thomas.forbriger committed
125
    "              signal trial [trial ...]" "\n"
thomas.forbriger's avatar
thomas.forbriger committed
126
127
128
129
130
131
132
133
134
    "   or: sigfit --help|-h" "\n"
  };

  // define full help text
  char help_text[]=
  {
    SIGFIT_CVSID "\n"
    "\n"
    "-v           be verbose" "\n"
thomas.forbriger's avatar
thomas.forbriger committed
135
136
    "-Sramp[=v]   add a ramp to the set of trial signals with amplitude v" "\n"
    "-Sconst[=v]  add a constant (of value v) to the set of trial signals" "\n"
thomas.forbriger's avatar
thomas.forbriger committed
137
138
139
    "-Sexp[=v]    add a exponential decay (with time constant v relative" "\n"
    "             to the length of the time series) to the set of trial" "\n"
    "             signals" "\n"
thomas.forbriger's avatar
thomas.forbriger committed
140
141
    "-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
142
    "-truncate    truncate all series to the same number of samples" "\n"
thomas.forbriger's avatar
thomas.forbriger committed
143
    "-residual f  write waveforms containing residual to file \"f\"" "\n"   
thomas.forbriger's avatar
thomas.forbriger committed
144
145
146
147
148
    "-searchrange[=r] set search range to stabilize the fit in cases" "\n"
    "             where a trial series does not (or almost not)" "\n"
    "             contribute to the signal" "\n"
    "             the search range r is given in units of" "\n"
    "             rms(signal)/rms(trial signals)" "\n"
149
    "-equalsearch use same search range for each trial signal" "\n"
thomas.forbriger's avatar
thomas.forbriger committed
150
    "-skip n      skip n seconds at beginning of each trace" "\n"
thomas.forbriger's avatar
thomas.forbriger committed
151
152
153
154
155
156
157
158
159
160
161
162
163
164
    "\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
165
    // 2: date tolerance
thomas.forbriger's avatar
thomas.forbriger committed
166
    {"Tdate",arg_yes,"0."},
thomas.forbriger's avatar
thomas.forbriger committed
167
168
    // 3: truncate to a common number of samples
    {"truncate",arg_no,"-"},
thomas.forbriger's avatar
thomas.forbriger committed
169
170
    // 4: output SFF file
    {"residual",arg_yes,"-"},
thomas.forbriger's avatar
thomas.forbriger committed
171
172
    // 5: search range
    {"searchrange",arg_opt,"100."},
173
174
    // 6: equal search ranges
    {"equalsearch",arg_no,"-"},
175
    // 7: fit a trend
thomas.forbriger's avatar
thomas.forbriger committed
176
    {"Sramp",arg_opt,"1."},
177
    // 8: fit an offset
thomas.forbriger's avatar
thomas.forbriger committed
178
    {"Sconst",arg_opt,"1."},
thomas.forbriger's avatar
thomas.forbriger committed
179
180
    // 9: fit an offset
    {"skip",arg_yes,"0"},
thomas.forbriger's avatar
thomas.forbriger committed
181
182
    // 10: fit an offset
    {"Sexp",arg_opt,"1."},
thomas.forbriger's avatar
thomas.forbriger committed
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
    {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);
  }

  // read options
  Options opt;
  opt.verbose=cmdline.optset(1);
thomas.forbriger's avatar
thomas.forbriger committed
207
  opt.Tdate=cmdline.double_arg(2);
thomas.forbriger's avatar
thomas.forbriger committed
208
  opt.truncate=cmdline.optset(3);
thomas.forbriger's avatar
thomas.forbriger committed
209
210
  opt.writeresidual=cmdline.optset(4);
  opt.residualname=cmdline.string_arg(4);
thomas.forbriger's avatar
thomas.forbriger committed
211
212
  opt.usesearchrange=cmdline.optset(5);
  opt.searchrange=cmdline.double_arg(5);
213
  opt.equalsearch=cmdline.optset(6);
214
215
  opt.fittrend=cmdline.optset(7);
  opt.fitoffset=cmdline.optset(8);
thomas.forbriger's avatar
thomas.forbriger committed
216
217
  opt.amptrend=cmdline.double_arg(7);
  opt.ampoffset=cmdline.double_arg(8);
thomas.forbriger's avatar
thomas.forbriger committed
218
219
  opt.skip=cmdline.double_arg(9);
  opt.doskip=cmdline.optset(9);
thomas.forbriger's avatar
thomas.forbriger committed
220
221
  opt.fitexp=cmdline.optset(10);
  opt.tcexp=cmdline.double_arg(10);
thomas.forbriger's avatar
thomas.forbriger committed
222

thomas.forbriger's avatar
thomas.forbriger committed
223
224
  /*----------------------------------------------------------------------*/

thomas.forbriger's avatar
thomas.forbriger committed
225
  // read file names
thomas.forbriger's avatar
thomas.forbriger committed
226
227
228
  TFXX_assert(cmdline.extra(),"ERROR: missing signal file name!");
  std::string signalname=cmdline.next();

229
230
  TFXX_assert((cmdline.extra()||opt.fittrend||opt.fitoffset),
              "ERROR: missing trial signal file name!");
thomas.forbriger's avatar
thomas.forbriger committed
231
232
233
  Tnamelist namelist;
  while (cmdline.extra()) { namelist.push_back(cmdline.next()); }

thomas.forbriger's avatar
thomas.forbriger committed
234
235
236
237
238
239
240
241
242
  // 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; }
thomas.forbriger's avatar
thomas.forbriger committed
243
244
245
246
247
248
249
250
    sff::InputWaveform<Tseries> inputwaveform;
    try {
      inputwaveform.read(is);
    } 
    catch(...) {
      cerr << "ERROR (sigfit): read failed!" << endl;
      throw;
    }
thomas.forbriger's avatar
thomas.forbriger committed
251
252
    sff::TraceHeader traceheader=inputwaveform.header();
    signal.header=traceheader.wid2();
253
    signal=inputwaveform.series();
thomas.forbriger's avatar
thomas.forbriger committed
254
255
256
257
258
259
260
261
262
263
264
265
266
267
    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
268
  if (opt.verbose)
thomas.forbriger's avatar
thomas.forbriger committed
269
  { cout << namelist.size() << " trial data files to read" << endl; }
thomas.forbriger's avatar
thomas.forbriger committed
270
271
272
273
  Tbundlevec bundlevec;
  for (Tnamelist::const_iterator i=namelist.begin(); i!=namelist.end(); i++)
  {
    if (opt.verbose)
thomas.forbriger's avatar
thomas.forbriger committed
274
    { cout << "read trial data file \"" << *i << "\"" << endl; }
thomas.forbriger's avatar
thomas.forbriger committed
275
276
277
278
279
280
    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
281
      { cout << "  read trial signal #" << bundlevec.size()+1 << endl; }
thomas.forbriger's avatar
thomas.forbriger committed
282
283
284
285
      sff::InputWaveform<Tseries> inputwaveform(is);
      sff::TraceHeader traceheader=inputwaveform.header();
      Tbundle bundle;
      bundle.header=traceheader.wid2();
286
      bundle=inputwaveform.series();
thomas.forbriger's avatar
thomas.forbriger committed
287
288
289
290
291
292
293
      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
294
  unsigned int ntracesread=bundlevec.size();
thomas.forbriger's avatar
thomas.forbriger committed
295

296
297
298
299
300
301
302
  // add synthetic trial files
  if (opt.fitoffset)
  {
    if (opt.verbose) { cout << "add offset to trial signals" << endl; }
    Tbundle synsignal;
    synsignal.header=signal.header;
    synsignal=Tseries(signal.shape());
thomas.forbriger's avatar
thomas.forbriger committed
303
    synsignal.series()=opt.ampoffset;
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
    synsignal.header.station="NSP";
    synsignal.header.channel="off";
    synsignal.header.auxid="NSP";
    synsignal.header.instype="NSP";
    bundlevec.push_back(synsignal);
  }

  if (opt.fittrend)
  {
    if (opt.verbose) { cout << "add trend to trial signals" << endl; }
    Tbundle synsignal;
    synsignal.header=signal.header;
    synsignal=Tseries(signal.shape());
    for (int i=synsignal.first(); i<=synsignal.last(); ++i)
    {
thomas.forbriger's avatar
thomas.forbriger committed
319
320
      synsignal(i)=2.*opt.amptrend*
        (i-(0.5*(synsignal.first()+synsignal.last())))/
321
322
323
324
325
326
327
328
329
        synsignal.size();
    }
    synsignal.header.station="NSP";
    synsignal.header.channel="ramp";
    synsignal.header.auxid="NSP";
    synsignal.header.instype="NSP";
    bundlevec.push_back(synsignal);
  }

thomas.forbriger's avatar
thomas.forbriger committed
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
  if (opt.fitexp)
  {
    if (opt.verbose) 
    {
      cout << "add exponential decay to trial signals" << endl; 
    }
    Tbundle synsignal;
    synsignal.header=signal.header;
    synsignal=Tseries(signal.shape());
    for (int i=synsignal.first(); i<=synsignal.last(); ++i)
    {
      double argval=-1.*double(i-synsignal.first())/
        (double(synsignal.size())*opt.tcexp);
      synsignal(i)=std::exp(argval);
    }
    synsignal.header.station="NSP";
    synsignal.header.channel="exp";
    synsignal.header.auxid="NSP";
    synsignal.header.instype="NSP";
    bundlevec.push_back(synsignal);
  }

thomas.forbriger's avatar
thomas.forbriger committed
352
353
  /*----------------------------------------------------------------------*/

thomas.forbriger's avatar
thomas.forbriger committed
354
355
356
357
  // truncate length if requested
  if (opt.truncate)
  {
    if (opt.verbose) { cout << "truncate signals if necessary..." << endl; }
358
    long int n=signal.last();
thomas.forbriger's avatar
thomas.forbriger committed
359
360
    for (Tbundlevec::const_iterator i=bundlevec.begin(); 
         i!=bundlevec.end(); i++)
361
362
    { n=n<i->last() ? n : i->last(); }
    if (signal.last()>n)
thomas.forbriger's avatar
thomas.forbriger committed
363
    {
364
365
      signal.setlastindex(n);
      signal.header.nsamples=signal.size();
thomas.forbriger's avatar
thomas.forbriger committed
366
367
368
369
370
371
      if (opt.verbose) { cout << "truncate signal to "
        << signal.header.nsamples << " samples" << endl; }
    }
    for (Tbundlevec::iterator i=bundlevec.begin(); 
         i!=bundlevec.end(); i++)
    { 
372
      if (i->last()>n)
thomas.forbriger's avatar
thomas.forbriger committed
373
      {
374
375
        i->setlastindex(n);
        i->header.nsamples=signal.size();
thomas.forbriger's avatar
thomas.forbriger committed
376
377
378
379
380
381
382
383
384
        if (opt.verbose) { cout << "truncate trial signal" << endl
          << i->header.line() << endl; }
      }
    }
  }

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

  // check header consistency
thomas.forbriger's avatar
thomas.forbriger committed
385
386
  sff::WID2compare compare(sff::Fnsamples | sff::Fdt | sff::Fdate);
  compare.setdatetolerance(opt.Tdate);
thomas.forbriger's avatar
thomas.forbriger committed
387
388
389
390
391
392
393
394
395
396
397
398
399
  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
400

thomas.forbriger's avatar
thomas.forbriger committed
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
  /*----------------------------------------------------------------------*/

  // skip samples if requested
  if (opt.doskip) 
  {
    if (opt.verbose)
    {
      cout << "skip " << opt.skip << " seconds" << endl;
    }
    int nskip=int(floor(opt.skip/signal.header.dt));
    if (nskip>0)
    {
      libtime::TRelativeTime add=nskip*libtime::double2time(signal.header.dt);
      signal.setfirstindex(signal.first()+nskip);
      signal.header.date+=add;
      for (Tbundlevec::iterator i=bundlevec.begin(); 
           i!=bundlevec.end(); i++)
      {
        i->setfirstindex(i->first() + nskip);
        i->header.date += add;
      }
      if (opt.verbose)
      {
        cout << "skipped " << nskip << " samples ("
          << add.timestring() << ")" << endl;
      }
    }
    else
    {
      if (opt.verbose)
      {
        cout << "NOTICE: nothing to skip..." << endl;
        cout << "  seems like " 
          << opt.skip << " seconds to skip < " 
          << signal.header.dt << " seconds sampling interval" << endl;
      }
    }
  }

thomas.forbriger's avatar
thomas.forbriger committed
440
  /*----------------------------------------------------------------------*/
thomas.forbriger's avatar
thomas.forbriger committed
441
  double signalrms=ts::rms(signal);
thomas.forbriger's avatar
thomas.forbriger committed
442
443
444
445
446
  
  // 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; }
thomas.forbriger's avatar
thomas.forbriger committed
447
  Tmatrix Matrix(N,N), rhs(N), coeff(N), normproduct(N), trialrms(N);
448
  double fulltrialrms=0.;
thomas.forbriger's avatar
thomas.forbriger committed
449
450
451
452
  for (int i=1; i<=N; ++i)
  {
    for (int k=i; k<=N; ++k)
    {
453
      Matrix(i,k)=ts::innerproduct(bundlevec[i-1], bundlevec[k-1]);
thomas.forbriger's avatar
thomas.forbriger committed
454
455
      Matrix(k,i)=Matrix(i,k);
    }
456
    rhs(i)=ts::innerproduct(bundlevec[i-1], signal);
thomas.forbriger's avatar
thomas.forbriger committed
457
    trialrms(i)=ts::rms(bundlevec[i-1]);
458
    fulltrialrms+=(trialrms(i)*trialrms(i));
thomas.forbriger's avatar
thomas.forbriger committed
459
460
    normproduct(i)=rhs(i)/(signalrms*trialrms(i)*signal.size());
  }
461
  fulltrialrms=sqrt(fulltrialrms/double(N));
thomas.forbriger's avatar
thomas.forbriger committed
462
463
464
465
466
467
468
469
470

  // add stabilization
  if (opt.usesearchrange)
  {
    if (opt.verbose) 
    {
      cout << "stabilize fit:" << endl; 
      cout << "  relative search range: " << opt.searchrange << endl;
    }
471
    double searchfac=signal.size()*signalrms;
thomas.forbriger's avatar
thomas.forbriger committed
472
473
    for (int i=1; i<=N; ++i)
    {
474
475
476
477
478
      double range;
      if (opt.equalsearch)
      { range=opt.searchrange*signalrms/fulltrialrms; }
      else
      { range=opt.searchrange*signalrms/trialrms(i); }
thomas.forbriger's avatar
thomas.forbriger committed
479
480
481
482
483
      if (opt.verbose)
      {
        cout << "  search range for trial series " << i 
          << " is " << range << endl;
      }
484
      Matrix(i,i)+=pow((searchfac/range),2);
thomas.forbriger's avatar
thomas.forbriger committed
485
    }
thomas.forbriger's avatar
thomas.forbriger committed
486
  }
487

488
  // solve
489
  coeff=linear::lapack::dposv(Matrix,rhs);
490
491
492

  // set up synthetics
  Tbundle synthetics;
493
494
  synthetics=Tseries(signal.shape());
  synthetics.series()=0.;
495
  for (int i=1; i<=N; ++i)
496
  { synthetics += coeff(i) * bundlevec[i-1]; }
497
498
499
  synthetics.header=signal.header;
  synthetics.header.channel="synt";

thomas.forbriger's avatar
thomas.forbriger committed
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
  // set up correction
  Tbundle correction;
  correction=Tseries(signal.shape());
  correction.series()=0.;
  if (ntracesread<bundlevec.size())
  {
    for (int i=ntracesread; i<N; ++i)
    { correction += coeff(i+1) * bundlevec[i]; }
  }
  correction.header=signal.header;
  correction.header.channel="corr";
  Tbundle corrsignal;
  corrsignal.header=signal.header;
  corrsignal.header.channel="csig";
  corrsignal=signal-correction;
  Tbundle corrsynthetics;
  corrsynthetics.header=signal.header;
  corrsynthetics.header.channel="csyn";
  corrsynthetics=synthetics-correction;

520
521
  // set up residual
  Tbundle residual;
522
523
  residual=Tseries(signal.shape());
  residual=signal-synthetics;
524
525
  residual.header=signal.header;
  residual.header.channel="diff";
thomas.forbriger's avatar
thomas.forbriger committed
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542

  // 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
543
544
545
546
547
548
549
550
551
    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
552
553
554
    os << corrsignal;
    os << corrsynthetics;
    os << correction;
thomas.forbriger's avatar
thomas.forbriger committed
555
  }
556
557
  
  // calculate rms values
558
  double residualrms=ts::rms(residual);
thomas.forbriger's avatar
thomas.forbriger committed
559
  double explained=1.-(residualrms/signalrms);
thomas.forbriger's avatar
thomas.forbriger committed
560
561
562
563
564
565
566
567
568

  // calculate direction cosines
  Tmatrix normcoeff(N);
  double length=0.;
  for (int i=1; i<=N; ++i)
  { length += coeff(i)*coeff(i); }
  length=sqrt(length);
  for (int i=1; i<=N; ++i)
  { normcoeff(i) = coeff(i)/length; }
569
  
thomas.forbriger's avatar
thomas.forbriger committed
570
571
572
573
  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
574
  cout << "   instrument: " << signal.header.instype << endl;
thomas.forbriger's avatar
thomas.forbriger committed
575
576
577
578
  cout << "    signalrms: " << signalrms << endl;
  cout << "  residualrms: " << residualrms << endl;
  cout << "explained rms: " << explained << endl;
  cout << " coefficients: ";
579
580
581
  for (int i=1; i<=N; ++i)
  { cout << coeff(i) << "  "; }
  cout << endl;
thomas.forbriger's avatar
thomas.forbriger committed
582
583
584
585
586
  cout << " rms coeffic.: " << length << endl;
  cout << " normalized coefficients: ";
  for (int i=1; i<=N; ++i)
  { cout << normcoeff(i) << "  "; }
  cout << endl;
thomas.forbriger's avatar
thomas.forbriger committed
587
588
589
  cout << " coefficients x rms: ";
  for (int i=1; i<=N; ++i)
  { 
thomas.forbriger's avatar
thomas.forbriger committed
590
591
592
593
594
595
596
    cout << coeff(i)*trialrms(i) << "  "; 
  }
  cout << endl;
  cout << " normalized scalar product: ";
  for (int i=1; i<=N; ++i)
  { 
    cout << normproduct(i) << "  "; 
thomas.forbriger's avatar
thomas.forbriger committed
597
598
  }
  cout << endl;
599

thomas.forbriger's avatar
thomas.forbriger committed
600
601
602
  // reportline
  cout << signal.header.station << " "
    << signal.header.channel << " "
thomas.forbriger's avatar
thomas.forbriger committed
603
    << signal.header.instype << " "
thomas.forbriger's avatar
thomas.forbriger committed
604
605
606
607
608
609
610
611
    << 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
612
613
614
615
616
617
618
619
620
  cout << formatfloat(length);
  for (int i=1; i<=N; ++i)
  { 
    cout << formatfloat(normcoeff(i));
  }
  for (int i=1; i<=N; ++i)
  { 
    cout << formatfloat(coeff(i)*trialrms(i)); 
  }
thomas.forbriger's avatar
thomas.forbriger committed
621
622
  for (int i=1; i<=N; ++i)
  { 
thomas.forbriger's avatar
thomas.forbriger committed
623
    cout << formatfloat(normproduct(i)); 
thomas.forbriger's avatar
thomas.forbriger committed
624
  }
thomas.forbriger's avatar
thomas.forbriger committed
625
  cout << endl;
thomas.forbriger's avatar
thomas.forbriger committed
626
627
628
}

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