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}.