reconstitute.cpp

Go to the documentation of this file.
00001 #include <iostream>
00002 #include "CTsSpectrum.h"
00003 #include "CMtuFormat.h"
00004 #include "CMttData.h"
00005 #include "types.h"
00006 using namespace std;
00007 int main()
00008 {
00009         CMtuFormat      OriginalData;
00010         CMtuFormat      NewData;
00011         CMtuFormat      DiffData;
00012         CMttData        ZInput;
00013         CTsSpectrum ExSpectrum;
00014         CTsSpectrum EySpectrum;
00015         CTsSpectrum HxSpectrum;
00016         CTsSpectrum HySpectrum;
00017         string tsfilename;
00018         string zfilename;
00019         int seglength = 0;
00020         int nsegs = 0;
00021         double exmean, eymean, hxmean, hymean;
00022         double freqstep;
00023         
00024         
00025         cout << "Time Series File: ";
00026         cin >> tsfilename;
00027         OriginalData.GetData(tsfilename.c_str());
00028         NewData.GetData(tsfilename.c_str());
00029         DiffData.GetData(tsfilename.c_str());
00030         
00031         NewData.Ex.data.clear();
00032         NewData.Ey.data.clear();
00033         DiffData.Ex.data.clear();
00034         DiffData.Ey.data.clear();
00035         
00036         
00037         cout << "Impedance File: ";
00038         cin >> zfilename;
00039         ZInput.GetData(zfilename.c_str());
00040         seglength = ZInput.frequency.size()*2;
00041         cout << "Seglength : " << seglength << endl;
00042         
00043         ExSpectrum.TsData = new ttsdata;
00044         EySpectrum.TsData = new ttsdata;
00045         HxSpectrum.TsData = &OriginalData.Hx.data;
00046         HySpectrum.TsData = &OriginalData.Hy.data;
00047         
00048         ExSpectrum.FreqData = new tcompdata;
00049         EySpectrum.FreqData = new tcompdata;
00050         HxSpectrum.FreqData = new tcompdata;
00051         HySpectrum.FreqData = new tcompdata;
00052         
00053         ExSpectrum.MultiCalc = false;
00054         EySpectrum.MultiCalc = false;
00055         HxSpectrum.MultiCalc = false;
00056         HySpectrum.MultiCalc = false;
00057         
00058         
00059         nsegs = HxSpectrum.TsData->size()/seglength-1;
00060         freqstep = 1./(seglength * OriginalData.samplerate);
00061         
00062         for (int j = 0; j < nsegs; ++j)
00063         {
00064                 exmean = 0;
00065                 eymean = 0;
00066                 hxmean = 0;
00067                 hymean = 0;
00068                 /*for (int k = j*seglength; k < (j+1)*seglength; ++k)
00069                 {
00070                         HxSpectrum.TsData->at(k) -= HxSpectrum.TsData->at(j*seglength) 
00071                                 +(HxSpectrum.TsData->at((j+1)*seglength) - HxSpectrum.TsData->at(j*seglength))* (k - j*seglength)/seglength;
00072                         HySpectrum.TsData->at(k) -= HySpectrum.TsData->at(j*seglength) 
00073                                 +(HySpectrum.TsData->at((j+1)*seglength) - HySpectrum.TsData->at(j*seglength))* (k - j*seglength)/seglength;
00074                         
00075                 }*/
00076                 /*for (int k = j*seglength; k < (j+1)*seglength; ++k)
00077                 {
00078                         hxmean += HxSpectrum.TsData->at(k);
00079                         hymean += HySpectrum.TsData->at(k);
00080                 }
00081                 for (int k = j*seglength; k < (j+1)*seglength; ++k)
00082                 {
00083                         HxSpectrum.TsData->at(k) -= hxmean;
00084                         HySpectrum.TsData->at(k) -= hymean;
00085                         HxSpectrum.TsData->at(k) *= 0.5 * (1 - cos(2 * PI *(k-seglength) /seglength));
00086                         HySpectrum.TsData->at(k) *= 0.5 * (1 - cos(2 * PI *(k-seglength) /seglength));
00087                 }*/
00088                 HxSpectrum.CalcSpectrum(j*seglength,(j+1)*seglength);
00089                 HySpectrum.CalcSpectrum(j*seglength,(j+1)*seglength);
00090                 ExSpectrum.FreqData->assign(seglength/2,0);
00091                 EySpectrum.FreqData->assign(seglength/2,0);
00092                 ExSpectrum.TsData->assign(seglength,0);
00093                 EySpectrum.TsData->assign(seglength,0);
00094                 for (int k = 0; k < HxSpectrum.FreqData->size()-1; ++k)
00095                 {
00096                         ExSpectrum.FreqData->at(k) = HxSpectrum.FreqData->at(k) * ZInput.DataXX.Z.at(k)
00097                                 + HySpectrum.FreqData->at(k) * ZInput.DataXY.Z.at(k);
00098                         EySpectrum.FreqData->at(k) = HxSpectrum.FreqData->at(k) * ZInput.DataYX.Z.at(k)
00099                                 + HySpectrum.FreqData->at(k) * ZInput.DataYY.Z.at(k);
00100                 }
00101                 ExSpectrum.CalcTimeSeries();
00102                 EySpectrum.CalcTimeSeries();
00103                 for (int k = 0; k < ExSpectrum.TsData->size(); ++k)
00104                 {
00105                         NewData.Ex.data.push_back(ExSpectrum.TsData->at(k));
00106                         NewData.Ey.data.push_back(EySpectrum.TsData->at(k));
00107                                 
00108                 }
00109         }       
00110         NewData.Ex.data.insert(NewData.Ex.data.end(),OriginalData.Hx.data.size() - NewData.Ex.data.size(),0);
00111         NewData.Ey.data.insert(NewData.Ey.data.end(),OriginalData.Hy.data.size() - NewData.Ey.data.size(),0);
00112         for (int k = 0; k < NewData.Ex.data.size(); ++k)
00113         {
00114                 DiffData.Ex.data.push_back(NewData.Ex.data.at(k)-OriginalData.Ex.data.at(k));
00115                 DiffData.Ey.data.push_back(NewData.Ey.data.at(k)-OriginalData.Ey.data.at(k));
00116         } 
00117         NewData.WriteData("out.ts");
00118         DiffData.WriteData("diff.ts");
00119 }