seispca.cpp

Go to the documentation of this file.

#include "PCA.h"
#include "statutils.h"
#include "CUniformRNG.h"
#include "VecMat.h"
#include "SeismicDataComp.h"
#include <iostream>
#include <fstream>
#include <vector>

using namespace gplib;

int main()
  {
    SeismicDataComp Radial, Tangential, Vertical;
    std::string radname, tangname, vername;

    std::cout << "Radial Component: ";
    std::cin >> radname;
    std::cout << "Tangential Component: ";
    std::cin >> tangname;
    std::cout << "Vertical Component: ";
    std::cin >> vername;

    Radial.ReadData(radname);
    Tangential.ReadData(tangname);
    Vertical.ReadData(vername);

    const size_t nobs = Radial.GetData().size();
    const size_t nchan = 3;
    gplib::rmat input(nchan, nobs);
    gplib::cmat evec(nchan, nchan);
    gplib::cvec eval(nchan);

    for (int i = 0; i < nobs; ++i)
      {
        input(0, i) = Radial.GetData().at(i);
        input(1, i) = Tangential.GetData().at(i);
        input(2, i) = Vertical.GetData().at(i);
      }

    PCA(input, evec, eval);

    gplib::cmat wmat(WhiteMat(evec, eval));
    std::cout << "pca evec: " << evec << std::endl;
    std::cout << "pca eval: " << eval << std::endl;
    std::cout << "pca WhM: " << wmat << std::endl;
    std::cout << "pca DeWhM: " << DeWhiteMat(evec, eval) << std::endl;

    gplib::cmat output(nchan, nobs);
    axpy_prod(wmat, input, output);
    SeismicDataComp RadOut(Radial), TangOut(Tangential), VerOut(Vertical);
    for (int i = 0; i < nobs; ++i)
      {
        RadOut.GetData().at(i) = output(0, i).real();
        TangOut.GetData().at(i) = output(1, i).real();
        VerOut.GetData().at(i) = output(2, i).real();
      }
    RadOut.WriteAsSac("pc1.out");
    TangOut.WriteAsSac("pc2.out");
    VerOut.WriteAsSac("pc3.out");

  }