foutra.cc 40.7 KB
Newer Older
thomas.forbriger's avatar
thomas.forbriger committed
1
2
3
4
5
6
7
8
9
10
11
/*! \file foutra.cc
 * \brief Fourier transforms
 * 
 * ----------------------------------------------------------------------------
 * 
 * \author Thomas Forbriger
 * \date 25/07/2006
 * 
 * Fourier transforms
 * 
 * Copyright (c) 2006 by Thomas Forbriger (BFO Schiltach) 
12
13
 *
 * ----
thomas.forbriger's avatar
thomas.forbriger committed
14
 * FOUTRA is free software; you can redistribute it and/or modify
15
16
17
18
19
20
21
22
23
24
25
26
27
 * 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
28
29
30
 * 
 * REVISIONS and CHANGES 
 *  - 25/07/2006   V1.0   Thomas Forbriger
thomas.forbriger's avatar
thomas.forbriger committed
31
 *  - 12/09/2007   V1.1   this version produces correct spectra
thomas.forbriger's avatar
thomas.forbriger committed
32
33
34
 *  - 13/09/2007   V1.2   The program is tested against Walter's code
 *                        now provides processing details in verbose mode and
 *                        in trace FREE block
thomas.forbriger's avatar
thomas.forbriger committed
35
 *  - 17/09/2007   V1.4   provides:
thomas.forbriger's avatar
thomas.forbriger committed
36
37
38
 *                        - demean and detrend
 *                        - scaling to average in relative bandwidth
 *                        - ASCII output
thomas.forbriger's avatar
thomas.forbriger committed
39
 *                        - output ASCII table on logarithmic frequency scale
thomas.forbriger's avatar
thomas.forbriger committed
40
 *                        - average only ASCII output if requested
thomas.forbriger's avatar
thomas.forbriger committed
41
42
 *  - 26/06/2009   V1.5   additional feature:
 *                        - adjust number of samples to make FFT more efficient
43
 *  - 08/10/2010   V1.6   start testing foutra scaling
44
 *  - 12/12/2010   V1.7   implemented amplitude scaling for harmonic signal
thomas.forbriger's avatar
thomas.forbriger committed
45
 *  - 14/12/2010   V1.8   implemented segment averaging
46
47
 *                        distinguish between input series, output spectrum
 *                        and segment (analysis) series
thomas.forbriger's avatar
thomas.forbriger committed
48
49
 *  - 10/01/2011   V1.8b  corrected Fourier transformation formula in
 *                        doxygen documentation
50
 *  - 03/12/2014   V1.9   provide output file format selector
51
52
 *  - 08/01/2015   V1.10  FIX: when using -scalerbw take bandwidth from
 *                        opt.scaledecades rather than opt.decades
53
54
 *  - 29/01/2015   V1.11  FEATURE: provide calculation of n-th derivative
 *                        of time series
55
 *  - 02/04/2019          use new header file interface to libtsioxx
56
 *
57
 *  \note 08/01/2010:
58
59
60
61
 *  Scaling for foutra power spectrum was tested against theory.
 *  The integral over the power spectral density calculated by foutra over the
 *  whole frequency band (up to Nyquist frequency) provides the total power of
 *  the signal (i.e. the variance, i.e. the square of the rms value).
thomas.forbriger's avatar
thomas.forbriger committed
62
63
 *
 *  \sa \ref page_foutra
thomas.forbriger's avatar
thomas.forbriger committed
64
 * 
65
66
67
68
69
 *  \todo
 *  The use of series containers is messy. The same container is used for
 *  different purposes, once being used as a reference, than as a copy.
 *  This should be sorted out.
 *
thomas.forbriger's avatar
thomas.forbriger committed
70
71
72
 * ============================================================================
 */
#define FOUTRA_VERSION \
73
  "FOUTRA   V1.10   Fourier transforms"
thomas.forbriger's avatar
thomas.forbriger committed
74
75
76
77

#include <iostream>
#include <tfxx/commandline.h>
#include <aff/series.h>
thomas.forbriger's avatar
thomas.forbriger committed
78
79
#include <aff/iterator.h>
#include <aff/dump.h>
thomas.forbriger's avatar
thomas.forbriger committed
80
81
#include <aff/seriesoperators.h>
#include <aff/functions/avg.h>
thomas.forbriger's avatar
thomas.forbriger committed
82
83
#include <aff/functions/rms.h>
#include <aff/functions/sqrsum.h>
thomas.forbriger's avatar
thomas.forbriger committed
84
#include <aff/subarray.h>
thomas.forbriger's avatar
thomas.forbriger committed
85
#include <tsxx/tsxx.h>
thomas.forbriger's avatar
thomas.forbriger committed
86
#include <tsxx/tapers.h>
thomas.forbriger's avatar
thomas.forbriger committed
87
#include <tsxx/filter.h>
88
#include <tsxx/wid2timeseries.h>
89
#include <tsioxx/tsioxx.h>
thomas.forbriger's avatar
thomas.forbriger committed
90
91
92
93
94
95
#include <fstream>
#include <tfxx/error.h>
#include <tfxx/rangestring.h>
#include <tfxx/xcmdline.h>
#include <tfxx/misc.h>
#include <tfxx/handle.h>
96
#include <tfxx/seitosh.h>
97
#include <datrwxx/readany.h>
98
#include <datrwxx/writeany.h>
thomas.forbriger's avatar
thomas.forbriger committed
99
#include <fourier/fftwaff.h>
thomas.forbriger's avatar
thomas.forbriger committed
100
#include <sstream>
thomas.forbriger's avatar
thomas.forbriger committed
101

102
103
104
105
#include "foutra_options_brief_help.h"
#include "foutra_options_help.h"
#include "foutra_usage_help.h"

thomas.forbriger's avatar
thomas.forbriger committed
106
107
108
109
110
using std::cout;
using std::cerr;
using std::endl;

struct Options {
thomas.forbriger's avatar
thomas.forbriger committed
111
  bool verbose, overwrite, debug, asciiout, logascii;
112
  std::string inputformat, asciibase, outputformat;
thomas.forbriger's avatar
thomas.forbriger committed
113
  bool amplitudespectrum, powerspectrum, boxcartaper;
thomas.forbriger's avatar
thomas.forbriger committed
114
  bool avgconstbw, avgrelbw, avgasciionly;
thomas.forbriger's avatar
thomas.forbriger committed
115
  bool demean, detrend, scalerbw, adjustdivisor;
116
  bool derivative;
thomas.forbriger's avatar
thomas.forbriger committed
117
  int ndemean, ndetrend, divisor;
118
  double decades, scaledecades, asciidecades, nderivative;
thomas.forbriger's avatar
thomas.forbriger committed
119
  int avgsamples;
120
  bool reportrms, harmonicsignal, padzeroes;
121
  int padfactor, nsegments;
thomas.forbriger's avatar
thomas.forbriger committed
122
123
124
125
126
127
128
129
130
131
132
}; // struct Options

// values type to be used for samples
typedef double Tvalue;

// time series
typedef aff::Series<Tvalue> Tseries;

// full featured time series file
typedef ts::sff::File<Tseries> Tfile;

thomas.forbriger's avatar
thomas.forbriger committed
133
134
135
typedef ts::TDsfftimeseries Ttimeseries;
typedef Ttimeseries::Tseries Tseries;

thomas.forbriger's avatar
thomas.forbriger committed
136
137
typedef fourier::fft::DRFFTWAFF Tfft;

thomas.forbriger's avatar
thomas.forbriger committed
138
139
/*======================================================================*/

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

  // define usage information
144
  char version_text[]=
thomas.forbriger's avatar
thomas.forbriger committed
145
146
147
148
149
150
151
152
153
154
155
156
  {
    FOUTRA_VERSION "\n"
  };

  // 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
157
158
159
160
161
162
    // 2: overwrite mode
    {"o",arg_no,"-"},
    // 3: input format
    {"type",arg_yes,"sff"},
    // 4: debug mode
    {"D",arg_no,"-"},
thomas.forbriger's avatar
thomas.forbriger committed
163
164
165
166
    // 5: amplitude spectrum
    {"amplitude",arg_no,"-"},
    // 6: power spectrum
    {"power",arg_no,"-"},
thomas.forbriger's avatar
thomas.forbriger committed
167
168
    // 7: apply boxcar taper
    {"boxcar",arg_no,"-"},
thomas.forbriger's avatar
thomas.forbriger committed
169
170
171
172
    // 8: average over constant number of samples
    {"avg",arg_opt,"21"},
    // 9: average over constant relative bandwidth (in decades)
    {"rbw",arg_opt,"0.167"},
thomas.forbriger's avatar
thomas.forbriger committed
173
174
175
176
177
178
179
180
    // 10: remove average
    {"demean",arg_opt,"0"},
    // 11: remove trend
    {"detrend",arg_opt,"0"},
    // 12: scale to mean in relative bandwidth
    {"scalerbw",arg_opt,"0.167"},
    // 13: write result to ASCII files
    {"ASCII",arg_opt,"spectrum"},
thomas.forbriger's avatar
thomas.forbriger committed
181
182
    // 14: write result to ASCII files
    {"logascii",arg_opt,"0.167"},
thomas.forbriger's avatar
thomas.forbriger committed
183
184
    // 15: average ASCII output only
    {"avgascii",arg_no,"-"},
thomas.forbriger's avatar
thomas.forbriger committed
185
186
    // 16: divisor for number of samples
    {"divisor",arg_opt,"100"},
187
188
    // 17: report rms value
    {"rms",arg_no,"-"},
189
190
    // 18: report rms value
    {"harmonic",arg_no,"-"},
191
192
    // 19: pad time series
    {"pad",arg_yes,"1"},
193
194
    // 20: subdivide input series
    {"nsegments",arg_yes,"1"},
195
196
    // 21: extended help
    {"xhelp",arg_no,"-"},
197
198
    // 22: input format
    {"Type",arg_yes,"sff"},
199
200
    // 23: input format
    {"derivative",arg_opt,"1"},
thomas.forbriger's avatar
thomas.forbriger committed
201
202
203
    {NULL}
  };

thomas.forbriger's avatar
thomas.forbriger committed
204
205
206
207
208
  static const char tracekey[]="t";

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

thomas.forbriger's avatar
thomas.forbriger committed
209
210
211
  // no arguments? print usage...
  if (iargc<2) 
  {
212
213
    cerr << version_text << endl;
    cerr << foutra_options_brief_help << endl;
214
    cerr << tfxx::seitosh::repository_reference << endl;
thomas.forbriger's avatar
thomas.forbriger committed
215
216
217
218
219
220
221
    exit(0);
  }

  // collect options from commandline
  Commandline cmdline(iargc, argv, options);

  // help requested? print full help text...
222
  if (cmdline.optset(0) || cmdline.optset(21))
thomas.forbriger's avatar
thomas.forbriger committed
223
  {
224
225
226
227
    cerr << version_text << endl;
    cerr << foutra_options_brief_help << endl;
    cerr << foutra_options_help << endl;
    cerr << foutra_usage_help << endl;
228
    datrw::supported_data_types(cerr);
229
    if (cmdline.optset(21)) { datrw::online_help(cerr); }
230
    cerr << endl << tfxx::seitosh::repository_reference << endl;
thomas.forbriger's avatar
thomas.forbriger committed
231
232
233
    exit(0);
  }

thomas.forbriger's avatar
thomas.forbriger committed
234
235
236
237
238
  Options opt;
  opt.verbose=cmdline.optset(1);
  opt.overwrite=cmdline.optset(2);
  opt.inputformat=cmdline.string_arg(3);
  opt.debug=cmdline.optset(4);
thomas.forbriger's avatar
thomas.forbriger committed
239
  opt.amplitudespectrum=cmdline.optset(5);
thomas.forbriger's avatar
thomas.forbriger committed
240
241
  if (!opt.amplitudespectrum)
  { opt.powerspectrum=cmdline.optset(6); }
thomas.forbriger's avatar
thomas.forbriger committed
242
  opt.boxcartaper=cmdline.optset(7);
thomas.forbriger's avatar
thomas.forbriger committed
243
244
245
246
  opt.avgconstbw=cmdline.optset(8);
  opt.avgsamples=cmdline.int_arg(8);
  opt.avgrelbw=cmdline.optset(9);
  opt.decades=cmdline.double_arg(9);
thomas.forbriger's avatar
thomas.forbriger committed
247
248
249
250
251
252
253
254
  opt.demean=cmdline.optset(10);
  opt.ndemean=cmdline.int_arg(10);
  opt.detrend=cmdline.optset(11);
  opt.ndetrend=cmdline.int_arg(11);
  opt.scalerbw=cmdline.optset(12);
  opt.scaledecades=cmdline.double_arg(12);
  opt.asciiout=cmdline.optset(13);
  opt.asciibase=cmdline.string_arg(13);
thomas.forbriger's avatar
thomas.forbriger committed
255
256
  opt.logascii=cmdline.optset(14);
  opt.asciidecades=cmdline.double_arg(14);
thomas.forbriger's avatar
thomas.forbriger committed
257
  opt.avgasciionly=cmdline.optset(15);
thomas.forbriger's avatar
thomas.forbriger committed
258
259
  opt.adjustdivisor=cmdline.optset(16);
  opt.divisor=cmdline.int_arg(16);
260
  opt.reportrms=cmdline.optset(17);
261
  opt.harmonicsignal=cmdline.optset(18);
262
263
  opt.padzeroes=cmdline.optset(19);
  opt.padfactor=cmdline.int_arg(19);
264
  opt.nsegments=cmdline.int_arg(20);
265
  opt.outputformat=cmdline.string_arg(22);
266
267
  opt.derivative=cmdline.optset(23);
  opt.nderivative=cmdline.double_arg(23);
thomas.forbriger's avatar
thomas.forbriger committed
268
269

  TFXX_assert((opt.divisor > 0), "illegal value for argument divisor");
270
271
  TFXX_assert((opt.nsegments > 0), "illegal value for argument nsegments");
  TFXX_assert((opt.padfactor > 0), "illegal value for argument pad");
thomas.forbriger's avatar
thomas.forbriger committed
272
273

  if (opt.verbose)
274
  { cout << FOUTRA_VERSION << endl; }
thomas.forbriger's avatar
thomas.forbriger committed
275
276
277
278
279
280
281

  // extract commandline arguments
  TFXX_assert(cmdline.extra(), "missing output file");
  std::string outfile=cmdline.next();
  TFXX_assert(cmdline.extra(), "missing input file");
  tfxx::cmdline::Tparsed arguments=parse_cmdline(cmdline, cmdlinekeys);
  if ((arguments.size()>1) && opt.verbose)
thomas.forbriger's avatar
thomas.forbriger committed
282
  {
thomas.forbriger's avatar
thomas.forbriger committed
283
284
285
    cout << "NOTICE: file specific information (SRCE line and file FREE)" <<
      endl
      <<    "        will be taken from first file only!" << endl;
thomas.forbriger's avatar
thomas.forbriger committed
286
287
  }

288
289
290
291
  /*======================================================================*/
  // do not apply segmentation if not PSD or harmonic signal analysis
  if (!(opt.powerspectrum || opt.harmonicsignal)) { opt.nsegments=1; }

thomas.forbriger's avatar
thomas.forbriger committed
292
293
294
  /*======================================================================*/
  // create processing description
  sff::FREE processfree;
295
296
297
298
299
300
301
302
303
  if (opt.nsegments>1)
  {
    std::ostringstream freeline;
    freeline << "Input series is subdivided into "
      << opt.nsegments << " segments with 50% overlap";
    processfree.append(freeline.str());
    processfree.append("The analysis is done individually for all segments");
    processfree.append("The result is the average over all segments");
  }
thomas.forbriger's avatar
thomas.forbriger committed
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
  if (opt.demean)
  {
      std::ostringstream freeline;
      freeline << "The average over ";
      if (opt.ndemean >0)
      {
        freeline << opt.ndemean;
      }
      else
      {
        freeline << "all";
      }
      freeline << " samples is removed from all samples";
      processfree.append(freeline.str());
  }
  if (opt.detrend)
  {
      std::ostringstream freeline;
      freeline << "The trend over ";
      if (opt.ndetrend >0)
      {
        freeline << opt.ndetrend;
      }
      else
      {
        freeline << "all";
      }
      freeline << " samples is removed from all samples";
      processfree.append(freeline.str());
  }
thomas.forbriger's avatar
thomas.forbriger committed
334
335
336
337
338
339
340
341
  if (opt.boxcartaper)
  {
    processfree.append("A boxcar taper (i.e no taper) is applied to the data");
  }
  else
  {
    processfree.append("A Hanning taper is applied to the data");
  }
342
  if (opt.padzeroes) 
343
344
  {
    std::ostringstream freeline;
345
    freeline << "pad with zeroes increasing the total "
346
347
348
      << "number of samples by a factor of " << opt.padfactor;
    processfree.append(freeline.str());
  }
349
  if (opt.harmonicsignal || opt.amplitudespectrum || opt.powerspectrum)
thomas.forbriger's avatar
thomas.forbriger committed
350
351
  {
    processfree.append("An appropriately scaled FFT (libdrfftw) is applied");
352
353
354
355
356
357
358
359
360
    if (opt.derivative)
    {
      processfree.append("Calculate values for derivative with respect to");
      processfree.append("  time by multiplication of the Fourier coefficients");
      std::ostringstream freeline;
      freeline << "  with the angular frequency to the power of " 
        << opt.nderivative;
      processfree.append(freeline.str());
    }
thomas.forbriger's avatar
thomas.forbriger committed
361
362
363
364
365
366
367
368
369
  }
  if (opt.amplitudespectrum)
  {
    processfree.append("The result are coefficients of the amplitude spectrum");
    processfree.append("Units are: input units / Hz");
  }
  else if (opt.powerspectrum)
  {
    processfree.append("The result are coefficients of the power spectrum");
thomas.forbriger's avatar
thomas.forbriger committed
370
371
372
373
374
375
376
377
378
379
380
381
    if (opt.scalerbw) 
    {
      std::ostringstream freeline;
      freeline << "The result is the average power in " 
        << opt.scaledecades << " decades";
      processfree.append(freeline.str());
      processfree.append("Units are: input units squared");
    }
    else
    {
      processfree.append("Units are: input units squared / Hz");
    }
thomas.forbriger's avatar
thomas.forbriger committed
382
    if ((opt.avgconstbw || opt.avgrelbw) && !opt.avgasciionly)
thomas.forbriger's avatar
thomas.forbriger committed
383
    {
thomas.forbriger's avatar
thomas.forbriger committed
384
      processfree.append("The spectrum is smoothed with a boxcar over");
thomas.forbriger's avatar
thomas.forbriger committed
385
386
387
388
389
390
391
392
393
394
395
396
      if (opt.avgconstbw)
      {
        std::ostringstream freeline;
        freeline << opt.avgsamples << " samples";
        processfree.append(freeline.str());
      }
      else if (opt.avgrelbw)
      {
        std::ostringstream freeline;
        freeline << opt.decades << " decades";
        processfree.append(freeline.str());
      }
thomas.forbriger's avatar
thomas.forbriger committed
397
398
399
400
    }
    else
    {
      processfree.append("The spectrum is not smoothed");
thomas.forbriger's avatar
thomas.forbriger committed
401
      if (opt.avgasciionly && opt.asciiout)
thomas.forbriger's avatar
thomas.forbriger committed
402
403
404
405
406
407
      {
        processfree.append("The ASCII output is smoothed with a boxcar over");
        std::ostringstream freeline;
        freeline << opt.decades << " decades";
        processfree.append(freeline.str());
      }
thomas.forbriger's avatar
thomas.forbriger committed
408
409
    }
  }
410
411
412
413
414
415
  else if (opt.harmonicsignal)
  {
    processfree.append("The input is understood to contain harmonic signals.");
    processfree.append("The result are signal amplitudes for harmonic peaks.");
    processfree.append("Units are: input units");
  }
thomas.forbriger's avatar
thomas.forbriger committed
416
417
418
419
  else
  {
    processfree.append("The output is a tapered version of the time series");
  }
420
  if (opt.amplitudespectrum || opt.powerspectrum || opt.harmonicsignal)
thomas.forbriger's avatar
thomas.forbriger committed
421
422
423
424
425
426
427
428
  {
    processfree.append("The sampling interval provided in the WID2 line");
    processfree.append("specifies the sampling interval of the spectrum");
    processfree.append("in Hz. The first coefficient is that at 0 Hz.");
  }

  if (opt.verbose)
  {
429
    cout << "processing of each trace will take place as follows:" << endl;
thomas.forbriger's avatar
thomas.forbriger committed
430
431
432
433
434
    for(sff::FREE::Tlines::const_iterator I=processfree.lines.begin(); 
        I != processfree.lines.end(); ++I) 
    { cout << "  " << *I << std::endl; }
  }

thomas.forbriger's avatar
thomas.forbriger committed
435
436
  /*======================================================================*/
  // start processing
thomas.forbriger's avatar
thomas.forbriger committed
437
438
439
  
  // create taper
  ts::tapers::Hanning taper;
thomas.forbriger's avatar
thomas.forbriger committed
440
441
  // create FFTW processor
  Tfft fft;
thomas.forbriger's avatar
thomas.forbriger committed
442
443
444
445
446
447
448
449
450
451
452

  // open output file
  // ----------------
  if (opt.verbose) { cout << "open output file " << outfile << endl; }
  // check if output file exists and open
  if (!opt.overwrite)
  {
    std::ifstream file(outfile.c_str(),std::ios_base::in);
    TFXX_assert((!file.good()),"ERROR: output file exists!");
  }
  std::ofstream ofs(outfile.c_str());
453
  datrw::oanystream os(ofs, opt.outputformat);
thomas.forbriger's avatar
thomas.forbriger committed
454
455
456
457
458
459

  // prepare file FREE block
  sff::FREE filefree;
  filefree.append(FOUTRA_VERSION);
  // set flag to process header of first input file
  bool firstfile=true;
thomas.forbriger's avatar
thomas.forbriger committed
460
461
  // count output traces
  int otrace=0;
thomas.forbriger's avatar
thomas.forbriger committed
462
463
464
465
466
467
468
469
  // cycle through all input files
  // -----------------------------
  tfxx::cmdline::Tparsed::const_iterator infile=arguments.begin();
  while (infile != arguments.end())
  {
    // open input file
    if (opt.verbose) { cout << "open input file " << infile->name << endl; }
    std::ifstream ifs(infile->name.c_str());
470
    datrw::ianystream is(ifs, opt.inputformat);
thomas.forbriger's avatar
thomas.forbriger committed
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
    // handle file header
    if (firstfile)
    {
      if (is.hasfree()) 
      { 
        sff::FREE infilefree;
        is >> infilefree;
        filefree.append("block read from first input file:");
        filefree.append(infilefree);
      }
      os << filefree;
      if (is.hassrce())
      {
        sff::SRCE insrceline;
        is >> insrceline;
        os << insrceline;
      }
    }

    // cycle through traces of input file
    // ----------------------------------
    // setup trace selection
    typedef tfxx::RangeList<int> Trangelist;
    bool doselect=infile->haskey(tracekey);
    Trangelist traceranges=
      tfxx::string::rangelist<Trangelist::Tvalue>(infile->value(tracekey));
    int itrace=0;
    while (is.good())
    {
      ++itrace;
      if ((!doselect) || traceranges.contains(itrace))
      {
        TFXX_debug(opt.debug, "main", "process trace #" << itrace );
        if (opt.verbose)
        { std::cout << "  process trace #" << itrace << std::endl; }
506
507
        Tseries inputseries;
        is >> inputseries;
thomas.forbriger's avatar
thomas.forbriger committed
508
509
        sff::WID2 wid2;
        is >> wid2;
thomas.forbriger's avatar
thomas.forbriger committed
510
        TFXX_debug(opt.debug, "main", "  series and WID2 are read");
thomas.forbriger's avatar
thomas.forbriger committed
511

thomas.forbriger's avatar
thomas.forbriger committed
512
  /*----------------------------------------------------------------------*/
513
514
515
516
  // report rms value if requested
        if (opt.reportrms)
        {
          std::cout << "    rms-value of input time series: "
517
            << aff::func::rms(inputseries) << std::endl;
518
519
        }
  /*----------------------------------------------------------------------*/
thomas.forbriger's avatar
thomas.forbriger committed
520
  // process data
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
  //
  // inputseries: full series read from file
  // series: series for one segment to be analysed 
  //         the spectral representation of a single segment is also stored in
  //         this container
  // spectrum: resulting spectral representation
  //
  // in the previous (pre segmentation) version of foutra only "series" was
  // used; we will make use of the aff::Series reference semantics to use the
  // old code with no more modifications than necessary; such "series" will
  // be used for different things
  //
  // wid2.dt remains valid throughout the processing
  // wid2.nsamples is only correct for the inputseries
  //   use series.size() for other request for number of samples

        // prepare segmentation
        // --------------------
        int nsegsamples=inputseries.size();
        if (opt.nsegments>1)
        {
          nsegsamples=2*inputseries.size()/(opt.nsegments+1); 
        }
thomas.forbriger's avatar
thomas.forbriger committed
544
545
546
547
548
549
550

        // adjust length of time series
        if (opt.adjustdivisor)
        {
          if (opt.verbose)
          {
            std::cout << "    adjust divisor to " << opt.divisor << std::endl;
551
552
553
            std::cout << "    number of input samples:      " << wid2.nsamples
              <<std::endl;
            std::cout << "    number of segment samples:    " << nsegsamples
thomas.forbriger's avatar
thomas.forbriger committed
554
555
              <<std::endl;
          }
556
557
          int rest=nsegsamples % opt.divisor;
          nsegsamples=nsegsamples-rest;
thomas.forbriger's avatar
thomas.forbriger committed
558
559
          if (opt.verbose)
          {
560
            std::cout << "    number of processing samples: " << nsegsamples
thomas.forbriger's avatar
thomas.forbriger committed
561
562
563
              <<std::endl;
          }
        }
thomas.forbriger's avatar
thomas.forbriger committed
564

565
566
        // segment stride
        int segstride=inputseries.size()/(opt.nsegments+1);
thomas.forbriger's avatar
thomas.forbriger committed
567

568
569
570
571
572
573
574
575
576
577
578
579
        if (opt.verbose && (opt.nsegments>1))
        {
           cout << "    number of input samples:           " 
             << inputseries.size() << endl;
           cout << "    number of segments:                " 
             << opt.nsegments << endl;
           cout << "    number of samples in each segment: " 
             << nsegsamples << endl;
           cout << "    stride between segments:           "
             << segstride << endl;
           cout << "    overlap of proximate segments:     "
             << 100.*(nsegsamples-segstride)/nsegsamples << "%" << endl;
thomas.forbriger's avatar
thomas.forbriger committed
580
581
        }

582
583
584
585
586
587
        // length of time series segment
        double T=wid2.dt*nsegsamples;
        // length of taper window
        double Tw=T;
        // frequency sampling interval
        double df=1./T;
thomas.forbriger's avatar
thomas.forbriger committed
588

589
        if (opt.padzeroes)
590
591
592
593
594
        {
          T *= opt.padfactor;
          df /= opt.padfactor;
        }

595
        if (opt.verbose)
thomas.forbriger's avatar
thomas.forbriger committed
596
        {
597
598
599
600
601
602
          cout << "    duration T of each segment:     " 
            << T << "s" << endl;
          cout << "    duration Tw of signal window:   "
            << Tw << "s" << endl;
          cout << "    frequency sampling interval df: "
            << df << "Hz" << endl;
thomas.forbriger's avatar
thomas.forbriger committed
603
        }
604
605
606
607
608
609
610
611
612
613
            
        // the result will be collected in
        Tseries spectrum;

        // segment counter
        int isegment=0;

        // start actual processing
        // -----------------------
        while (isegment<opt.nsegments)
thomas.forbriger's avatar
thomas.forbriger committed
614
        {
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
          int firstindex=inputseries.first()+isegment*segstride;
          int lastindex=firstindex+nsegsamples-1;
          TFXX_debug(opt.debug, "main", "  segment index range: "
                     << firstindex << "-" << lastindex);
          TFXX_assert((firstindex >= inputseries.first()) &&
                      (firstindex <= inputseries.last()) &&
                      (lastindex >= inputseries.first()) &&
                      (lastindex <= inputseries.last()),
                      "ERROR: index out of range; program design error");
          Tseries segseries=aff::subarray(inputseries)(firstindex,lastindex);
          Tseries series=segseries;
          if (opt.nsegments>1) { series=segseries.copyout(); }
          series.shift(-series.first());

          // demean
          if (opt.demean) { 
            TFXX_debug(opt.debug, "main", "  demean");
            ts::filter::RemoveAverage demean(opt.ndemean);
            demean(ts::filter::Ttimeseries(series, wid2.dt));
thomas.forbriger's avatar
thomas.forbriger committed
634
          }
thomas.forbriger's avatar
thomas.forbriger committed
635

636
637
638
639
640
          // detrend
          if (opt.detrend) { 
            TFXX_debug(opt.debug, "main", "  detrend");
            ts::filter::RemoveTrend detrend(opt.ndetrend);
            detrend(ts::filter::Ttimeseries(series, wid2.dt));
thomas.forbriger's avatar
thomas.forbriger committed
641
          }
642
643
644
645
646

          // apply taper
          if (!opt.boxcartaper) { 
            TFXX_debug(opt.debug, "main", "  apply taper");
            taper.apply(series); 
thomas.forbriger's avatar
thomas.forbriger committed
647
648
          }

649
650
651
652
653
654
655
656
657
658
659
          // pad time series
          // Notice: padding the signals requires correct scaling of amplitudes
          // this is only implemented for harmonic signals up to now
          if (opt.padzeroes)
          {
            Tseries newseries(series.size()*opt.padfactor);
            newseries=0.;
            Tseries subseries(aff::subarray(newseries)(series.f(),series.l()));
            subseries.copyin(series);
            series=newseries;
          }
thomas.forbriger's avatar
thomas.forbriger committed
660

661
662
          // call FFT for amplitude or power spectrum or harmonic signal
          if (opt.amplitudespectrum || opt.powerspectrum || opt.harmonicsignal)
thomas.forbriger's avatar
thomas.forbriger committed
663
          {
664
665
666
667
668
669
670
671
            TFXX_debug(opt.debug, "main", "call FFT");
            Tfft::Tspectrum coeff=fft(series, wid2.dt);
            TFXX_debug(opt.debug, "main",
                       "returned from FFT; create output series");
            series=Tseries(coeff.size());
            TFXX_debug(opt.debug, "main", "create iterators");
            aff::Iterator<Tseries> S(series);
            aff::Browser<Tfft::Tspectrum> C(coeff);
672
673
674
675
676
677
678
679
680
681
682
683
684

            // take derivative if selected
            if (opt.derivative)
            {
              for (int i=coeff.f()+1; i<=coeff.l(); ++i)
              {
                double frequency=df*(i-coeff.f())*2*M_PI;
                double thepower=opt.nderivative;
                double factor=std::pow(frequency,thepower);
                coeff(i) *= factor;
              }
            }

685
686
            TFXX_debug(opt.debug, "main", "calculate square of modulus and copy");
            while(S.valid() && C.valid())
thomas.forbriger's avatar
thomas.forbriger committed
687
            {
688
689
              *S = C->real()*C->real()+C->imag()*C->imag();
              ++S; ++C;
thomas.forbriger's avatar
thomas.forbriger committed
690
            }
691
692
693
694
          }
            
          // take sqrt for amplitude spectrum
          if (opt.amplitudespectrum || opt.harmonicsignal)
thomas.forbriger's avatar
thomas.forbriger committed
695
          {
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
            TFXX_debug(opt.debug, "main", "calculate sqrt and scale");
            aff::Iterator<Tseries> S(series);
            while(S.valid())
            {
              *S = std::sqrt(*S);
              ++S;
            }
          }
          
          // apply scaling appropriate for harmonic signals
          if (opt.harmonicsignal)
          {
            // amplitude of Hanning taper instrument function equals 1
            // amplitude of Boxcar taper instrument function equals 2
            double fscaling;
            if (opt.boxcartaper) 
            { fscaling = 2./Tw; }
            else
            { fscaling = 4./Tw; }
            series *= fscaling;
          }
thomas.forbriger's avatar
thomas.forbriger committed
717

718
719
          // apply appropriate scaling and averaging for power spectrum
          if (opt.powerspectrum)
thomas.forbriger's avatar
thomas.forbriger committed
720
          {
721
722
723
            // scaling factor to adjust for taper effect
            double tapscaling=1.;
            if (opt.boxcartaper)
thomas.forbriger's avatar
thomas.forbriger committed
724
            {
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
              tapscaling=1.;
            }
            else
            {
              // scaling factor for Hanning
              // see below for derivation of this factor
              tapscaling=8./3.;
            }

            // we have an energy spectrum so far
            // adjust scaling factor to obtain signal power
            double scalingfactor=2.*tapscaling/Tw;

            // scale to relative bandwidth if requested
            if (opt.scalerbw)
            {
              TFXX_debug(opt.debug, "main", 
                         "scale to average in relative bandwidth");
743
744
              double bwfactor=(std::pow(10.,opt.scaledecades)-1.)/
                std::pow(10.,opt.scaledecades/2.);
745
746
747
748
              // cout << "bwfactor: " << bwfactor << endl;
              // cout << "index range: " 
              //   << series.f() << "-" << series.l() << endl;
              for (int k=series.f(); k<=series.l(); ++k)
thomas.forbriger's avatar
thomas.forbriger committed
749
              {
750
751
752
                // center frequency
                double fm=k*df;
                series(k) *= scalingfactor*bwfactor*fm;
thomas.forbriger's avatar
thomas.forbriger committed
753
              }
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
            } // if (opt.scalerbw)
            else
            {
              TFXX_debug(opt.debug, "main", "scale");
              series *= scalingfactor;
            } // if (opt.scalerbw) else

            // apply smoothing
            if (opt.avgconstbw && !opt.avgasciionly)
            {
              TFXX_debug(opt.debug, "main", "average over constant bandwidth");
              // create a copy
              Tseries p(series.size());
              p.copyin(series);
              int d=opt.avgsamples/2+1;
              for (int k=p.f(); k<=p.l(); ++k)
              {
                int k1=k-d;
                int k2=k+d;
                k1 = k1 < p.f() ? p.f() : k1;
                k2 = k2 > p.l() ? p.l() : k2;
                // cout << k1 << " " << k << " " << k2 << endl;
                TFXX_assert(k2>k1, "negative bandwidth for averaging")
                Tseries sub=aff::subarray(p)(k1,k2);
                series(k)=aff::func::avg(sub);
                /*
                series(k)=0;
                for (int j=sub.f(); j<=sub.l(); ++j)
                {
                  series(k) += sub(j);
                }
                series(k) /= double(sub.size());
                */
              } // for (int k=p.f(); k<=p.l(); ++k)
            } // if (opt.avgconstbw && !opt.avgasciionly)
            else if (opt.avgrelbw && !opt.avgasciionly)
            {
              TFXX_debug(opt.debug, "main", "average over relative bandwidth");
              // create a copy
              TFXX_debug(opt.debug, "main", "  create a copy");
              Tseries p(series.size());
              TFXX_debug(opt.debug, "main", "  copy in data");
              p.copyin(series);
              double bwfactor=std::sqrt(std::pow(10.,opt.decades));
              // cout << bwfactor << endl;
              TFXX_debug(opt.debug, "main", "  cycle over data");
              for (int k=p.f(); k<=p.l(); ++k)
              {
                // center frequency
                double fm=k*df;
                double f1=fm/bwfactor;
                double f2=fm*bwfactor;
                // cout << fm << " " << f1 << " " << f2 << endl;
                int k1=int(f1/df);
                int k2=int(f2/df)+1;
                k1 = k1 < p.f() ? p.f() : k1;
                k2 = k2 > p.l() ? p.l() : k2;
                TFXX_assert(k2>k1, "negative bandwidth for averaging")
                Tseries sub=aff::subarray(p)(k1,k2);
                series(k)=aff::func::avg(sub);
              } // for (int k=p.f(); k<=p.l(); ++k)
            } // else if (opt.avgrelbw && !opt.avgasciionly)
          } // if (opt.powerspectrum)

          // collect result
          if (isegment==0)
thomas.forbriger's avatar
thomas.forbriger committed
820
          {
821
822
823
824
825
826
827
828
829
            spectrum=series.copyout();
            TFXX_debug(opt.debug, "main", "copy result");
          }
          else
          {
            // spectrum += series;
            for (int i=spectrum.f(),
                 j=series.f();
                 i<=spectrum.l(); ++i, ++j)
thomas.forbriger's avatar
thomas.forbriger committed
830
            {
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
              TFXX_assert(j<=series.l(),
                          "ERROR: program design error");
              spectrum(i) += series(j);
            }
            TFXX_debug(opt.debug, "main", "added result");
          }

          ++isegment;
          TFXX_debug(opt.debug, "main",
                     "finished segment #" << isegment);
        } // while (isegment<opt.nsegments)

        if (opt.nsegments>1) { spectrum /= opt.nsegments; }
        // use reference semantics to make results available under previous
        // name
        Tseries series=spectrum;

  /*----------------------------------------------------------------------*/
  // end of analysis section
  // begin of output section
  /*----------------------------------------------------------------------*/
thomas.forbriger's avatar
thomas.forbriger committed
852

thomas.forbriger's avatar
thomas.forbriger committed
853
        // adjust sampling interval to frequency interval
854
855
        if (opt.amplitudespectrum || opt.powerspectrum || opt.harmonicsignal)
        { wid2.dt=df; }
thomas.forbriger's avatar
thomas.forbriger committed
856
857
        wid2.nsamples=series.size();

thomas.forbriger's avatar
thomas.forbriger committed
858
        os << wid2;
thomas.forbriger's avatar
thomas.forbriger committed
859
        TFXX_debug(opt.debug, "main", "  series and WID are written");
thomas.forbriger's avatar
thomas.forbriger committed
860
861
862
863
864
865
        if (is.hasinfo()) { sff::INFO info; is >> info; os << info; }
        if (is.hasfree() || true) 
        {
          sff::FREE tracefree;
          is >> tracefree;
          tracefree.append(FOUTRA_VERSION);
thomas.forbriger's avatar
thomas.forbriger committed
866
867
          tracefree.append("data read from file " + infile->name);
          tracefree.append(processfree);
thomas.forbriger's avatar
thomas.forbriger committed
868
869
          os << tracefree;
        }
870
        os << series;
thomas.forbriger's avatar
thomas.forbriger committed
871
872
873
874

        // write ASCII output if requested
        if (opt.asciiout)
        {
thomas.forbriger's avatar
thomas.forbriger committed
875
876
          TFXX_debug(opt.debug, "main", "write to ASCII file");
          ++otrace;
thomas.forbriger's avatar
thomas.forbriger committed
877
878
879
880
881
882
883
          std::ostringstream outname;
          outname << opt.asciibase << ".";
          outname.width(3);
          outname.fill('0');
          outname << otrace << ".asc";
          std::string filename(outname.str());
          std::ofstream aos(filename.c_str());
thomas.forbriger's avatar
thomas.forbriger committed
884
885
886

          Tseries f;
          Tseries p;
thomas.forbriger's avatar
thomas.forbriger committed
887
          if (opt.logascii)
thomas.forbriger's avatar
thomas.forbriger committed
888
          {
thomas.forbriger's avatar
thomas.forbriger committed
889
890
            TFXX_debug(opt.debug, "main", 
                       "prepare data with logarithmic scale");
thomas.forbriger's avatar
thomas.forbriger committed
891
892
893
894
895
896
            // f_k = f_min * pow(10,k*ndecades)
            // k = log10( f_k/f_min ) / ndecades
            int kmin=0;
            // f_min = df
            // f_max = df*series.l();
            int kmax=int(std::log10(double(series.l()))/opt.asciidecades);
thomas.forbriger's avatar
thomas.forbriger committed
897
898
899
900
            TFXX_debug(opt.debug, "main", 
                       "  output index range: " << kmin << "-" << kmax);
            f=Tseries(kmin,kmax);
            p=Tseries(kmin,kmax);
thomas.forbriger's avatar
thomas.forbriger committed
901
            double bwfactor=std::sqrt(std::pow(10.,opt.decades));
thomas.forbriger's avatar
thomas.forbriger committed
902
            for (int k=f.f(); k<=f.l(); ++k)
thomas.forbriger's avatar
thomas.forbriger committed
903
904
            {
              // center frequency
thomas.forbriger's avatar
thomas.forbriger committed
905
              double fm=df*std::pow(10.,double(k*opt.asciidecades));
thomas.forbriger's avatar
thomas.forbriger committed
906
              int l=int(0.5+fm/df);
thomas.forbriger's avatar
thomas.forbriger committed
907
908
909
              f(k)=l*df;
              if (opt.avgasciionly)
              {
thomas.forbriger's avatar
thomas.forbriger committed
910
911
912
913
914
915
                // center frequency
                double fm=f(k);
                double f1=fm/bwfactor;
                double f2=fm*bwfactor;
                int l1=int(f1/df);
                int l2=int(f2/df)+1;
thomas.forbriger's avatar
thomas.forbriger committed
916
917
                l1 = l1 < series.f() ? series.f() : l1;
                l2 = l2 > series.l() ? series.l() : l2;
thomas.forbriger's avatar
thomas.forbriger committed
918
919
920
921
922
923
924
925
926
                /*
                TFXX_debug(opt.debug, "main", 
                           "  center frequency: " << fm << " Hz");
                TFXX_debug(opt.debug, "main", 
                           "  average: " << f1 << " Hz - " << f2 << " Hz");
                TFXX_debug(opt.debug, "main", 
                           "  index range: " << l1 << " - " << l2);
                */
                TFXX_assert(l2>=l1, "negative bandwidth for averaging");
thomas.forbriger's avatar
thomas.forbriger committed
927
928
                Tseries sub=aff::subarray(series)(l1,l2);
                p(k)=aff::func::avg(sub);
929
              } // if (opt.avgasciionly)
thomas.forbriger's avatar
thomas.forbriger committed
930
931
932
              else
              {
                p(k)=series(l);
933
934
935
              } // if (opt.avgasciionly) else
            } // for (int k=f.f(); k<=f.l(); ++k)
          } // if (opt.logascii)
thomas.forbriger's avatar
thomas.forbriger committed
936
937
          else
          {
thomas.forbriger's avatar
thomas.forbriger committed
938
939
940
941
942
            // just copy and prepare frequency vector
            f=Tseries(series.size());
            p=Tseries(series.size());
            p.copyin(series);
            for (int k=f.f(); k<=f.l(); ++k)
thomas.forbriger's avatar
thomas.forbriger committed
943
            {
thomas.forbriger's avatar
thomas.forbriger committed
944
              f(k) = k*df;
thomas.forbriger's avatar
thomas.forbriger committed
945
            }
946
          } // if (opt.logascii) else
thomas.forbriger's avatar
thomas.forbriger committed
947
          for (int k=f.f(); k<=f.l(); ++k)
thomas.forbriger's avatar
thomas.forbriger committed
948
949
950
951
952
953
954
955
956
          {
            // center frequency
            aos.setf(std::ios_base::scientific,std::ios_base::floatfield);
            aos.precision(10);
            aos << f(k) << " ";
            aos.setf(std::ios_base::scientific,std::ios_base::floatfield);
            aos.precision(10);
            aos << p(k) << endl;
          }
957
        } // if (opt.asciiout)
thomas.forbriger's avatar
thomas.forbriger committed
958

thomas.forbriger's avatar
thomas.forbriger committed
959
960
        TFXX_debug(opt.debug, "main", 
                   "trace #" << itrace << " successfully processed");
961
      } // if ((!doselect) || traceranges.contains(itrace))
thomas.forbriger's avatar
thomas.forbriger committed
962
963
964
965
966
967
      else
      {
        TFXX_debug(opt.debug, "main", "skip trace #" << itrace );
        if (opt.verbose)
        { std::cout << "  skip trace #" << itrace << std::endl; }
        is.skipseries();
968
      } // if ((!doselect) || traceranges.contains(itrace)) else
969
    } // while (is.good())
thomas.forbriger's avatar
thomas.forbriger committed
970
971
972
973
    
    // go to next file
    firstfile=false;
    ++infile;
974
  } // while (infile != arguments.end())
thomas.forbriger's avatar
thomas.forbriger committed
975

thomas.forbriger's avatar
thomas.forbriger committed
976
} // main
thomas.forbriger's avatar
thomas.forbriger committed
977

thomas.forbriger's avatar
thomas.forbriger committed
978
979
980
981
982
983
984
985
986
987
988
/*======================================================================*/

/*! \page page_foutra Spectral analysis (foutra.cc)
 *
 * Sections in this page:
 *   - \ref sec_foutra_fourier_scaling
 *   - \ref sec_foutra_hanning_PSD_scaling
 *   - \ref sec_foutra_scaling_harmonics
 *     - \ref subsec_foutra_scaling_harmonics_harmsig
 *     - \ref subsec_foutra_scaling_harmonics_boxcar
 *     - \ref subsec_foutra_scaling_harmonics_hanning
989
 *     - \ref subsec_foutra_scaling_harmonics_tests
thomas.forbriger's avatar
thomas.forbriger committed
990
991
992
 *
 * <HR>
 * \section sec_foutra_fourier_scaling Scaling of Fourier transforms
993
994
995
996
 *
 * The Fourier transform presented in libfourier is scaled appropriately to
 * represent coefficients of the Fourier integral transform (see libfourier
 * documentation).
thomas.forbriger's avatar
thomas.forbriger committed
997
 *
thomas.forbriger's avatar
thomas.forbriger committed
998
999
1000
1001
1002
 * The maximum of the transform of the Boxcar taper is T,
 * where T is the duration of the window. Similarly the maximum
 * of the Fourier transform of the Hanning taper is 0.5*T, where the duration of
 * the window is T.
 *
thomas.forbriger's avatar
thomas.forbriger committed
1003
1004
1005
1006
1007
1008
1009
1010
 * \date 12/2010 \author thof
 */

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

/*! \page page_foutra
 * <HR>
 * \section sec_foutra_hanning_PSD_scaling Scaling factor for Hanning taper when applied in PSD calculation
thomas.forbriger's avatar
thomas.forbriger committed
1011
1012
 *
 * Tapering stationary noise with a Hanning taper reduces total signal energy
1013
 * an thus the signal power obtained through an FFT. In his tutorial Agnew
thomas.forbriger's avatar
thomas.forbriger committed
1014
1015
 * recommends to scale each sample \f$s_k\f$ of the series by \f$w_k/W\f$,
 * where \f$w_k\f$ is the \f$k\f$-th sample of the taper and 
1016
 *
thomas.forbriger's avatar
thomas.forbriger committed
1017
1018
1019
 * \f[
 *   W^2 = \frac{1}{N} \sum\limits_0^{N-1} w_k^2
 * \f]
1020
1021
1022
1023
1024
 *
 * is a measure for the loss in total signal power due to application of the
 * taper. This applies only for staionary noise.
 * 
 * From Walter's Matlab
thomas.forbriger's avatar
thomas.forbriger committed
1025
1026
1027
1028
1029
 * scripts I took:
 *
 * If data is the time series vector and y is the Hanning taper of appropriate
 * length, Walter calculates
 *   
thomas.forbriger's avatar
thomas.forbriger committed
1030
 * \code
thomas.forbriger's avatar
thomas.forbriger committed
1031
 *   w = sqrt(sum(y.*y)/ndata);
thomas.forbriger's avatar
thomas.forbriger committed
1032
 * \endcode
thomas.forbriger's avatar
thomas.forbriger committed
1033
1034
1035
 *
 * and
 *   
thomas.forbriger's avatar
thomas.forbriger committed
1036
 * \code
thomas.forbriger's avatar
thomas.forbriger committed
1037
 *   y=y/w;
thomas.forbriger's avatar
thomas.forbriger committed
1038
 * \endcode
thomas.forbriger's avatar
thomas.forbriger committed
1039
1040
1041
 *
 * where ndata is the length of data and y.
 *
thomas.forbriger's avatar
thomas.forbriger committed
1042
1043
1044
1045
 * For a Hanning taper
 * \f[
 *   w_k = \sin^2(k \pi/(N-1)).
 * \f]
thomas.forbriger's avatar
thomas.forbriger committed
1046
 *
1047
 * Thus
thomas.forbriger's avatar
thomas.forbriger committed
1048
 *
thomas.forbriger's avatar
thomas.forbriger committed
1049
1050
1051
 * \f[
 *   W^2 = \frac{1}{N} \sum\limits_0^{N-1} \sin^4(k\pi/(N-1))
 * \f]
thomas.forbriger's avatar
thomas.forbriger committed
1052
 *
1053
 * From Gradshteyn and Ryzhik (eq. 1.321 3) I find
thomas.forbriger's avatar
thomas.forbriger committed
1054
 *
thomas.forbriger's avatar
thomas.forbriger committed
1055
1056
 * \f[
 *   \sin^4(k\pi/(N-1))
thomas.forbriger's avatar
thomas.forbriger committed
1057
 *
thomas.forbriger's avatar
thomas.forbriger committed
1058
1059
1060
1061
 *     = \frac{1}{8} \left( \cos(4 k\pi/(N-1))
 *                 - 4 \cos(2 k\pi/(N-1))
 *                 + 3 \right).
 * \f]
thomas.forbriger's avatar
thomas.forbriger committed
1062
 *
1063
1064
 * Within the sum the contribution of both cos-terms will vanish, since both
 * are averaged over one and two periods, respectively. Thus
thomas.forbriger's avatar
thomas.forbriger committed
1065
 *
thomas.forbriger's avatar
thomas.forbriger committed
1066
1067
1068
 * \f[
 *   W^2 = \frac{1}{N} \sum\limits_0^{N-1} \frac{3}{8} = \frac{3}{8}.
 * \f]
thomas.forbriger's avatar
thomas.forbriger committed
1069
 *
1070
1071
 * Since foutra is not scaling the taper but scaling the power spectrum, we
 * have to apply the factor 8/3 to the result of power spectrum calculation.
thomas.forbriger's avatar
thomas.forbriger committed
1072
 *
thomas.forbriger's avatar
thomas.forbriger committed
1073
 * This factor 8/3=2.66667 was tested against the value for \f$W^2\f$, when
1074
 * explicitely derived from a Hanning taper time series by the above formula.
thomas.forbriger's avatar
thomas.forbriger committed
1075
 *
thomas.forbriger's avatar
thomas.forbriger committed
1076
1077
1078
1079
1080
1081
1082
 * In the Makefile you will find a section with test code. 
 * This offers instantaneous testing of PSD scaling in foutra.
 * \dontinclude Makefile
 * \skip # test foutra PSD scaling
 * \until #----------------------------------------------------------------------
 *
 * \date 13/9/2007 \author thof
1083
 *
thomas.forbriger's avatar
thomas.forbriger committed
1084
1085
1086
1087
1088
1089
1090
 */

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

/*! \page page_foutra
 * <HR>
 * \section sec_foutra_scaling_harmonics Scaling for harmonic signals
1091
1092
1093
 *
 * For harmonic signals the FFT normalized to the duration of the available
 * time window has an amplitude of the value of the maximum of the Fourier
thomas.forbriger's avatar
thomas.forbriger committed
1094
 * transform of the window function times half the amplitude of the time domain
1095
 * signal. To obtain a spectral representation with peaks of the amplitude of
thomas.forbriger's avatar
thomas.forbriger committed
1096
 * the time domain signal, the FFT must be scaled accordingly.
1097
 *
thomas.forbriger's avatar
thomas.forbriger committed
1098
1099
1100
 * We understand 
 * \f[
 *   \tilde{f}(\omega)=\int\limits_{-\infty}^{+\infty} 
thomas.forbriger's avatar
thomas.forbriger committed
1101
 *     f(t) \,e^{-i\omega t} \textrm{d}{t}.
thomas.forbriger's avatar
thomas.forbriger committed
1102
1103
1104
1105
 * \f]
 * as the Fourier transformation of the time domain signal
 * \f[
 *   f(t)=\int\limits_{-\infty}^{+\infty} 
thomas.forbriger's avatar
thomas.forbriger committed
1106
 *     \tilde{f}(\omega) \,e^{i\omega t} \frac{\textrm{d}{\omega}}{2\pi}.
thomas.forbriger's avatar
thomas.forbriger committed
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
 * \f]
 * If we apply a time domain window function \f$w(t)\f$ to the function
 * \f$f(t)\f$, we obtain the tapered function
 * \f[
 *   g(t)=f(t) w(t)
 * \f]
 * and its Fourier transform
 * \f[
 *   \tilde{g}(\omega)=\int\limits_{-\infty}^{+\infty}
 *     \tilde{f}(\omega')\,\tilde{w}(\omega-\omega')\,\textrm{d}\omega.
 * \f]
1118
 *
thomas.forbriger's avatar
thomas.forbriger committed
1119
 * \subsection subsec_foutra_scaling_harmonics_harmsig Application to harmonic signals 
1120
 *
thomas.forbriger's avatar
thomas.forbriger committed
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
 * Let
 * \f[
 *   f(t)=A\,\cos(\omega_0 t+\phi)
 * \f]
 * be the harmonic signal under investigation.
 * Then
 * \f[
 *   f(t)=A\left\{\cos(\omega_0 t)\cos(\phi)-
 *     \sin(\omega_0 t)\sin(\phi)\right\}
 * \f]
 * and
 * \f[
 *   \tilde{f}(\omega) =
 *     \frac{A}{2}\left\{
 *       \cos(\phi)
 *         \left[
 *           \delta(\omega-\omega_0) 
 *           +
 *           \delta(\omega+\omega_0) 
 *         \right]
 *       +i\sin(\phi)
 *         \left[
 *           \delta(\omega-\omega_0) 
 *           -
 *           \delta(\omega+\omega_0) 
 *         \right]
 *     \right\},
 * \f]
 * where \f$\delta(\omega)\f$ is Dirac's delta function with
 * \f[
 *   \delta(\omega)=\left\{
 *     \begin{array}{ll}
 *       \infty & \textrm{if $\omega=0$ and}\\
 *       0 & \textrm{otherwise}
 *     \end{array}
 *   \right.
 * \f]
 * and
 * \f[
 *   \int\limits_{-\infty}^{+\infty}
1161
 *     \delta(\omega)\,\textrm{d}\omega=1
thomas.forbriger's avatar
thomas.forbriger committed
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
 * \f]
 * such that
 * \f[
 *   \tilde{f}(\omega)=\int\limits_{-\infty}^{+\infty}
 *     \tilde{f}(\omega')\,\delta(\omega-\omega')\,\textrm{d}\omega.
 * \f]
 * This way I obtain
 * \f[
 *   \tilde{g}(\omega)=
 *     \frac{A}{2}\left\{
 *       e^{i\phi}\tilde{w}(\omega-\omega_0)
 *       +
 *       e^{-i\phi}\tilde{w}(\omega+\omega_0)
 *     \right\}
 * \f]
 * for the Fourier transform of the tapered function.
 * If we ignore interference with side-lobes from the negative frequency
 * \f$-\omega_0\f$ and side-lobes of potential other harmonics at nearby
 * frequencies, we can approximate
 * \f[
 *   \tilde{g}(\omega_0)\approx\frac{A}{2}e^{i\phi}\tilde{w}(0)
 * \f]
 * and
 * \f[
 *   \left|\tilde{g}(\omega_0)\right|\approx
 *     \frac{A}{2}\left|\tilde{w}(0)\right|.
 * \f]
1189
 *
thomas.forbriger's avatar
thomas.forbriger committed
1190
 * \subsection subsec_foutra_scaling_harmonics_boxcar Boxcar taper function
1191
 *
thomas.forbriger's avatar
thomas.forbriger committed
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
 * The boxcar taper is defined as
 * \f[
 *   w(t) = 
 *   \left\{
 *     \begin{array}{ll}
 *       1 & \textrm{if $|t|\leq T/2$ and}\\
 *       0 & \textrm{otherwise}
 *     \end{array}
 *   \right.
 * \f]
 * with
 * \f[
 *   \tilde{w}(\omega)
  *    = T \frac{\sin(\omega T/2)}{\omega T/2}
 * \f]
 * and
 * \f[
 *   w_{\textrm{max}}=w(0)=T.
 * \f]
1211
 *
thomas.forbriger's avatar
thomas.forbriger committed
1212
 * \subsection subsec_foutra_scaling_harmonics_hanning Hanning taper function
1213
 *
1214
 * The Hanning taper is defined as
thomas.forbriger's avatar
thomas.forbriger committed
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
 * \f[
 *   w(t) = 
 *   \left\{
 *     \begin{array}{ll}
 *       \cos^2(\pi t / T)) 
 *       = \frac{1}{2}\left[
 *           1 + \cos(2\pi t / T)
 *         \right]
 *       & \textrm{if $|t|\leq T/2$ and}\\
 *       0 & \textrm{otherwise}
 *     \end{array}
 *   \right.
 * \f]
 * with
 * \f[
 *   \tilde{w}(\omega)
 *    = T/2 \frac{\sin(\omega T/2)}{\omega T/2}
 *    + T/4 \frac{\sin(\omega T/2+\pi)}{\omega T/2+\pi}
 *    + T/4 \frac{\sin(\omega T/2-\pi)}{\omega T/2-\pi}
 * \f]
 * (see Blackman, R.B. and Tukey, J.W. 1958.
 *      The measurement of power spectra.
 *      Dover Publications, Inc., New York.
 *      Section B.5)
 * and
 * \f[
 *   w_{\textrm{max}}=w(0)=\frac{T}{2}.
 * \f]
1243
 *
1244
1245
1246
1247
1248
1249
1250
1251
 * \subsection subsec_foutra_scaling_harmonics_tests Instant tests
 *
 * In the Makefile you will find a section with test code. 
 * This offers instantaneous testing of harmonic signal scaling in foutra.
 * \dontinclude Makefile
 * \skip # test foutra scaling for harmonic signals
 * \until #----------------------------------------------------------------------
 *
thomas.forbriger's avatar
thomas.forbriger committed
1252
 * \date 10.01.2011 \author thof
thomas.forbriger's avatar
thomas.forbriger committed
1253
1254
 */

thomas.forbriger's avatar
thomas.forbriger committed
1255
/* ----- END OF foutra.cc ----- */