corrselect.cpp

Go to the documentation of this file.

#include "SeismicStationList.h"
#include "VecMat.h"
#include "miscfunc.h"
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include "Util.h"

using namespace std;
using namespace gplib;

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

gplib::rmat CalcCorrMatrix(const SeismicStationList::tseiscompvector &List)
  {
    const unsigned int nrf = List.size();
    gplib::rmat CorrMatrix(nrf, nrf);
    CorrMatrix *= 0.0;
    for (unsigned int i = 0; i < nrf; ++i)
      {
        for (unsigned int j = i; j < nrf; ++j)
          {
            CorrMatrix(i, j) = Cross(List.at(i)->GetData(),
                List.at(j)->GetData(), 0, min(List.at(i)->GetData().size(),
                    List.at(j)->GetData().size()));
            CorrMatrix(j, i) = CorrMatrix(i, j);
          }
      }
    return CorrMatrix;
  }

int main()
  {
    cout
        << " This is corrselect: Select receiver functions based on their correlation"
        << endl; // write some info
    cout << " Reads in receiver function in SAC format" << endl;
    cout << " This is Version: " << version << endl << endl;

    SeismicStationList RF;

    string reclistname = AskFilename("File with list of filenames: ");
    cin >> reclistname;

    RF.ReadList(reclistname);

    unsigned int nrf = RF.GetList().size();
    cout << "Read " << nrf << " receiver functions." << endl;
    string cont = "y";
    while (cont != "n")
      {
        gplib::rmat CorrMatrix = CalcCorrMatrix(RF.GetList());
        nrf = RF.GetList().size();
        for (unsigned int j = 0; j < nrf; ++j)
          cout << j << " " << CorrMatrix(0, j) << endl;
        ublas::scalar_vector<double> mul(nrf, 1. / nrf);
        gplib::rvec AvgCorr = prod(CorrMatrix, mul);

        cout << "Average Correlation: ";
        copy(AvgCorr.begin(), AvgCorr.end(), ostream_iterator<double> (cout,
            " "));
        cout << endl;
        cout << "Threshold: ";
        double threshold = 0.0;
        cin >> threshold;
        SeismicStationList::tseiscompvector::iterator it = RF.GetList().begin();
        for (unsigned int j = 0; j < nrf; ++j, ++it)
          {
            if (AvgCorr(j) < threshold)
              RF.GetList().erase(it);
          }
        cout << "Continue (y/n)";
        cin >> cont;
      }
    nrf = RF.GetList().size();
    for (unsigned int j = 0; j < nrf; ++j)
      {
        cout << RF.GetList().at(j)->GetName() << endl;
      }
  }