seisnoise.cpp

Go to the documentation of this file.

#include <iostream>
#include <string>
#include "SeismicDataComp.h"
#include <algorithm>
#include "Util.h"
#include <boost/random/lagged_fibonacci.hpp>
#include <boost/random/normal_distribution.hpp>
#include <boost/random/variate_generator.hpp>
#include <boost/bind.hpp>

using namespace std;
using namespace gplib;

/*! \file 
 * This program adds random gaussian noise to a single sac file. The noise
 * level is specified relative to the maximum amplitude in the time series. The noise
 * is completely random and uncorrelated between datapoints. This might not be very
 * typical for seismic noise that is often signal generated and has characteristic 
 * frequencies, but seems to be sufficient to perform basic tests on synthetic data.
 */
int main(int argc, char* argv[])
  {
    SeismicDataComp Data;
    string datafilename, outfilename;
    double noiselevel;

    string version = "$Id: seisnoise.cpp 1816 2009-09-07 11:28:35Z mmoorkamp $";
    cout << endl << endl;
    cout << "Program " << version << endl;
    cout << " Add random noise to a seismogram" << endl;
    cout
        << " Input is a single component file, output a sac file with added noise"
        << endl;
    cout << " The noise level is specified relative to the maximum amplitude"
        << endl;
    cout << endl << endl;
    if (argc > 3)
      {
        datafilename = argv[1];
        outfilename = argv[2];
        noiselevel = atof(argv[3]);
      }
    else
      {
        datafilename = AskFilename("Input filename: ");
        outfilename = AskFilename("Output filename: ");
        cout << "Noise level: ";
        cin >> noiselevel;
      }
    Data.ReadData(datafilename);
    const double maxelem = *max_element(Data.GetData().begin(), Data.GetData().end());
    boost::lagged_fibonacci607
        generator(static_cast<unsigned int>(std::time(0)));
    boost::normal_distribution<> noise_dist(0.0, noiselevel*maxelem);
    boost::variate_generator<boost::lagged_fibonacci607&, boost::normal_distribution<> >
        noise(generator, noise_dist);

    //add noise to each point of the time series
    transform(Data.GetData().begin(), Data.GetData().end(), Data.GetData().begin(), boost::bind(plus<double>(), _1, boost::bind(noise)));
    Data.WriteAsSac(outfilename);
  }