seisadaptive.cpp

Go to the documentation of this file.

#include <fstream>
#include <iostream>
#include <vector>
#include <string>
#include <algorithm>
#include <boost/shared_ptr.hpp>
#include <boost/random/lagged_fibonacci.hpp>
#include <boost/random/normal_distribution.hpp>
#include <boost/random/variate_generator.hpp>
#include <boost/date_time/posix_time/posix_time.hpp>
#include "LMSCanceller.h"
#include "RLSCanceller.h"
#include "AMRLSCanceller.h"
#include "SeisAdaptiveConfig.h"
#include "SeismicDataComp.h"


using namespace std;
string version = "$Id: mtuadaptive.cpp 942 2007-04-08 16:13:01Z max $";



int main()
{
        cout << "This is mtuadaptive: Apply an adaptive filter (LMS or RLS) to cancel noise in MT time series" << endl<< endl;
        cout  << " The program will ask for reference and input filename, further settings are read from 'mtuadaptive.conf' " << endl;
        cout << " Output will be 2 Phoenix format files with ending '.clean' and '.eps' where the chosen channel is overwritten" << endl<< endl;
        cout << " This is Version: " << version << endl << endl;
        
        SeismicDataComp InputData, NoiseData;
        
        string tsfilename, noisefilename;
        SeisAdaptiveConfig Configuration;
        Configuration.GetData("mtuadaptive.conf");
        AdaptiveFilter *Filter;
        int filterchoice = -1;
        cout << "Which type of adaptive filter do you want ?" << endl;
        cout << " 1: LMS" << endl;
        cout << " 2: RLS" << endl;
        cout << " 3: AMRLS" << endl;
        cin >> filterchoice;
        switch (filterchoice)
        {
                case 1:
                        Filter = new LMSCanceller(Configuration.filterlength,Configuration.mu);
                break;
                case 2:
                        Filter = new RLSCanceller(Configuration.filterlength,Configuration.delta,Configuration.lambda);
                break;
                case 3:
                        Filter = new AMRLSCanceller(Configuration.filterlength,Configuration.delta,Configuration.lambda,Configuration.alpha);
                break;
                default:
                        cerr << "Invalid number input, don't know what to do and will exit ! ";
                        return 100;
                break;
        }
        
        cout << "Noise Data: ";
        cin >> noisefilename;
        NoiseData.GetData(noisefilename);
        
        const double orignoisemax = *max_element(NoiseData.GetData().begin(),NoiseData.GetData().end());
        boost::lagged_fibonacci607 RandomGenerator(static_cast<unsigned int>(std::time(0)));
        boost::normal_distribution<> Noisedist(0.0,fabs(orignoisemax/10.0));
        boost::variate_generator<boost::lagged_fibonacci607&, boost::normal_distribution<> > NoiseGenerator(RandomGenerator,Noisedist);

        const size_t noiselength = NoiseData.GetData().size();
        
        const size_t inputsize = Configuration.filterlength;
        const size_t outputsize = 1;
        gplib::rvec currentinput(inputsize), currdesired(outputsize), currout(outputsize);
        
        for (size_t i = 0; i < noiselength; ++i)
        {
                copy(NoiseData.GetData().begin()+i,NoiseData.GetData().begin()+inputsize+i,
                                currentinput.begin());
                currdesired(0) = NoiseGenerator();
                Filter->CalcOutput(currentinput,currout);
                Filter->AdaptFilter(currentinput,currdesired);  
        }       
        
        cout << "Radial Time Series Filename: ";
        cin >> tsfilename;
        InputData.GetData(tsfilename);
        const size_t inputlength = InputData.GetData().size();
        gplib::rvec newout(inputlength);
        for (size_t i = 0; i < inputlength; ++i)
        {
                copy(InputData.GetData().begin()+i,InputData.GetData().begin()+inputsize+i,
                                currentinput.begin());
                Filter->CalcOutput(currentinput,currout);
                newout(i) = currout(0);
        }
        copy(newout.begin(),newout.end(),InputData.GetData().begin());
        InputData.WriteAsSac("filter.r");
        ofstream rweightfile("filter.r.weights");
        copy(Filter->GetWeightsAsVector().begin(),Filter->GetWeightsAsVector().end(),ostream_iterator<double>(rweightfile,"\n"));
        
        cout << "Vertical Time Series Filename: ";
        cin >> tsfilename;
        InputData.GetData(tsfilename);
        
        for (size_t i = 0; i < inputlength; ++i)
        {
                copy(InputData.GetData().begin()+i,InputData.GetData().begin()+inputsize+i,
                                currentinput.begin());
                Filter->CalcOutput(currentinput,currout);
                newout(i) = currout(0);
        }
        copy(newout.begin(),newout.end(),InputData.GetData().begin());
        InputData.WriteAsSac("filter.z");
        ofstream zweightfile("filter.z.weights");
        copy(Filter->GetWeightsAsVector().begin(),Filter->GetWeightsAsVector().end(),ostream_iterator<double>(zweightfile,"\n"));
}

---