simple_processing.cpp

Go to the documentation of this file.

#include <iostream>
#include <iterator>
#include <fstream>
#include "TsSpectrum.h"
#include "TimeSeriesData.h"
#include "statutils.h"
#include "MTStation.h"
#include "types.h"
#include "WFunc.h"
#include "TimeFrequency.h"
#include "VecMat.h"
#include "MtuFilter.h"
#include "Util.h"

using namespace std;
using namespace gplib;

string version =
    "$Id: simple_processing.cpp 1816 2009-09-07 11:28:35Z mmoorkamp $";
int main(int argc, char *argv[])
  {
    cout
        << "This is simpleproc: Simple most estimation of MT transfer function, spectral stacking"
        << endl << endl;
    cout << " Usage   simpleproc  inputfilename";
    cout
        << " Output will have the same name as MTU Input with '.mtt' appended "
        << endl << endl;
    cout << " This is Version: " << version << endl << endl;

    TimeSeriesData MtuData;
    MTStation Output;
    string filename, basefilename;
    const unsigned int seglength = 2400;

    if (argc == 2)
      {
        filename = argv[1];
      }
    else
      {
        filename = AskFilename("Time Series Filename: ");
      }
    string datafilename = filename;
    size_t dotpos = filename.find('.', 0);
    if (dotpos != string::npos)
      basefilename = filename.erase(dotpos);
    string exfiltername = basefilename + "_1.cts";
    string eyfiltername = basefilename + "_2.cts";
    string hxfiltername = basefilename + "_3.cts";
    string hyfiltername = basefilename + "_4.cts";
    MtuData.GetData(datafilename.c_str());

    const unsigned int nfreqs = seglength / 2 + 1;
    tcompdata ExHxCorr(nfreqs, 0);
    tcompdata ExHyCorr(nfreqs, 0);
    tcompdata EyHxCorr(nfreqs, 0);
    tcompdata EyHyCorr(nfreqs, 0);
    tcompdata HxHxCorr(nfreqs, 0);
    tcompdata HxHyCorr(nfreqs, 0);
    tcompdata HyHyCorr(nfreqs, 0);
    tcompdata Hdet(nfreqs, 0);
    tcompdata Zxx(nfreqs, 0);
    tcompdata Zxy(nfreqs, 0);
    tcompdata Zyx(nfreqs, 0);
    tcompdata Zyy(nfreqs, 0);

    /*if (MtuData.GetData().GetSamplerate() == 1)
     {
     double rate;
     cout << "Enter sample rate in Hz: ";
     cin >> rate;
     MtuData.GetData().SetSamplerate(rate);
     }*/

    double freqstep = MtuData.GetData().GetSamplerate() / (seglength);
    MtuFilter ExFilter(seglength, freqstep), EyFilter(seglength, freqstep),
        HxFilter(seglength, freqstep), HyFilter(seglength, freqstep);

    ExFilter.GetData(exfiltername);
    EyFilter.GetData(eyfiltername);
    HxFilter.GetData(hxfiltername);
    HyFilter.GetData(hyfiltername);

    gplib::cmat ExTimeFrequency, EyTimeFrequency, HxTimeFrequency,
        HyTimeFrequency;

    ExTimeFrequency = TimeFrequency(
        MtuData.GetData().GetEx().GetData().begin(),
        MtuData.GetData().GetEx().GetData().end(), seglength, Hamming());
    EyTimeFrequency = TimeFrequency(
        MtuData.GetData().GetEy().GetData().begin(),
        MtuData.GetData().GetEy().GetData().end(), seglength, Hamming());
    HxTimeFrequency = TimeFrequency(
        MtuData.GetData().GetHx().GetData().begin(),
        MtuData.GetData().GetHx().GetData().end(), seglength, Hamming());
    HyTimeFrequency = TimeFrequency(
        MtuData.GetData().GetHy().GetData().begin(),
        MtuData.GetData().GetHy().GetData().end(), seglength, Hamming());

    const unsigned int nsegs = ExTimeFrequency.size1();
    ofstream logfile((datafilename + ".log").c_str());
    for (size_t i = 0; i < Zxx.size(); ++i)
      {
        for (size_t j = 0; j < nsegs; ++j)
          {
            ExTimeFrequency(j, i) /= ExFilter.GetFilterCoeff().at(i);
            EyTimeFrequency(j, i) /= EyFilter.GetFilterCoeff().at(i);
            HxTimeFrequency(j, i) /= HxFilter.GetFilterCoeff().at(i);
            HyTimeFrequency(j, i) /= HyFilter.GetFilterCoeff().at(i);
            ExHxCorr.at(i) += (ExTimeFrequency(j, i) * conj(HxTimeFrequency(j,
                i)));
            ExHyCorr.at(i) += (ExTimeFrequency(j, i) * conj(HyTimeFrequency(j,
                i)));
            EyHxCorr.at(i) += (EyTimeFrequency(j, i) * conj(HxTimeFrequency(j,
                i)));
            EyHyCorr.at(i) += (EyTimeFrequency(j, i) * conj(HyTimeFrequency(j,
                i)));
            HxHyCorr.at(i) += (HxTimeFrequency(j, i) * conj(HyTimeFrequency(j,
                i)));
            HxHxCorr.at(i) += (HxTimeFrequency(j, i) * conj(HxTimeFrequency(j,
                i)));
            HyHyCorr.at(i) += (HyTimeFrequency(j, i) * conj(HyTimeFrequency(j,
                i)));

          }

        Hdet.at(i) = (HxHxCorr.at(i) * HyHyCorr.at(i) - HxHyCorr.at(i) * conj(
            HxHyCorr.at(i)));
        logfile << "Hdet: " << Hdet.at(i) << endl;
        Zxx.at(i) = ((ExHxCorr.at(i) * HyHyCorr.at(i) - ExHyCorr.at(i) * conj(
            HxHyCorr.at(i))) / Hdet.at(i));
        Zxy.at(i) = ((ExHyCorr.at(i) * HxHxCorr.at(i) - ExHxCorr.at(i)
            * HxHyCorr.at(i)) / Hdet.at(i));
        Zyx.at(i) = ((EyHxCorr.at(i) * HyHyCorr.at(i) - EyHyCorr.at(i) * conj(
            HxHyCorr.at(i))) / Hdet.at(i));
        Zyy.at(i) = ((EyHyCorr.at(i) * HxHxCorr.at(i) - EyHxCorr.at(i)
            * HxHyCorr.at(i)) / Hdet.at(i));
      }
    Output.AssignAll(Zxx.size());
    for (size_t i = 0; i < Zxx.size(); ++i)
      {
        MTTensor
            CurrZ(Zxx.at(i), Zxy.at(i), Zyx.at(i), Zyy.at(i), i * freqstep);
        Output.SetMTData().at(i) = CurrZ;
      }
    vector<double> transferfunc(Zxx.size() * 2 - 1, 0.0);
    TsSpectrum().CalcTimeSeries(Zxy.begin(), Zxy.end(), transferfunc.begin(),
        transferfunc.end());
    ofstream tffile((datafilename + ".tf").c_str());
    copy(transferfunc.begin(), transferfunc.end(), ostream_iterator<double> (
        tffile, "\n"));

    Output.WriteAsMtt(datafilename.c_str());
  }