mtulmscancel.cpp

Go to the documentation of this file.

#include <fstream>
#include <iostream>
#include <vector>
#include <string>
#include <algorithm>
#include "LMSCanceller.h"
#include "RLSCanceller.h"
#include "TimeSeriesData.h"
#include "CLmsConfig.h"
#include "ApplyFilter.h"
#include "netcdfcpp.h"

using namespace std;
string version = "$Id: mtulmscancel.cpp 865 2007-03-02 17:47:39Z max $";


void WriteNetCDF(const string name,const ublas::matrix<double,boost::numeric::ublas::column_major> &Matrix)
{
        NcFile mtrescdf(name.c_str(),NcFile::Replace);
        NcDim* xd = mtrescdf.add_dim("x",Matrix.size1());
        NcDim* yd = mtrescdf.add_dim("y",Matrix.size2());
        NcVar* x = mtrescdf.add_var("x",ncFloat,xd);
        NcVar* y = mtrescdf.add_var("y",ncFloat,yd);
        NcVar* z = mtrescdf.add_var("z",ncFloat,xd,yd);
        float *xvals = new float[xd->size()];
        float *yvals = new float[yd->size()];
        float *zvals = new float[xd->size()*yd->size()];
        for (size_t i = 0; i < Matrix.size1(); ++i)
        {
                xvals[i] = i;
                for (size_t j = 0; j < Matrix.size2(); ++j)
                {
                        yvals[j] = j;
                        zvals[(Matrix.size1()-1-i)*Matrix.size2() +j] = Matrix(i,j);
                }
        }       
        x->put(xvals,xd->size());
        y->put(yvals,yd->size());
        z->put(zvals,z->edges());
}


int main()
{
        cout << "This is lmcancel: Use an lms-adaptive filter to cancel noise in MT time series" << endl<< endl;
        cout  << " The program will ask for reference and input filename, further settings are read from 'lms.conf' " << endl;
        cout << " Output will be 2 Phoenix format files with ending '.clean' and '.eps'" << endl<< endl;
        cout << " This is Version: " << version << endl << endl;
        
        CLmsConfig Configuration;
        TimeSeriesData InputData,ReferenceData;
        
        int index;
        string tsfilename, noisefilename;
        
        int ntimeseries;
        Configuration.GetData("lms.conf");
        LMSCanceller Filter(Configuration.filterlength); 
        Filter.SetMu(Configuration.mu);
        
        ApplyFilter Canceller(Filter,true);
        
        
        cout << "Reference Data: ";
        cin >> noisefilename;
        ReferenceData.GetData(noisefilename);
        cout << "Component Index (Hx,Hy,Ex,Ey): ";
        cin >> index;
        TimeSeriesComponent *RefComp;
        switch (index)
        {
                case 1:
                        RefComp = &ReferenceData.GetData().GetHx();
                break;
                case 2:
                        RefComp = &ReferenceData.GetData().GetHy();
                break;
                case 3:
                        RefComp = &ReferenceData.GetData().GetEx();
                break;
                case 4:
                        RefComp = &ReferenceData.GetData().GetEy();
                break;
                default:
                        cerr << "Component index not valid !";
                        return 100;
                break;
        }
        Canceller.AddReferenceChannel(*RefComp);
        cout << "Time Series Filename: ";
        cin >> tsfilename;
        InputData.GetData(tsfilename);
        cout << "Number of InputData: ";
        cin >> ntimeseries;
        
        for (int i=0; i < ntimeseries; ++i)
        {
                cout << "Component Index (Hx,Hy,Ex,Ey): ";
                cin >> index;
                switch (index)
                {
                        case 1:
                                Canceller.AddInputChannel(InputData.GetData().GetHx());
                        break;
                        case 2:
                                Canceller.AddInputChannel(InputData.GetData().GetHy());
                        break;
                        case 3:
                                Canceller.AddInputChannel(InputData.GetData().GetEx());
                        break;
                        case 4:
                                Canceller.AddInputChannel(InputData.GetData().GetEy());
                        break;
                        default:
                                cerr << "Component index not valid !";
                                return 100;
                        break;
                }
        }
        Canceller.SetWeightSaveIntervall(1000);
        Canceller.SetShift(Configuration.shift);
        Canceller.FilterData();
        Canceller.FilterData();
        
        //*RefComp = *Canceller.GetOutChannels().front();
        copy(Canceller.GetOutChannels().front()->GetData().begin(),Canceller.GetOutChannels().front()->GetData().end(),RefComp->GetData().begin());
        ReferenceData.WriteAsMtu(noisefilename+".clean");
        ofstream cleanfile((noisefilename+".clean").c_str());
        copy(Canceller.GetOutChannels().front()->GetData().begin(),Canceller.GetOutChannels().front()->GetData().end(),ostream_iterator<double>(cleanfile,"\n"));
        copy(Canceller.GetEpsValues().front().begin(),Canceller.GetEpsValues().front().end(),RefComp->GetData().begin());
        
        ReferenceData.WriteAsMtu(noisefilename+".eps");
        ofstream epsfile((noisefilename+".eps").c_str());
        copy(Canceller.GetEpsValues().front().begin(),Canceller.GetEpsValues().front().end(),ostream_iterator<double>(epsfile,"\n"));
        
        //cout << "Last Weights: " << endl;
        //Filter.PrintWeights(cout);
        //cout << endl;
        ofstream weightfile((noisefilename+".weights").c_str());
        copy(Filter.GetWeights().begin(),Filter.GetWeights().end(),ostream_iterator<double>(weightfile,"\n"));
        WriteNetCDF(noisefilename+".weights.nc",Canceller.GetWeightHistory());
}