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 relnoise, absnoise;

    string version = "$Id: seisnoise.cpp 1882 2014-06-12 08:15:20Z 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 and as a minimum absolute value."
        << endl;
    cout << endl << endl;
    if (argc > 4)
      {
        datafilename = argv[1];
        outfilename = argv[2];
        relnoise = atof(argv[3]);
        absnoise = atof(argv[4]);
      }
    else
      {
        datafilename = AskFilename("Input filename: ");
        outfilename = AskFilename("Output filename: ");
        cout << "Relative Noise level: ";
        cin >> relnoise;
        cout << "Absolute Noise level: ";
        cin >> absnoise;
      }
    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, std::max(relnoise * maxelem, absnoise));
    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);
  }