sigfit.cc 23.8 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$
thomas.forbriger's avatar
thomas.forbriger committed
7
8
9
10
11
12
 * \author Thomas Forbriger
 * \date 28/01/2004
 * 
 * fit signal by trial-signals
 * 
 * Copyright (c) 2004 by Thomas Forbriger (BFO Schiltach) 
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
 *
 * ----
 * 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. 
 * 
 * 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, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 * ----
thomas.forbriger's avatar
thomas.forbriger committed
29
30
31
 * 
 * REVISIONS and CHANGES 
 *  - 28/01/2004   V1.0   Thomas Forbriger
thomas.forbriger's avatar
thomas.forbriger committed
32
33
 *  - 24/02/2004   V1.1   introduced search range stabilization
 *                        more output values
34
 *                 V1.2   correction: searchfac is a factor!
35
 *  - 16/03/2004   V1.3   add trend and offset fit
thomas.forbriger's avatar
thomas.forbriger committed
36
37
38
39
40
 *  - 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
41
 *                        - provide exponential decay
thomas.forbriger's avatar
thomas.forbriger committed
42
 *  - 02/08/2007   V1.5   explain results line
43
 *  - 11/02/2011   V1.6   add definition of traces
Daniel Armbruster's avatar
Daniel Armbruster committed
44
 *  - 20/02/2012   V1.7   read and write any file format (damb)
thomas.forbriger's avatar
thomas.forbriger committed
45
46
47
48
 * 
 * ============================================================================
 */
#define SIGFIT_VERSION \
49
  "SIGFIT   V1.6   fit signal by trial-signals"
thomas.forbriger's avatar
thomas.forbriger committed
50
#define SIGFIT_CVSID \
thomas.forbriger's avatar
thomas.forbriger committed
51
  "$Id$"
thomas.forbriger's avatar
thomas.forbriger committed
52
53
54
55
56

#include <fstream>
#include <iostream>
#include <list>
#include <vector>
thomas.forbriger's avatar
thomas.forbriger committed
57
#include <cmath>
thomas.forbriger's avatar
thomas.forbriger committed
58
#include <tfxx/commandline.h>
Daniel Armbruster's avatar
Daniel Armbruster committed
59
60
61
#include <tfxx/xcmdline.h>
#include <tfxx/rangelist.h>
#include <tfxx/rangestring.h>
thomas.forbriger's avatar
thomas.forbriger committed
62
#include <tfxx/error.h>
Daniel Armbruster's avatar
Daniel Armbruster committed
63
64
#include <datrwxx/readany.h>
#include <datrwxx/writeany.h>
65
#include <tsxx/tsxx.h>
thomas.forbriger's avatar
thomas.forbriger committed
66
67
68
#include <tsxx/innerproduct.h>
#include <aff/array.h>
#include <aff/dump.h>
69
#include <aff/seriesoperators.h>
70
#include <linearxx/lapackxx.h>
thomas.forbriger's avatar
thomas.forbriger committed
71

Daniel Armbruster's avatar
Daniel Armbruster committed
72
73
//typedef std::list<std::string> Tnamelist;
// double precision time series bundle which includes a WID2 header
74
typedef ts::TDsfftimeseries Tbundle;
Daniel Armbruster's avatar
Daniel Armbruster committed
75
typedef datrw::Tdseries Tseries;
thomas.forbriger's avatar
thomas.forbriger committed
76
typedef std::vector<Tbundle> Tbundlevec;
thomas.forbriger's avatar
thomas.forbriger committed
77
typedef aff::Array<Tbundle::Tseries::Tvalue> Tmatrix;
thomas.forbriger's avatar
thomas.forbriger committed
78
79
80
81
82
83
84
85

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
86
  double Tdate;
thomas.forbriger's avatar
thomas.forbriger committed
87
  bool truncate;
thomas.forbriger's avatar
thomas.forbriger committed
88
89
  std::string residualname;
  bool writeresidual;
thomas.forbriger's avatar
thomas.forbriger committed
90
91
  double searchrange;
  bool usesearchrange;
92
  bool equalsearch;
thomas.forbriger's avatar
thomas.forbriger committed
93
94
  bool fittrend, fitoffset, doskip, fitexp;
  double amptrend, ampoffset, skip, tcexp;
Daniel Armbruster's avatar
Daniel Armbruster committed
95
  std::string itype, otype;
thomas.forbriger's avatar
thomas.forbriger committed
96
97
}; // struct Options

thomas.forbriger's avatar
unclear    
thomas.forbriger committed
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
/*======================================================================*/
/* 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)
{
thomas.forbriger's avatar
thomas.forbriger committed
113
  TFXX_abort("removetre does not work");
thomas.forbriger's avatar
unclear    
thomas.forbriger committed
114
115
116
117
118
119
120
121
122
  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
123
  //double a,b;
thomas.forbriger's avatar
unclear    
thomas.forbriger committed
124
125
126
127
  sum /= s.size();
  for (int i=s.first(); i<=s.last(); ++i) { s(i) -= sum; }
}

thomas.forbriger's avatar
thomas.forbriger committed
128
129
130
131
132
133
134
135
// 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
136
137
138
/*======================================================================*/
/* main progran */

thomas.forbriger's avatar
thomas.forbriger committed
139
140
141
142
143
144
145
int main(int iargc, char* argv[])
{

  // define usage information
  char usage_text[]=
  {
    SIGFIT_VERSION "\n"
thomas.forbriger's avatar
thomas.forbriger committed
146
    "usage: sigfit [-v] [-Tdate v] [-truncate]" "\n"
thomas.forbriger's avatar
thomas.forbriger committed
147
148
    "              [-Sramp[=v]] [-Sconst[=v]] [-Sexp[=v]]" "\n"
    "              [-residual f]" "\n"
thomas.forbriger's avatar
thomas.forbriger committed
149
    "              [-searchrange[=r]] [-equalsearch] [-skip n]" "\n"
Daniel Armbruster's avatar
Daniel Armbruster committed
150
151
    "              [-itype F] [-otype F]" "\n"
    "              signal [t:T f:F] trial [t:T f:F] [trial [t:T f:F]...]" "\n"
thomas.forbriger's avatar
thomas.forbriger committed
152
153
154
155
156
157
158
159
160
    "   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
161
162
    "-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
163
164
165
    "-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
166
167
    "-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
168
    "-truncate    truncate all series to the same number of samples" "\n"
thomas.forbriger's avatar
thomas.forbriger committed
169
    "-residual f  write waveforms containing residual to file \"f\"" "\n"   
thomas.forbriger's avatar
thomas.forbriger committed
170
171
172
173
174
    "-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"
175
    "-equalsearch use same search range for each trial signal" "\n"
thomas.forbriger's avatar
thomas.forbriger committed
176
    "-skip n      skip n seconds at beginning of each trace" "\n"
Daniel Armbruster's avatar
Daniel Armbruster committed
177
178
    "-itype F     set input file formats to F (default: sff)" "\n"
    "-otype F     set output (residual) file format to F (default: sff)" "\n"
thomas.forbriger's avatar
thomas.forbriger committed
179
180
181
    "\n"
    "signal       signal to by explained by a linear combination" "\n"
    "             of trial signals" "\n"
Daniel Armbruster's avatar
Daniel Armbruster committed
182
183
184
185
186
187
188
    "trial        any set of trial signals" "\n"
    "\n"
    "File specific options:" "\n"
    "t:T          select specific traces from input file" "\n"
    "             T can be a list of traces like '1,4,5' or" "\n"
    "             a range like '6-19' or mixed like '5,8,12-17,20'" "\n"
    "f:F          specific file format (overrides -itype setting)" "\n"
thomas.forbriger's avatar
thomas.forbriger committed
189
190
191
    "\n"
    "The last line of the output is a summary of all results." "\n"
    "This line contains the following information separated by blanks:" "\n"
thomas.forbriger's avatar
thomas.forbriger committed
192
    "  channel station instrument julian_day date time" "\n"
thomas.forbriger's avatar
thomas.forbriger committed
193
194
195
196
197
198
199
    "  explained_rms signal_rms residual_rms" "\n"
    "  for each trial signal: regression_coefficient" "\n"
    "  rms_of_all_coefficients" "\n"
    "  for each trial signal: normalized_regression_coefficient" "\n"
    "  for each trial signal: regression_coefficient*trial_signal_rms" "\n"
    "  for each trial signal:" "\n"
    "    normalized_scalar_product_between_trial_signal_and_signal" "\n"
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
    "\n"
    "The output file (residual argument on the command line) will contain\n"
    "several traces (N test signals including DC offset and linear trend):\n"
    "  trace #1              signal, channel name as in input file\n"
    "  trace #2              \"syn\", synthetic trace, superposition of all\n"
    "                        test signals\n"
    "  trace #3              \"dif\", residual dif=signal-syn\n"
    "  traces #4 - #4+N-1    raw test signals\n"
    "  traces #4+N - #4+2N-1 scaled test signals in units of the final\n"
    "                        displaying their contribution to \"syn\"\n"
    "  trace #4+2N           \"csi\": csi=signal-cor\n"
    "  trace #4+2N+1         \"cys\": csy=syn-cor\n"
    "  trace #4+2N+2         \"cor\" internal corretion, i.e. the\n"
    "                        contributions of the DC offset and the linear\n"
    "                        trend to \"syn\"\n"
thomas.forbriger's avatar
thomas.forbriger committed
215
216
217
218
219
220
221
222
223
224
  };

  // 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
225
    // 2: date tolerance
thomas.forbriger's avatar
thomas.forbriger committed
226
    {"Tdate",arg_yes,"0."},
thomas.forbriger's avatar
thomas.forbriger committed
227
228
    // 3: truncate to a common number of samples
    {"truncate",arg_no,"-"},
thomas.forbriger's avatar
thomas.forbriger committed
229
230
    // 4: output SFF file
    {"residual",arg_yes,"-"},
thomas.forbriger's avatar
thomas.forbriger committed
231
232
    // 5: search range
    {"searchrange",arg_opt,"100."},
233
234
    // 6: equal search ranges
    {"equalsearch",arg_no,"-"},
235
    // 7: fit a trend
thomas.forbriger's avatar
thomas.forbriger committed
236
    {"Sramp",arg_opt,"1."},
237
    // 8: fit an offset
thomas.forbriger's avatar
thomas.forbriger committed
238
    {"Sconst",arg_opt,"1."},
thomas.forbriger's avatar
thomas.forbriger committed
239
240
    // 9: fit an offset
    {"skip",arg_yes,"0"},
thomas.forbriger's avatar
thomas.forbriger committed
241
242
    // 10: fit an offset
    {"Sexp",arg_opt,"1."},
Daniel Armbruster's avatar
Daniel Armbruster committed
243
244
245
246
    // 11: input file format
    {"itype",arg_yes,"sff"},
    // 12: output file format
    {"otype",arg_yes,"sff"},
thomas.forbriger's avatar
thomas.forbriger committed
247
248
249
    {NULL}
  };

Daniel Armbruster's avatar
Daniel Armbruster committed
250
251
252
253
254
255
256
257
258
259
260
261
262
  //! key to select traces
  const char* const tracekey="t";
  //! key to select file format
  const char* const formatkey="f";
  //! list of keys for filename specific parameters
  static const char* keys[]={
    //! select traces
    tracekey,
    //! set file format
    formatkey,
    0
  }; // const char* keys[]

thomas.forbriger's avatar
thomas.forbriger committed
263
264
265
266
267
268
269
270
271
272
273
274
275
  // 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))
  {
Daniel Armbruster's avatar
Daniel Armbruster committed
276
277
    cout << usage_text << endl;
    cout << help_text << endl;
thomas.forbriger's avatar
thomas.forbriger committed
278
279
280
281
282
283
    exit(0);
  }

  // read options
  Options opt;
  opt.verbose=cmdline.optset(1);
thomas.forbriger's avatar
thomas.forbriger committed
284
  opt.Tdate=cmdline.double_arg(2);
thomas.forbriger's avatar
thomas.forbriger committed
285
  opt.truncate=cmdline.optset(3);
thomas.forbriger's avatar
thomas.forbriger committed
286
287
  opt.writeresidual=cmdline.optset(4);
  opt.residualname=cmdline.string_arg(4);
thomas.forbriger's avatar
thomas.forbriger committed
288
289
  opt.usesearchrange=cmdline.optset(5);
  opt.searchrange=cmdline.double_arg(5);
290
  opt.equalsearch=cmdline.optset(6);
291
292
  opt.fittrend=cmdline.optset(7);
  opt.fitoffset=cmdline.optset(8);
thomas.forbriger's avatar
thomas.forbriger committed
293
294
  opt.amptrend=cmdline.double_arg(7);
  opt.ampoffset=cmdline.double_arg(8);
thomas.forbriger's avatar
thomas.forbriger committed
295
296
  opt.skip=cmdline.double_arg(9);
  opt.doskip=cmdline.optset(9);
thomas.forbriger's avatar
thomas.forbriger committed
297
298
  opt.fitexp=cmdline.optset(10);
  opt.tcexp=cmdline.double_arg(10);
Daniel Armbruster's avatar
Daniel Armbruster committed
299
300
  opt.itype=cmdline.string_arg(11);
  opt.otype=cmdline.string_arg(12);
thomas.forbriger's avatar
thomas.forbriger committed
301

thomas.forbriger's avatar
thomas.forbriger committed
302
303
  /*----------------------------------------------------------------------*/

thomas.forbriger's avatar
thomas.forbriger committed
304
  // read file names
thomas.forbriger's avatar
thomas.forbriger committed
305
  TFXX_assert(cmdline.extra(),"ERROR: missing signal file name!");
Daniel Armbruster's avatar
Daniel Armbruster committed
306
307
308
309
310
311
  //std::string signalname=cmdline.next();
  tfxx::cmdline::Tparsed infiles(tfxx::cmdline::parse_cmdline(cmdline, keys));
  //TFXX_assert((cmdline.extra()||opt.fittrend||opt.fitoffset),
  //            "ERROR: missing trial signal file name!");
  //Tnamelist namelist;
  //while (cmdline.extra()) { namelist.push_back(cmdline.next()); }
thomas.forbriger's avatar
thomas.forbriger committed
312

thomas.forbriger's avatar
thomas.forbriger committed
313
  // read signal file
Daniel Armbruster's avatar
Daniel Armbruster committed
314
315
316
317
318
319
320
321
322
323
324
325
  tfxx::cmdline::Tparsed::const_iterator infile(infiles.begin());
  TFXX_assert((infile!=infiles.end()||opt.fittrend||opt.fitoffset),
    "ERROR (sigfit): missing trial signal file name");
  std::string tracelist("1");
  if (infile->haskey(tracekey)) { tracelist=infile->value(tracekey); }
  tfxx::RangeListStepper<int> rls(tfxx::string::rangelist<int>(tracelist));
  TFXX_assert(rls.valid(), "Illegal tracelist");
  std::string signalname(infile->name);
  int signaltrace = rls.current();
  std::string signalformat(opt.itype);
  if (infile->haskey(formatkey)) { signalformat=infile->value(formatkey); }
 
thomas.forbriger's avatar
thomas.forbriger committed
326
  Tbundle signal;
Daniel Armbruster's avatar
Daniel Armbruster committed
327
328
  sff::FREE signalfilefree, signaltracefree;
  bool hasSignalFileFree;
thomas.forbriger's avatar
thomas.forbriger committed
329
  {
Daniel Armbruster's avatar
Daniel Armbruster committed
330
331
332
333
    if (opt.verbose) 
    { 
      cout << "sigfit: open signal file " << signalname << " of format " 
        << signalformat << endl;
thomas.forbriger's avatar
thomas.forbriger committed
334
    }
Daniel Armbruster's avatar
Daniel Armbruster committed
335
336
337
338
    std::ifstream ifs(signalname.c_str(), 
      datrw::ianystream::openmode(signalformat));
    datrw::ianystream is(ifs, signalformat);
    if (is.hasfree()) 
thomas.forbriger's avatar
thomas.forbriger committed
339
    { 
Daniel Armbruster's avatar
Daniel Armbruster committed
340
341
      is >> signalfilefree;
      hasSignalFileFree = true;
thomas.forbriger's avatar
thomas.forbriger committed
342
    }
Daniel Armbruster's avatar
Daniel Armbruster committed
343
344
345
346
347
348
349
    int itrace=1;
    while (itrace!=signaltrace) { is.skipseries(); ++itrace; }
    TFXX_assert(is.good(), "Illegal trace for signal input");
    if (opt.verbose) { cout << "read trace" << endl; }
    is >> signal;
    is >> signal.header;
    if (is.hasfree()) { is >> signaltracefree; }
thomas.forbriger's avatar
thomas.forbriger committed
350
  }
Daniel Armbruster's avatar
Daniel Armbruster committed
351
352
  ++infile;
 
thomas.forbriger's avatar
thomas.forbriger committed
353
  // read trial files
thomas.forbriger's avatar
thomas.forbriger committed
354
  Tbundlevec bundlevec;
Daniel Armbruster's avatar
Daniel Armbruster committed
355
  while (infile != infiles.end())
thomas.forbriger's avatar
thomas.forbriger committed
356
  {
Daniel Armbruster's avatar
Daniel Armbruster committed
357
358
359
360
361
362
363
364
365
366
    std::string itype(opt.itype);
    if (infile->haskey(formatkey)) { itype=infile->value(formatkey); }
    if (opt.verbose) 
    { 
      cout << "sigfit: Open trial signal file " << infile->name << " of format "
        << itype << endl;
    }
    bool selectedTraces=infile->haskey(tracekey);
    // no traces of file were selected -> read the entire file
    if (! selectedTraces)
thomas.forbriger's avatar
thomas.forbriger committed
367
    {
Daniel Armbruster's avatar
Daniel Armbruster committed
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
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
      std::ifstream ifs(infile->name.c_str(),
                          datrw::ianystream::openmode(itype));
      datrw::ianystream is(ifs, itype);
      int traceNum = 1;
      while (is.good())
      {
        if (opt.verbose)
        {
          cout << "sigfit: Reading trace #" << traceNum << " ..." << endl;
        }
        Tbundle bundle;
        is >> bundle;
        is >> bundle.header;
        bundlevec.push_back(bundle);
        if (opt.verbose) 
        { 
          cout << "  " << bundle.header.line().substr(0,69) << endl;
        }
        ++traceNum;
      }
    } else
    // traces were selected by rangelist -> read only selected traces
    {
      tfxx::RangeListStepper<int> rls(tfxx::string::rangelist<int>(
        infile->value(tracekey)));
      while (rls.valid())
      {
        std::ifstream ifs(infile->name.c_str(),
                            datrw::ianystream::openmode(itype));
        datrw::ianystream is(ifs, itype);
        int traceNum = 1;
        // spool to selected trace
        while (traceNum!=rls.current()) { is.skipseries(); ++traceNum; }
        TFXX_assert(is.good(), "Illegal trace selected.");

        if (opt.verbose)
        {
          cout << "sigfit: Reading trace #" << traceNum << " ..." << endl;
        }
        Tbundle bundle;
        is >> bundle;
        is >> bundle.header;
        bundlevec.push_back(bundle);
        if (opt.verbose) 
        { 
          cout << "  " << bundle.header.line().substr(0,69) << endl;
        }
        ++rls;
      }
thomas.forbriger's avatar
thomas.forbriger committed
417
    }
Daniel Armbruster's avatar
Daniel Armbruster committed
418
    ++infile;
thomas.forbriger's avatar
thomas.forbriger committed
419
  }
Daniel Armbruster's avatar
Daniel Armbruster committed
420

thomas.forbriger's avatar
thomas.forbriger committed
421
  unsigned int ntracesread=bundlevec.size();
thomas.forbriger's avatar
thomas.forbriger committed
422

423
424
425
426
427
428
429
  // 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
430
    synsignal.series()=opt.ampoffset;
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
    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
446
447
      synsignal(i)=2.*opt.amptrend*
        (i-(0.5*(synsignal.first()+synsignal.last())))/
448
449
450
451
452
453
454
455
456
        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
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
  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
479
480
  /*----------------------------------------------------------------------*/

thomas.forbriger's avatar
thomas.forbriger committed
481
482
483
484
  // truncate length if requested
  if (opt.truncate)
  {
    if (opt.verbose) { cout << "truncate signals if necessary..." << endl; }
485
    long int n=signal.last();
thomas.forbriger's avatar
thomas.forbriger committed
486
487
    for (Tbundlevec::const_iterator i=bundlevec.begin(); 
         i!=bundlevec.end(); i++)
488
489
    { n=n<i->last() ? n : i->last(); }
    if (signal.last()>n)
thomas.forbriger's avatar
thomas.forbriger committed
490
    {
491
492
      signal.setlastindex(n);
      signal.header.nsamples=signal.size();
thomas.forbriger's avatar
thomas.forbriger committed
493
494
495
496
497
498
      if (opt.verbose) { cout << "truncate signal to "
        << signal.header.nsamples << " samples" << endl; }
    }
    for (Tbundlevec::iterator i=bundlevec.begin(); 
         i!=bundlevec.end(); i++)
    { 
499
      if (i->last()>n)
thomas.forbriger's avatar
thomas.forbriger committed
500
      {
501
502
        i->setlastindex(n);
        i->header.nsamples=signal.size();
thomas.forbriger's avatar
thomas.forbriger committed
503
504
505
506
507
508
509
510
511
        if (opt.verbose) { cout << "truncate trial signal" << endl
          << i->header.line() << endl; }
      }
    }
  }

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

  // check header consistency
thomas.forbriger's avatar
thomas.forbriger committed
512
513
  sff::WID2compare compare(sff::Fnsamples | sff::Fdt | sff::Fdate);
  compare.setdatetolerance(opt.Tdate);
thomas.forbriger's avatar
thomas.forbriger committed
514
515
516
517
518
519
520
521
522
523
524
525
526
  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
527

thomas.forbriger's avatar
thomas.forbriger committed
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
  /*----------------------------------------------------------------------*/

  // 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
567
  /*----------------------------------------------------------------------*/
thomas.forbriger's avatar
thomas.forbriger committed
568
  double signalrms=ts::rms(signal);
thomas.forbriger's avatar
thomas.forbriger committed
569
570
571
572
573
  
  // 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
574
  Tmatrix Matrix(N,N), rhs(N), coeff(N), normproduct(N), trialrms(N);
575
  double fulltrialrms=0.;
thomas.forbriger's avatar
thomas.forbriger committed
576
577
578
579
  for (int i=1; i<=N; ++i)
  {
    for (int k=i; k<=N; ++k)
    {
580
      Matrix(i,k)=ts::innerproduct(bundlevec[i-1], bundlevec[k-1]);
thomas.forbriger's avatar
thomas.forbriger committed
581
582
      Matrix(k,i)=Matrix(i,k);
    }
583
    rhs(i)=ts::innerproduct(bundlevec[i-1], signal);
thomas.forbriger's avatar
thomas.forbriger committed
584
    trialrms(i)=ts::rms(bundlevec[i-1]);
585
    fulltrialrms+=(trialrms(i)*trialrms(i));
thomas.forbriger's avatar
thomas.forbriger committed
586
587
    normproduct(i)=rhs(i)/(signalrms*trialrms(i)*signal.size());
  }
588
  fulltrialrms=sqrt(fulltrialrms/double(N));
thomas.forbriger's avatar
thomas.forbriger committed
589
590
591
592
593
594
595
596
597

  // add stabilization
  if (opt.usesearchrange)
  {
    if (opt.verbose) 
    {
      cout << "stabilize fit:" << endl; 
      cout << "  relative search range: " << opt.searchrange << endl;
    }
598
    double searchfac=signal.size()*signalrms;
thomas.forbriger's avatar
thomas.forbriger committed
599
600
    for (int i=1; i<=N; ++i)
    {
601
602
603
604
605
      double range;
      if (opt.equalsearch)
      { range=opt.searchrange*signalrms/fulltrialrms; }
      else
      { range=opt.searchrange*signalrms/trialrms(i); }
thomas.forbriger's avatar
thomas.forbriger committed
606
607
608
609
610
      if (opt.verbose)
      {
        cout << "  search range for trial series " << i 
          << " is " << range << endl;
      }
611
      Matrix(i,i)+=pow((searchfac/range),2);
thomas.forbriger's avatar
thomas.forbriger committed
612
    }
thomas.forbriger's avatar
thomas.forbriger committed
613
  }
614

615
  // solve
616
  coeff=linear::lapack::dposv(Matrix,rhs);
617
618
619

  // set up synthetics
  Tbundle synthetics;
620
621
  synthetics=Tseries(signal.shape());
  synthetics.series()=0.;
622
  for (int i=1; i<=N; ++i)
623
  { synthetics += coeff(i) * bundlevec[i-1]; }
624
625
626
  synthetics.header=signal.header;
  synthetics.header.channel="synt";

thomas.forbriger's avatar
thomas.forbriger committed
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
  // 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;

647
648
  // set up residual
  Tbundle residual;
649
650
  residual=Tseries(signal.shape());
  residual=signal-synthetics;
651
652
  residual.header=signal.header;
  residual.header.channel="diff";
thomas.forbriger's avatar
thomas.forbriger committed
653
654
655
656
657
658

  // write waveforms
  if (opt.writeresidual)
  {
    if (opt.verbose) 
    {
Daniel Armbruster's avatar
Daniel Armbruster committed
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
      cout << "write residual waveform to " << opt.residualname 
        << " in format " << opt.otype << endl;
    }
    std::ofstream ofs(opt.residualname.c_str(),
                        datrw::oanystream::openmode(opt.otype));
    datrw::oanystream os(ofs, opt.otype);
    if (os.handlesfilefree())
    {
      sff::FREE residualFileFree;
      residualFileFree.append(SIGFIT_VERSION); 
      residualFileFree.append(SIGFIT_CVSID); 
      if (hasSignalFileFree)
      {
        residualFileFree.append("Signal file free block:");
        residualFileFree.append(signalfilefree);
      }
      os << residualFileFree;
thomas.forbriger's avatar
thomas.forbriger committed
676
    }
Daniel Armbruster's avatar
Daniel Armbruster committed
677
678
679
680
681
682
    os << signal.header;
    os << signal.series();
    os << synthetics.header;
    os << synthetics.series();
    os << residual.header;
    os << residual.series();
thomas.forbriger's avatar
thomas.forbriger committed
683
684
    for (Tbundlevec::const_iterator i=bundlevec.begin(); 
         i!=bundlevec.end(); i++)
Daniel Armbruster's avatar
Daniel Armbruster committed
685
    { os << i->header; os << *i; }
thomas.forbriger's avatar
thomas.forbriger committed
686
687
688
689
690
691
692
    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);
Daniel Armbruster's avatar
Daniel Armbruster committed
693
      os << psyn.header; 
thomas.forbriger's avatar
thomas.forbriger committed
694
695
      os << psyn;
    }
Daniel Armbruster's avatar
Daniel Armbruster committed
696
697
698
699
700
701
    os << corrsignal.header; 
    os << corrsignal.series();
    os << corrsynthetics.header;
    os << corrsynthetics.series();
    os << correction.header;
    os << correction.series();
thomas.forbriger's avatar
thomas.forbriger committed
702
  }
703
704
  
  // calculate rms values
705
  double residualrms=ts::rms(residual);
thomas.forbriger's avatar
thomas.forbriger committed
706
  double explained=1.-(residualrms/signalrms);
thomas.forbriger's avatar
thomas.forbriger committed
707
708
709
710
711
712
713
714
715

  // 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; }
716
  
thomas.forbriger's avatar
thomas.forbriger committed
717
718
719
720
  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
721
  cout << "   instrument: " << signal.header.instype << endl;
thomas.forbriger's avatar
thomas.forbriger committed
722
723
724
725
  cout << "    signalrms: " << signalrms << endl;
  cout << "  residualrms: " << residualrms << endl;
  cout << "explained rms: " << explained << endl;
  cout << " coefficients: ";
726
727
728
  for (int i=1; i<=N; ++i)
  { cout << coeff(i) << "  "; }
  cout << endl;
thomas.forbriger's avatar
thomas.forbriger committed
729
730
731
732
733
  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
734
735
736
  cout << " coefficients x rms: ";
  for (int i=1; i<=N; ++i)
  { 
thomas.forbriger's avatar
thomas.forbriger committed
737
738
739
740
741
742
743
    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
744
745
  }
  cout << endl;
746

thomas.forbriger's avatar
thomas.forbriger committed
747
748
749
  // reportline
  cout << signal.header.station << " "
    << signal.header.channel << " "
thomas.forbriger's avatar
thomas.forbriger committed
750
    << signal.header.instype << " "
thomas.forbriger's avatar
thomas.forbriger committed
751
752
753
754
755
756
757
758
    << 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
759
760
761
762
763
764
765
766
767
  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
768
769
  for (int i=1; i<=N; ++i)
  { 
thomas.forbriger's avatar
thomas.forbriger committed
770
    cout << formatfloat(normproduct(i)); 
thomas.forbriger's avatar
thomas.forbriger committed
771
  }
thomas.forbriger's avatar
thomas.forbriger committed
772
  cout << endl;
thomas.forbriger's avatar
thomas.forbriger committed
773
774
775
}

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