reconstitute.cpp

Go to the documentation of this file.

#include <iostream>
#include "CTsSpectrum.h"
#include "CMtuFormat.h"
#include "CMttData.h"
#include "types.h"
using namespace std;
int main()
{
        CMtuFormat      OriginalData;
        CMtuFormat      NewData;
        CMtuFormat      DiffData;
        CMttData        ZInput;
        CTsSpectrum ExSpectrum;
        CTsSpectrum EySpectrum;
        CTsSpectrum HxSpectrum;
        CTsSpectrum HySpectrum;
        string tsfilename;
        string zfilename;
        int seglength = 0;
        int nsegs = 0;
        double exmean, eymean, hxmean, hymean;
        double freqstep;
        
        
        cout << "Time Series File: ";
        cin >> tsfilename;
        OriginalData.GetData(tsfilename.c_str());
        NewData.GetData(tsfilename.c_str());
        DiffData.GetData(tsfilename.c_str());
        
        NewData.Ex.data.clear();
        NewData.Ey.data.clear();
        DiffData.Ex.data.clear();
        DiffData.Ey.data.clear();
        
        
        cout << "Impedance File: ";
        cin >> zfilename;
        ZInput.GetData(zfilename.c_str());
        seglength = ZInput.frequency.size()*2;
        cout << "Seglength : " << seglength << endl;
        
        ExSpectrum.TsData = new ttsdata;
        EySpectrum.TsData = new ttsdata;
        HxSpectrum.TsData = &OriginalData.Hx.data;
        HySpectrum.TsData = &OriginalData.Hy.data;
        
        ExSpectrum.FreqData = new tcompdata;
        EySpectrum.FreqData = new tcompdata;
        HxSpectrum.FreqData = new tcompdata;
        HySpectrum.FreqData = new tcompdata;
        
        ExSpectrum.MultiCalc = false;
        EySpectrum.MultiCalc = false;
        HxSpectrum.MultiCalc = false;
        HySpectrum.MultiCalc = false;
        
        
        nsegs = HxSpectrum.TsData->size()/seglength-1;
        freqstep = 1./(seglength * OriginalData.samplerate);
        
        for (int j = 0; j < nsegs; ++j)
        {
                exmean = 0;
                eymean = 0;
                hxmean = 0;
                hymean = 0;
                /*for (int k = j*seglength; k < (j+1)*seglength; ++k)
                {
                        HxSpectrum.TsData->at(k) -= HxSpectrum.TsData->at(j*seglength) 
                                +(HxSpectrum.TsData->at((j+1)*seglength) - HxSpectrum.TsData->at(j*seglength))* (k - j*seglength)/seglength;
                        HySpectrum.TsData->at(k) -= HySpectrum.TsData->at(j*seglength) 
                                +(HySpectrum.TsData->at((j+1)*seglength) - HySpectrum.TsData->at(j*seglength))* (k - j*seglength)/seglength;
                        
                }*/
                /*for (int k = j*seglength; k < (j+1)*seglength; ++k)
                {
                        hxmean += HxSpectrum.TsData->at(k);
                        hymean += HySpectrum.TsData->at(k);
                }
                for (int k = j*seglength; k < (j+1)*seglength; ++k)
                {
                        HxSpectrum.TsData->at(k) -= hxmean;
                        HySpectrum.TsData->at(k) -= hymean;
                        HxSpectrum.TsData->at(k) *= 0.5 * (1 - cos(2 * PI *(k-seglength) /seglength));
                        HySpectrum.TsData->at(k) *= 0.5 * (1 - cos(2 * PI *(k-seglength) /seglength));
                }*/
                HxSpectrum.CalcSpectrum(j*seglength,(j+1)*seglength);
                HySpectrum.CalcSpectrum(j*seglength,(j+1)*seglength);
                ExSpectrum.FreqData->assign(seglength/2,0);
                EySpectrum.FreqData->assign(seglength/2,0);
                ExSpectrum.TsData->assign(seglength,0);
                EySpectrum.TsData->assign(seglength,0);
                for (int k = 0; k < HxSpectrum.FreqData->size()-1; ++k)
                {
                        ExSpectrum.FreqData->at(k) = HxSpectrum.FreqData->at(k) * ZInput.DataXX.Z.at(k)
                                + HySpectrum.FreqData->at(k) * ZInput.DataXY.Z.at(k);
                        EySpectrum.FreqData->at(k) = HxSpectrum.FreqData->at(k) * ZInput.DataYX.Z.at(k)
                                + HySpectrum.FreqData->at(k) * ZInput.DataYY.Z.at(k);
                }
                ExSpectrum.CalcTimeSeries();
                EySpectrum.CalcTimeSeries();
                for (int k = 0; k < ExSpectrum.TsData->size(); ++k)
                {
                        NewData.Ex.data.push_back(ExSpectrum.TsData->at(k));
                        NewData.Ey.data.push_back(EySpectrum.TsData->at(k));
                                
                }
        }       
        NewData.Ex.data.insert(NewData.Ex.data.end(),OriginalData.Hx.data.size() - NewData.Ex.data.size(),0);
        NewData.Ey.data.insert(NewData.Ey.data.end(),OriginalData.Hy.data.size() - NewData.Ey.data.size(),0);
        for (int k = 0; k < NewData.Ex.data.size(); ++k)
        {
                DiffData.Ex.data.push_back(NewData.Ex.data.at(k)-OriginalData.Ex.data.at(k));
                DiffData.Ey.data.push_back(NewData.Ey.data.at(k)-OriginalData.Ey.data.at(k));
        } 
        NewData.WriteData("out.ts");
        DiffData.WriteData("diff.ts");
}

---