foutra.cc

Go to the documentation of this file.

#define FOUTRA_VERSION \
  "FOUTRA   V1.10   Fourier transforms"

#include <iostream>
#include <tfxx/commandline.h>
#include <aff/series.h>
#include <aff/iterator.h>
#include <aff/dump.h>
#include <aff/seriesoperators.h>
#include <aff/functions/avg.h>
#include <aff/functions/rms.h>
#include <aff/functions/sqrsum.h>
#include <aff/subarray.h>
#include <tsxx/tsxx.h>
#include <tsxx/tapers.h>
#include <tsxx/filter.h>
#include <tsxx/wid2timeseries.h>
#include <tsioxx/tsioxx.h>
#include <fstream>
#include <tfxx/error.h>
#include <tfxx/rangestring.h>
#include <tfxx/xcmdline.h>
#include <tfxx/misc.h>
#include <tfxx/handle.h>
#include <tfxx/seitosh.h>
#include <datrwxx/readany.h>
#include <datrwxx/writeany.h>
#include <fourier/fftwaff.h>
#include <sstream>

#include "foutra_options_brief_help.h"
#include "foutra_options_help.h"
#include "foutra_usage_help.h"

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

struct Options {
  bool verbose, overwrite, debug, asciiout, logascii;
  std::string inputformat, asciibase, outputformat;
  bool amplitudespectrum, powerspectrum, boxcartaper;
  bool avgconstbw, avgrelbw, avgasciionly;
  bool demean, detrend, scalerbw, adjustdivisor;
  bool derivative;
  int ndemean, ndetrend, divisor;
  double decades, scaledecades, asciidecades, nderivative;
  int avgsamples;
  bool reportrms, harmonicsignal, padzeroes;
  int padfactor, nsegments;
}; // 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;

typedef ts::TDsfftimeseries Ttimeseries;
typedef Ttimeseries::Tseries Tseries;

typedef fourier::fft::DRFFTWAFF Tfft;

/*======================================================================*/

int main(int iargc, char* argv[])
{

  // define usage information
  char version_text[]=
  {
    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,"-"},
    // 2: overwrite mode
    {"o",arg_no,"-"},
    // 3: input format
    {"type",arg_yes,"sff"},
    // 4: debug mode
    {"D",arg_no,"-"},
    // 5: amplitude spectrum
    {"amplitude",arg_no,"-"},
    // 6: power spectrum
    {"power",arg_no,"-"},
    // 7: apply boxcar taper
    {"boxcar",arg_no,"-"},
    // 8: average over constant number of samples
    {"avg",arg_opt,"21"},
    // 9: average over constant relative bandwidth (in decades)
    {"rbw",arg_opt,"0.167"},
    // 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"},
    // 14: write result to ASCII files
    {"logascii",arg_opt,"0.167"},
    // 15: average ASCII output only
    {"avgascii",arg_no,"-"},
    // 16: divisor for number of samples
    {"divisor",arg_opt,"100"},
    // 17: report rms value
    {"rms",arg_no,"-"},
    // 18: report rms value
    {"harmonic",arg_no,"-"},
    // 19: pad time series
    {"pad",arg_yes,"1"},
    // 20: subdivide input series
    {"nsegments",arg_yes,"1"},
    // 21: extended help
    {"xhelp",arg_no,"-"},
    // 22: input format
    {"Type",arg_yes,"sff"},
    // 23: input format
    {"derivative",arg_opt,"1"},
    {NULL}
  };

  static const char tracekey[]="t";

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

  // no arguments? print usage...
  if (iargc<2) 
  {
    cerr << version_text << endl;
    cerr << foutra_options_brief_help << endl;
    cerr << tfxx::seitosh::repository_reference << endl;
    exit(0);
  }

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

  // help requested? print full help text...
  if (cmdline.optset(0) || cmdline.optset(21))
  {
    cerr << version_text << endl;
    cerr << foutra_options_brief_help << endl;
    cerr << foutra_options_help << endl;
    cerr << foutra_usage_help << endl;
    datrw::supported_data_types(cerr);
    if (cmdline.optset(21)) { datrw::online_help(cerr); }
    cerr << endl << tfxx::seitosh::repository_reference << endl;
    exit(0);
  }

  Options opt;
  opt.verbose=cmdline.optset(1);
  opt.overwrite=cmdline.optset(2);
  opt.inputformat=cmdline.string_arg(3);
  opt.debug=cmdline.optset(4);
  opt.amplitudespectrum=cmdline.optset(5);
  if (!opt.amplitudespectrum)
  { opt.powerspectrum=cmdline.optset(6); }
  opt.boxcartaper=cmdline.optset(7);
  opt.avgconstbw=cmdline.optset(8);
  opt.avgsamples=cmdline.int_arg(8);
  opt.avgrelbw=cmdline.optset(9);
  opt.decades=cmdline.double_arg(9);
  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);
  opt.logascii=cmdline.optset(14);
  opt.asciidecades=cmdline.double_arg(14);
  opt.avgasciionly=cmdline.optset(15);
  opt.adjustdivisor=cmdline.optset(16);
  opt.divisor=cmdline.int_arg(16);
  opt.reportrms=cmdline.optset(17);
  opt.harmonicsignal=cmdline.optset(18);
  opt.padzeroes=cmdline.optset(19);
  opt.padfactor=cmdline.int_arg(19);
  opt.nsegments=cmdline.int_arg(20);
  opt.outputformat=cmdline.string_arg(22);
  opt.derivative=cmdline.optset(23);
  opt.nderivative=cmdline.double_arg(23);

  TFXX_assert((opt.divisor > 0), "illegal value for argument divisor");
  TFXX_assert((opt.nsegments > 0), "illegal value for argument nsegments");
  TFXX_assert((opt.padfactor > 0), "illegal value for argument pad");

  if (opt.verbose)
  { cout << FOUTRA_VERSION << endl; }

  // 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)
  {
    cout << "NOTICE: file specific information (SRCE line and file FREE)" <<
      endl
      <<    "        will be taken from first file only!" << endl;
  }

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

  /*======================================================================*/
  // create processing description
  sff::FREE processfree;
  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");
  }
  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());
  }
  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");
  }
  if (opt.padzeroes) 
  {
    std::ostringstream freeline;
    freeline << "pad with zeroes increasing the total "
      << "number of samples by a factor of " << opt.padfactor;
    processfree.append(freeline.str());
  }
  if (opt.harmonicsignal || opt.amplitudespectrum || opt.powerspectrum)
  {
    processfree.append("An appropriately scaled FFT (libdrfftw) is applied");
    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());
    }
  }
  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");
    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");
    }
    if ((opt.avgconstbw || opt.avgrelbw) && !opt.avgasciionly)
    {
      processfree.append("The spectrum is smoothed with a boxcar over");
      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());
      }
    }
    else
    {
      processfree.append("The spectrum is not smoothed");
      if (opt.avgasciionly && opt.asciiout)
      {
        processfree.append("The ASCII output is smoothed with a boxcar over");
        std::ostringstream freeline;
        freeline << opt.decades << " decades";
        processfree.append(freeline.str());
      }
    }
  }
  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");
  }
  else
  {
    processfree.append("The output is a tapered version of the time series");
  }
  if (opt.amplitudespectrum || opt.powerspectrum || opt.harmonicsignal)
  {
    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)
  {
    cout << "processing of each trace will take place as follows:" << endl;
    for(sff::FREE::Tlines::const_iterator I=processfree.lines.begin(); 
        I != processfree.lines.end(); ++I) 
    { cout << "  " << *I << std::endl; }
  }

  /*======================================================================*/
  // start processing
  
  // create taper
  ts::tapers::Hanning taper;
  // create FFTW processor
  Tfft fft;

  // 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());
  datrw::oanystream os(ofs, opt.outputformat);

  // prepare file FREE block
  sff::FREE filefree;
  filefree.append(FOUTRA_VERSION);
  // set flag to process header of first input file
  bool firstfile=true;
  // count output traces
  int otrace=0;
  // 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());
    datrw::ianystream is(ifs, opt.inputformat);
    // 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; }
        Tseries inputseries;
        is >> inputseries;
        sff::WID2 wid2;
        is >> wid2;
        TFXX_debug(opt.debug, "main", "  series and WID2 are read");

  /*----------------------------------------------------------------------*/
  // report rms value if requested
        if (opt.reportrms)
        {
          std::cout << "    rms-value of input time series: "
            << aff::func::rms(inputseries) << std::endl;
        }
  /*----------------------------------------------------------------------*/
  // process data
  //
  // 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); 
        }

        // adjust length of time series
        if (opt.adjustdivisor)
        {
          if (opt.verbose)
          {
            std::cout << "    adjust divisor to " << opt.divisor << std::endl;
            std::cout << "    number of input samples:      " << wid2.nsamples
              <<std::endl;
            std::cout << "    number of segment samples:    " << nsegsamples
              <<std::endl;
          }
          int rest=nsegsamples % opt.divisor;
          nsegsamples=nsegsamples-rest;
          if (opt.verbose)
          {
            std::cout << "    number of processing samples: " << nsegsamples
              <<std::endl;
          }
        }

        // segment stride
        int segstride=inputseries.size()/(opt.nsegments+1);

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

        // length of time series segment
        double T=wid2.dt*nsegsamples;
        // length of taper window
        double Tw=T;
        // frequency sampling interval
        double df=1./T;

        if (opt.padzeroes)
        {
          T *= opt.padfactor;
          df /= opt.padfactor;
        }

        if (opt.verbose)
        {
          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;
        }
            
        // the result will be collected in
        Tseries spectrum;

        // segment counter
        int isegment=0;

        // start actual processing
        // -----------------------
        while (isegment<opt.nsegments)
        {
          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));
          }

          // detrend
          if (opt.detrend) { 
            TFXX_debug(opt.debug, "main", "  detrend");
            ts::filter::RemoveTrend detrend(opt.ndetrend);
            detrend(ts::filter::Ttimeseries(series, wid2.dt));
          }

          // apply taper
          if (!opt.boxcartaper) { 
            TFXX_debug(opt.debug, "main", "  apply taper");
            taper.apply(series); 
          }

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

          // call FFT for amplitude or power spectrum or harmonic signal
          if (opt.amplitudespectrum || opt.powerspectrum || opt.harmonicsignal)
          {
            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);

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

            TFXX_debug(opt.debug, "main", "calculate square of modulus and copy");
            while(S.valid() && C.valid())
            {
              *S = C->real()*C->real()+C->imag()*C->imag();
              ++S; ++C;
            }
          }
            
          // take sqrt for amplitude spectrum
          if (opt.amplitudespectrum || opt.harmonicsignal)
          {
            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;
          }

          // apply appropriate scaling and averaging for power spectrum
          if (opt.powerspectrum)
          {
            // scaling factor to adjust for taper effect
            double tapscaling=1.;
            if (opt.boxcartaper)
            {
              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
            //
            // This factor computes the one-sided PSD, i.e. twice the Fourier
            // transform of the normalized auto-correlation function.
            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");
              double bwfactor=(std::pow(10.,opt.scaledecades)-1.)/
                std::pow(10.,opt.scaledecades/2.);
              // cout << "bwfactor: " << bwfactor << endl;
              // cout << "index range: " 
              //   << series.f() << "-" << series.l() << endl;
              for (int k=series.f(); k<=series.l(); ++k)
              {
                // center frequency
                double fm=k*df;
                series(k) *= scalingfactor*bwfactor*fm;
              }
            } // 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)
          {
            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)
            {
              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
  /*----------------------------------------------------------------------*/

        // adjust sampling interval to frequency interval
        if (opt.amplitudespectrum || opt.powerspectrum || opt.harmonicsignal)
        { wid2.dt=df; }
        wid2.nsamples=series.size();

        os << wid2;
        TFXX_debug(opt.debug, "main", "  series and WID are written");
        if (is.hasinfo()) { sff::INFO info; is >> info; os << info; }
        if (is.hasfree() || true) 
        {
          sff::FREE tracefree;
          is >> tracefree;
          tracefree.append(FOUTRA_VERSION);
          tracefree.append("data read from file " + infile->name);
          tracefree.append(processfree);
          os << tracefree;
        }
        os << series;

        // write ASCII output if requested
        if (opt.asciiout)
        {
          TFXX_debug(opt.debug, "main", "write to ASCII file");
          ++otrace;
          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());

          Tseries f;
          Tseries p;
          if (opt.logascii)
          {
            TFXX_debug(opt.debug, "main", 
                       "prepare data with logarithmic scale");
            // 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);
            TFXX_debug(opt.debug, "main", 
                       "  output index range: " << kmin << "-" << kmax);
            f=Tseries(kmin,kmax);
            p=Tseries(kmin,kmax);
            double bwfactor=std::sqrt(std::pow(10.,opt.decades));
            for (int k=f.f(); k<=f.l(); ++k)
            {
              // center frequency
              double fm=df*std::pow(10.,double(k*opt.asciidecades));
              int l=int(0.5+fm/df);
              f(k)=l*df;
              if (opt.avgasciionly)
              {
                // 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;
                l1 = l1 < series.f() ? series.f() : l1;
                l2 = l2 > series.l() ? series.l() : l2;
                /*
                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");
                Tseries sub=aff::subarray(series)(l1,l2);
                p(k)=aff::func::avg(sub);
              } // if (opt.avgasciionly)
              else
              {
                p(k)=series(l);
              } // if (opt.avgasciionly) else
            } // for (int k=f.f(); k<=f.l(); ++k)
          } // if (opt.logascii)
          else
          {
            // 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)
            {
              f(k) = k*df;
            }
          } // if (opt.logascii) else
          for (int k=f.f(); k<=f.l(); ++k)
          {
            // 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;
          }
        } // if (opt.asciiout)

        TFXX_debug(opt.debug, "main", 
                   "trace #" << itrace << " successfully processed");
      } // if ((!doselect) || traceranges.contains(itrace))
      else
      {
        TFXX_debug(opt.debug, "main", "skip trace #" << itrace );
        if (opt.verbose)
        { std::cout << "  skip trace #" << itrace << std::endl; }
        is.skipseries();
      } // if ((!doselect) || traceranges.contains(itrace)) else
    } // while (is.good())
    
    // go to next file
    firstfile=false;
    ++infile;
  } // while (infile != arguments.end())

} // main

/*======================================================================*/

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

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

/* ----- END OF foutra.cc ----- */