mtupspec.cpp

Go to the documentation of this file.

#include <iostream>
#include <string>
#include <fstream>
#include <vector>
#include "StackedSpectrum.h"
#include "WFunc.h"
#include "TimeSeriesData.h"
#include "Util.h"

using namespace std;
using namespace gplib;

string version = "$Id: mtupspec.cpp 1816 2009-09-07 11:28:35Z mmoorkamp $";

/*!
 * \addtogroup UtilProgs Utility Programs
 *@{
 * \file
 * Calculate the power spectra for the 4 horizontal components in a MT time-series.
 * Output is in ascii format for plotting with xmgrace etc.
 */

//! Calculate the power spectrum and write it to an ascii file
void CalcPSpecAndWrite(ttsdata Data, const double samplerate,
    const int seglength, const string name)
  {
    vector<complex<double> > Spectrum(seglength / 2 + 1);
    StackedSpectrum(Data.begin(), Data.end(), Spectrum.begin(), seglength,
        Hanning());
    ofstream outfile(name.c_str());
    for (unsigned int i = 1; i < Spectrum.size(); ++i) //we do not output the static contribution (0 frequency)
      outfile << i * samplerate / seglength << " " << abs(Spectrum.at(i))
          << endl;
  }

int main(int argc, char *argv[])
  {
    TimeSeriesData TsData;

    string infilename;
    int seglength = 2400;
    cout
        << "This is mtupspec: Calculate power spectra from  Phoenix time series"
        << endl << endl;
    cout
        << " Usage      mtupspec filename     if no filename given, the program will ask for one"
        << endl;
    cout << " Output will be 4 ascii file with ending _specex, _specey etc."
        << endl << endl;
    cout << " This is Version: " << version << endl << endl;

    //get the filename information
    if (argc == 2)
      {
        infilename = argv[1];
      }
    else
      {
        infilename = AskFilename("Infilename: ");
      }

    //get the length of th individual segments for the fouries transform
    cout << "Length for each segment: ";
    cin >> seglength;
    // we need at least 2 points in each segment
    if (seglength < 2)
      {
        std::cerr << "Segment must be longer than 2 !";
        return 100;
      }
    // the segment length cannot be longer than the time series
    if (seglength > TsData.GetData().GetEx().GetData().size())
      {
        std::cerr << "Segment must shorter than the time series !";
        return 200;
      }
    TsData.GetData(infilename);

    const double samplerate = TsData.GetData().GetSamplerate();
    //calculate the power spectra for each component and write them to a file
    CalcPSpecAndWrite(TsData.GetData().GetEx().GetData(), samplerate,
        seglength, (infilename + "_specex"));
    CalcPSpecAndWrite(TsData.GetData().GetEy().GetData(), samplerate,
        seglength, (infilename + "_specey"));
    CalcPSpecAndWrite(TsData.GetData().GetHx().GetData(), samplerate,
        seglength, (infilename + "_spechx"));
    CalcPSpecAndWrite(TsData.GetData().GetHy().GetData(), samplerate,
        seglength, (infilename + "_spechy"));
    CalcPSpecAndWrite(TsData.GetData().GetHz().GetData(), samplerate,
        seglength, (infilename + "_spechz"));

  }
/*@}*/