mtubandpass.cpp

Go to the documentation of this file.

#include <iostream>
#include <fstream>
#include <string>
#include <numeric>
#include "Util.h"
#include "TimeSeriesData.h"
#include "statutils.h"
#include "FilterFunc.h"

using namespace std;
using namespace gplib;

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

/*!
 * \addtogroup UtilProgs Utility Programs
 *@{
 * \file
 * Apply a band pass filter to each component of the MT time series. The program
 * asks for the corner frequencies in Hz and the number of passes for the filter.
 * The output filename is determined by appending .fil to the input filename.
 */

int main(int argc, char *argv[])
  {
    string infilename;
    cout << "This is mtufilter: Apply a band pass filter to MT time series"
        << endl << endl;
    cout << " Usage:      mtufilter infilename " << endl;
    cout << " Ending '.fil'  will be automatically assigned to outfilename"
        << endl << endl;
    cout << " This is Version: " << version << endl << endl;
    if (argc == 2)
      {
        infilename = argv[1];
      }
    else
      {

        infilename = AskFilename(" Mtu-Filename: ");
      }

    TimeSeriesData Data;
    Data.GetData(infilename);
    // get the lower corner frequency
    double lowfreq;
    cout << "Lower Corner frequency [Hz]: ";
    cin >> lowfreq;

    // get the lower corner frequency
    double upfreq;
    cout << "Upper Corner frequency [Hz]: ";
    cin >> upfreq;

    // calculate the dimensionless corner frequency for the filtering class
    const double lowfilfreq = lowfreq * Data.GetData().GetEx().GetSamplerate();
    const double upfilfreq = upfreq * Data.GetData().GetEx().GetSamplerate();
    // the more passes we do with this filter, the stronger the attenuation
    size_t npass = 1;
    cout << "Number of passes: ";
    cin >> npass;

    // make sure the all components have zero mean
    SubMean(Data.GetData().GetEx().GetData().begin(),
        Data.GetData().GetEx().GetData().end());
    SubMean(Data.GetData().GetEy().GetData().begin(),
        Data.GetData().GetEy().GetData().end());
    SubMean(Data.GetData().GetHx().GetData().begin(),
        Data.GetData().GetHx().GetData().end());
    SubMean(Data.GetData().GetHy().GetData().begin(),
        Data.GetData().GetHy().GetData().end());
    SubMean(Data.GetData().GetHz().GetData().begin(),
        Data.GetData().GetHz().GetData().end());

    //we can do several passes of the filter to get strong attenuation of high frequencies
    for (size_t i = 0; i < npass; ++i)
      {
        // apply the filter to one component and store the result in that component
        transform(Data.GetData().GetEx().GetData().begin(),
            Data.GetData().GetEx().GetData().end(),
            Data.GetData().GetEx().GetData().begin(), SimpleBp(lowfilfreq,upfilfreq));
        transform(Data.GetData().GetEy().GetData().begin(),
            Data.GetData().GetEy().GetData().end(),
            Data.GetData().GetEy().GetData().begin(), SimpleBp(lowfilfreq,upfilfreq));
        transform(Data.GetData().GetHx().GetData().begin(),
            Data.GetData().GetHx().GetData().end(),
            Data.GetData().GetHx().GetData().begin(), SimpleBp(lowfilfreq,upfilfreq));
        transform(Data.GetData().GetHy().GetData().begin(),
            Data.GetData().GetHy().GetData().end(),
            Data.GetData().GetHy().GetData().begin(), SimpleBp(lowfilfreq,upfilfreq));
        transform(Data.GetData().GetHz().GetData().begin(),
            Data.GetData().GetHz().GetData().end(),
            Data.GetData().GetHz().GetData().begin(), SimpleBp(lowfilfreq,upfilfreq));
      }
    // write data in the same format that we read it in
    Data.WriteBack(infilename + ".fil");
  }
/*@}*/