GPLIB++
lmsprocessing.cpp
Go to the documentation of this file.
1 #include <fstream>
2 #include <iostream>
3 #include <vector>
4 #include <string>
5 #include <algorithm>
6 #include "CLMSCanceller.h"
7 #include "CTimeSeriesData.h"
8 #include "CStackedSpectrum.h"
9 #include "CMTStation.h"
10 #include "CMttData.h"
11 
12 using namespace std;
13 using namespace gplib;
14 
15 int main()
16  {
17  CLMSCanceller Canceller;
18  CTimeSeriesData TsData;
19  StackedSpectrum WeightSpectrum;
20  CMTStation MTData;
21  CMttData MttData;
22  //CMTStation Output;
23  MTData.MTData = &MttData;
24  const int filterlength = 111;
25  //const double mu = 1e-17;
26  const double mu = -1;
27  const int ntimeseries = 2;
28  const int shift = 0; // filterlength/2;
29  //const double freqstep = 150/filterlength;
30  double freqstep;
31  int i, j;
32  string tsfilename;
33  ifstream noisefile, tsfile;
34  ofstream outfile, weightfile, epsfile;
35  double temp;
36  vector<double> noise, epsilon, output, ts;
37 
38  cout << "Time Series Filename: ";
39  cin >> tsfilename;
40  TsData.GetData(tsfilename);
41  freqstep = TsData.Data->samplerate / filterlength;
42  Canceller.filterlength = filterlength;
43  Canceller.mu.assign(ntimeseries, -1);
44  //Canceller.input = &TsData.Data->Ex.data;
45  Canceller.FilterOutput = &(TsData.Data->Ex.data);
46  Canceller.Input.push_back(&TsData.Data->Hx.data);
47  Canceller.Input.push_back(&TsData.Data->Hy.data);
48  Canceller.Epsilon = &epsilon;
49  Canceller.Desired = &output;
50  Canceller.FilterData(shift);
51 
52  cout << "Mu: ";
53  for (i = 0; i < ntimeseries; ++i)
54  cout << Canceller.mu.at(i) << " ";
55  cout << endl;
56  epsfile.open((tsfilename + ".eps").c_str());
57  outfile.open((tsfilename + ".clean").c_str());
58  for (i = 0; i < output.size(); ++i)
59  {
60  outfile << output.at(i) << endl;
61  epsfile << Canceller.Epsilon->at(i) << endl;
62  }
63  outfile.close();
64  epsfile.close();
65 
66  WeightSpectrum.TimeSeries.assign(filterlength, 0);
67  copy(Canceller.Weights.begin() + filterlength, Canceller.Weights.end(),
68  WeightSpectrum.TimeSeries.begin());
69  WeightSpectrum.CalcSpectrum(filterlength);
70 
71  weightfile.open((tsfilename + ".weights").c_str());
72  for (i = 0; i < WeightSpectrum.Spectrum.size(); ++i)
73  weightfile << i << " " << abs(WeightSpectrum.Spectrum.at(i)) << " "
74  << arg(WeightSpectrum.Spectrum.at(i)) << endl;
75  weightfile.close();
76 
77  MTData.MTData->AssignAll(filterlength / 2);
78  cout << "Assigning Tensor values. " << endl << flush;
79  copy(WeightSpectrum.Spectrum.begin(), WeightSpectrum.Spectrum.begin()
80  + filterlength / 2, MTData.MTData->DataXY.Z.begin());
81  for (int i = 0; i < filterlength / 2; ++i)
82  MTData.MTData->frequency.at(i) = i * freqstep;
83  MTData.MTData->frequency.at(0) = 0.00000001;
84  cout << "Writing Data. " << endl << flush;
85  MTData.WriteAsMtt(tsfilename);
86  }
int main()
CMTStation MTData
Definition: cadijoint.cpp:15
void StackedSpectrum(InputIterator tsbegin, InputIterator tsend, OutputIterator freqbegin, const size_t seglength, WindowFunctype WFunc)
This template is used to calculate stacked spectra for example for power spectrum estimation...