mtuadaptive.cpp

Go to the documentation of this file.
00001 #include <fstream>
00002 #include <iostream>
00003 #include <vector>
00004 #include <string>
00005 #include <algorithm>
00006 #include <boost/shared_ptr.hpp>
00007 #include <boost/date_time/posix_time/posix_time.hpp>
00008 #include "LMSCanceller.h"
00009 #include "RLSCanceller.h"
00010 #include "AMRLSCanceller.h"
00011 #include "TimeSeriesData.h"
00012 #include "CMtuAdaptiveConfig.h"
00013 #include "ApplyFilter.h"
00014 #include "netcdfcpp.h"
00015 #include "NeuralNetwork.h"
00016 
00017 using namespace std;
00018 string version = "$Id: mtuadaptive.cpp 942 2007-04-08 16:13:01Z max $";
00019 
00020 void GetNNSetup(CMtuAdaptiveConfig &Configuration, const size_t ntimeseries, 
00021         NeuralNetwork::ttypeArray &NNLayers,double &NNmaxinit)
00022 {
00023         NeuralNetwork::ttypeVector typeVector(Configuration.filterlength,SigmoidalNeuron::bipolar); // we want filterlength number of bipolar neurons per hidden layer 
00024         for (int i = 0; i < Configuration.hiddenlayers; ++i) //intialize the type array for the hidden layers
00025         {
00026                 NNLayers.push_back(typeVector); // all layers are the same, so we copy the same vector there
00027         }
00028         typeVector.assign(1,SigmoidalNeuron::identity); 
00029         NNLayers.push_back(typeVector); // and then we add it to the type Array
00030 }
00031 
00032 void WriteNetCDF(const string name,const ublas::matrix<double,boost::numeric::ublas::column_major> &Matrix)
00033 {
00034         NcFile mtrescdf(name.c_str(),NcFile::Replace);
00035         NcDim* xd = mtrescdf.add_dim("x",Matrix.size1());
00036         NcDim* yd = mtrescdf.add_dim("y",Matrix.size2());
00037         NcVar* x = mtrescdf.add_var("x",ncFloat,xd);
00038         NcVar* y = mtrescdf.add_var("y",ncFloat,yd);
00039         NcVar* z = mtrescdf.add_var("z",ncFloat,xd,yd);
00040         float *xvals = new float[xd->size()];
00041         float *yvals = new float[yd->size()];
00042         float *zvals = new float[xd->size()*yd->size()];
00043         for (size_t i = 0; i < Matrix.size1(); ++i)
00044         {
00045                 xvals[i] = i;
00046                 for (size_t j = 0; j < Matrix.size2(); ++j)
00047                 {
00048                         yvals[j] = j;
00049                         zvals[(Matrix.size1()-1-i)*Matrix.size2() +j] = Matrix(i,j);
00050                 }
00051         }       
00052         x->put(xvals,xd->size());
00053         y->put(yvals,yd->size());
00054         z->put(zvals,z->edges());
00055 }
00056 
00057 
00058 int main()
00059 {
00060         cout << "This is mtuadaptive: Apply an adaptive filter (LMS or RLS) to cancel noise in MT time series" << endl<< endl;
00061         cout  << " The program will ask for reference and input filename, further settings are read from 'mtuadaptive.conf' " << endl;
00062         cout << " Output will be 2 Phoenix format files with ending '.clean' and '.eps' where the chosen channel is overwritten" << endl<< endl;
00063         cout << " This is Version: " << version << endl << endl;
00064         
00065         CMtuAdaptiveConfig Configuration;
00066         TimeSeriesData InputData,ReferenceData;
00067         
00068         int index;
00069         string tsfilename, noisefilename;
00070         
00071         int ntimeseries;
00072         Configuration.GetData("mtuadaptive.conf");
00073         AdaptiveFilter *Filter;
00074         int filterchoice = -1;
00075         cout << "Which type of adaptive filter do you want ?" << endl;
00076         cout << " 1: LMS" << endl;
00077         cout << " 2: RLS" << endl;
00078         cout << " 3: AMRLS" << endl;
00079         cout << " 4: NN" << endl;
00080         cin >> filterchoice;
00081         
00082         cout << "Reference Data: ";
00083         cin >> noisefilename;
00084         ReferenceData.GetData(noisefilename);
00085         cout << "Component Index (Hx,Hy,Ex,Ey): ";
00086         cin >> index;
00087         cout << "Input Time Series Filename: ";
00088         cin >> tsfilename;
00089         InputData.GetData(tsfilename);
00090         cout << "Input Start time: " << InputData.GetData().GetTime().front() << endl;
00091         cout << "Reference Start time: " << ReferenceData.GetData().GetTime().front() << endl;
00092         if (InputData.GetData().GetTime().front() != ReferenceData.GetData().GetTime().front())
00093         {
00094                 cerr << "Time series not synchronized !" << endl;
00095                 return 100;
00096         }
00097         int lengthdiff = ReferenceData.GetData().Size() - InputData.GetData().Size();
00098         if (lengthdiff > 0)
00099         {
00100                 cout << "Removing " << lengthdiff << " datapoints from Reference time-series." <<endl; 
00101                 ReferenceData.GetData().erase(ReferenceData.GetData().Size() - lengthdiff,ReferenceData.GetData().Size());
00102         }
00103         if (lengthdiff < 0)
00104         {
00105                 cout << "Removing " << lengthdiff << " datapoints from Input time-series." <<endl;
00106                 InputData.GetData().erase(InputData.GetData().Size() + lengthdiff,InputData.GetData().Size());
00107         }
00108         cout << "Input End time: " << InputData.GetData().GetTime().back() << endl;
00109         cout << "Reference End time: " << ReferenceData.GetData().GetTime().back() << endl;
00110                 
00111         TimeSeriesComponent *RefComp;
00112         switch (index)
00113         {
00114                 case 1:
00115                         RefComp = &ReferenceData.GetData().GetHx();
00116                 break;
00117                 case 2:
00118                         RefComp = &ReferenceData.GetData().GetHy();
00119                 break;
00120                 case 3:
00121                         RefComp = &ReferenceData.GetData().GetEx();
00122                 break;
00123                 case 4:
00124                         RefComp = &ReferenceData.GetData().GetEy();
00125                 break;
00126                 default:
00127                         cerr << "Component index not valid !";
00128                         return 100;
00129                 break;
00130         }
00131         
00132 
00133         
00134         cout << "Number of input channels: ";
00135         cin >> ntimeseries;
00136         
00137         NeuralNetwork::ttypeArray NNLayers;
00138         double NNmaxinit;
00139         switch (filterchoice)
00140         {
00141                 case 1:
00142                         Filter = new LMSCanceller(Configuration.filterlength,Configuration.mu);
00143                 break;
00144                 case 2:
00145                         Filter = new RLSCanceller(Configuration.filterlength,Configuration.delta,Configuration.lambda);
00146                 break;
00147                 case 3:
00148                         Filter = new AMRLSCanceller(Configuration.filterlength,Configuration.delta,Configuration.lambda,Configuration.alpha);
00149                 break;
00150                 case 4:
00151                         GetNNSetup(Configuration,ntimeseries,NNLayers,NNmaxinit);
00152                         Filter = new NeuralNetwork(Configuration.filterlength,1,Configuration.mu,NNLayers,NNmaxinit);
00153                 break;
00154                 default:
00155                         cerr << "Invalid number input, don't know what to do and will exit ! ";
00156                         return 100;
00157                 break;
00158         }
00159         
00160         ApplyFilter Canceller(*Filter,true);
00161         Canceller.AddReferenceChannel(*RefComp);
00162         for (int i=0; i < ntimeseries; ++i)
00163         {
00164                 cout << "Component Index (Hx,Hy,Ex,Ey): ";
00165                 cin >> index;
00166                 switch (index)
00167                 {
00168                         case 1:
00169                                 Canceller.AddInputChannel(InputData.GetData().GetHx());
00170                         break;
00171                         case 2:
00172                                 Canceller.AddInputChannel(InputData.GetData().GetHy());
00173                         break;
00174                         case 3:
00175                                 Canceller.AddInputChannel(InputData.GetData().GetEx());
00176                         break;
00177                         case 4:
00178                                 Canceller.AddInputChannel(InputData.GetData().GetEy());
00179                         break;
00180                         default:
00181                                 cerr << "Component index not valid !";
00182                                 return 100;
00183                         break;
00184                 }
00185         }
00186         
00187         
00188         
00189         
00190         Canceller.SetWeightSaveIntervall(1000);
00191         Canceller.SetShift(Configuration.shift);
00192         Canceller.ShowProgress(true);
00193         cout << " First iteration: " << endl << endl;
00194         
00195         Canceller.FilterData();
00196         
00197         cout << endl << endl << " Second iteration: " << endl << endl;
00198         
00199         Canceller.FilterData();
00200         
00201         copy(Canceller.GetOutChannels().front()->GetData().begin(),Canceller.GetOutChannels().front()->GetData().end(),RefComp->GetData().begin());
00202         ReferenceData.WriteAsMtu(noisefilename+".clean");
00203         ofstream cleanfile((noisefilename+".clean").c_str());
00204         copy(Canceller.GetOutChannels().front()->GetData().begin(),Canceller.GetOutChannels().front()->GetData().end(),ostream_iterator<double>(cleanfile,"\n"));
00205         copy(Canceller.GetEpsValues().front().begin(),Canceller.GetEpsValues().front().end(),RefComp->GetData().begin());
00206         
00207         ReferenceData.WriteAsMtu(noisefilename+".eps");
00208         ofstream epsfile((noisefilename+".eps").c_str());
00209         copy(Canceller.GetEpsValues().front().begin(),Canceller.GetEpsValues().front().end(),ostream_iterator<double>(epsfile,"\n"));
00210         
00211         WriteNetCDF(noisefilename+".weights.nc",Canceller.GetWeightHistory());
00212         delete Filter;
00213 }

Generated on Fri Jul 4 15:30:21 2008 for GPLIB++ by  doxygen 1.5.5