lmsprocessing.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 "CLMSCanceller.h"
00007 #include "CTimeSeriesData.h"
00008 #include "CStackedSpectrum.h"
00009 #include "CMTStation.h"
00010 #include "CMttData.h"
00011
00012 using namespace std;
00013 using namespace gplib;
00014
00015 int main()
00016 {
00017 CLMSCanceller Canceller;
00018 CTimeSeriesData TsData;
00019 StackedSpectrum WeightSpectrum;
00020 CMTStation MTData;
00021 CMttData MttData;
00022
00023 MTData.MTData = &MttData;
00024 const int filterlength = 111;
00025
00026 const double mu = -1;
00027 const int ntimeseries = 2;
00028 const int shift = 0;
00029
00030 double freqstep;
00031 int i, j;
00032 string tsfilename;
00033 ifstream noisefile, tsfile;
00034 ofstream outfile, weightfile, epsfile;
00035 double temp;
00036 vector<double> noise, epsilon, output, ts;
00037
00038 cout << "Time Series Filename: ";
00039 cin >> tsfilename;
00040 TsData.GetData(tsfilename);
00041 freqstep = TsData.Data->samplerate / filterlength;
00042 Canceller.filterlength = filterlength;
00043 Canceller.mu.assign(ntimeseries, -1);
00044
00045 Canceller.FilterOutput = &(TsData.Data->Ex.data);
00046 Canceller.Input.push_back(&TsData.Data->Hx.data);
00047 Canceller.Input.push_back(&TsData.Data->Hy.data);
00048 Canceller.Epsilon = ε
00049 Canceller.Desired = &output;
00050 Canceller.FilterData(shift);
00051
00052 cout << "Mu: ";
00053 for (i = 0; i < ntimeseries; ++i)
00054 cout << Canceller.mu.at(i) << " ";
00055 cout << endl;
00056 epsfile.open((tsfilename + ".eps").c_str());
00057 outfile.open((tsfilename + ".clean").c_str());
00058 for (i = 0; i < output.size(); ++i)
00059 {
00060 outfile << output.at(i) << endl;
00061 epsfile << Canceller.Epsilon->at(i) << endl;
00062 }
00063 outfile.close();
00064 epsfile.close();
00065
00066 WeightSpectrum.TimeSeries.assign(filterlength, 0);
00067 copy(Canceller.Weights.begin() + filterlength, Canceller.Weights.end(),
00068 WeightSpectrum.TimeSeries.begin());
00069 WeightSpectrum.CalcSpectrum(filterlength);
00070
00071 weightfile.open((tsfilename + ".weights").c_str());
00072 for (i = 0; i < WeightSpectrum.Spectrum.size(); ++i)
00073 weightfile << i << " " << abs(WeightSpectrum.Spectrum.at(i)) << " "
00074 << arg(WeightSpectrum.Spectrum.at(i)) << endl;
00075 weightfile.close();
00076
00077 MTData.MTData->AssignAll(filterlength / 2);
00078 cout << "Assigning Tensor values. " << endl << flush;
00079 copy(WeightSpectrum.Spectrum.begin(), WeightSpectrum.Spectrum.begin()
00080 + filterlength / 2, MTData.MTData->DataXY.Z.begin());
00081 for (int i = 0; i < filterlength / 2; ++i)
00082 MTData.MTData->frequency.at(i) = i * freqstep;
00083 MTData.MTData->frequency.at(0) = 0.00000001;
00084 cout << "Writing Data. " << endl << flush;
00085 MTData.WriteAsMtt(tsfilename);
00086 }