GPLIB++: mtufilter.cpp Source File

00001 #include <iostream>
00002 #include <algorithm>
00003 #include "CTsSpectrum.h"
00004 #include "CTimeSeriesData.h"
00005 #include "CMttData.h"
00006 #include "types.h"
00007 #include "CMtuFilter.h"
00008 
00009 using namespace std;
00010 
00011 int main()
00012 {
00013         CMtuFilter Filter;
00014         CTimeSeriesData         MtuData;
00015         CTsSpectrum ExSpectrum;
00016         CTsSpectrum EySpectrum;
00017         CTsSpectrum HxSpectrum;
00018         CTsSpectrum HySpectrum;
00019         tcompdata       ExFilterCoeff;
00020         tcompdata       EyFilterCoeff;
00021         tcompdata       HxFilterCoeff;
00022         tcompdata       HyFilterCoeff;
00023         tcompdata               Hdet;
00024         ttsdata ExSegment, EySegment, HxSegment, HySegment;
00025         string                  basefilename, exfiltername, eyfiltername, hxfiltername, hyfiltername, datafilename,filename;
00026         int                     seglength = 2400;
00027         int                     nsegs = 0;
00028         int                     dotpos = 0;
00029         double          exmean, eymean, hxmean, hymean;
00030         double          freqstep;
00031         double window;
00032         
00033         cout << "Filename: ";
00034         cin >> filename;
00035         datafilename = filename;
00036         dotpos = filename.find('.',0);
00037         if (dotpos != string::npos)
00038                 basefilename = filename.erase(dotpos);
00039         exfiltername = basefilename + "_1.cts";
00040         eyfiltername = basefilename + "_2.cts";
00041         hxfiltername = basefilename + "_3.cts";
00042         hyfiltername = basefilename + "_4.cts";
00043         MtuData.GetData(datafilename.c_str());
00044         
00045         ExSpectrum.TsData = &MtuData.Data->Ex.data;
00046         EySpectrum.TsData = &MtuData.Data->Ey.data;
00047         HxSpectrum.TsData = &MtuData.Data->Hx.data;
00048         HySpectrum.TsData = &MtuData.Data->Hy.data;
00049         
00050         ExSpectrum.FreqData = new tcompdata;
00051         EySpectrum.FreqData = new tcompdata;
00052         HxSpectrum.FreqData = new tcompdata;
00053         HySpectrum.FreqData = new tcompdata;
00054         
00055         ExSpectrum.MultiCalc = false;
00056         EySpectrum.MultiCalc = false;
00057         HxSpectrum.MultiCalc = false;
00058         HySpectrum.MultiCalc = false;
00059         
00060         nsegs = ExSpectrum.TsData->size()/seglength-1;
00061         freqstep = MtuData.Data->samplerate/(seglength);
00062         
00063         Filter.seglength = seglength;
00064         Filter.freqstep = freqstep;
00065         Filter.FilterCoeff = &ExFilterCoeff;
00066         Filter.GetData(exfiltername);
00067         Filter.WriteData(exfiltername);
00068         Filter.FilterCoeff = &EyFilterCoeff;
00069         Filter.GetData(eyfiltername);
00070         Filter.WriteData(eyfiltername);
00071         Filter.FilterCoeff = &HxFilterCoeff;
00072         Filter.GetData(hxfiltername);
00073         Filter.WriteData(hxfiltername);
00074         Filter.FilterCoeff = &HyFilterCoeff;
00075         Filter.GetData(hyfiltername);
00076         Filter.WriteData(hyfiltername);
00077         
00078         cout << "Filter calculated. " << flush << endl;
00079         for (int j = 0; j < nsegs; ++j)
00080         {
00081                 ExSpectrum.TsData = &MtuData.Data->Ex.data;
00082                 EySpectrum.TsData = &MtuData.Data->Ey.data;
00083                 HxSpectrum.TsData = &MtuData.Data->Hx.data;
00084                 HySpectrum.TsData = &MtuData.Data->Hy.data;
00085                 
00086                 exmean = 0;
00087                 eymean = 0;
00088                 hxmean = 0;
00089                 hymean = 0;
00090                 /*for (int k = j*seglength; k < (j+1)*seglength; ++k)
00091                 {
00092                         ExSpectrum.TsData->at(k) -= ExSpectrum.TsData->at(j*seglength) 
00093                                 +(ExSpectrum.TsData->at((j+1)*seglength) - ExSpectrum.TsData->at(j*seglength))* (k - j*seglength)/seglength;
00094                         EySpectrum.TsData->at(k) -= EySpectrum.TsData->at(j*seglength) 
00095                                 +(EySpectrum.TsData->at((j+1)*seglength) - EySpectrum.TsData->at(j*seglength))* (k - j*seglength)/seglength;
00096                         HxSpectrum.TsData->at(k) -= HxSpectrum.TsData->at(j*seglength) 
00097                                 +(HxSpectrum.TsData->at((j+1)*seglength) - HxSpectrum.TsData->at(j*seglength))* (k - j*seglength)/seglength;
00098                         HySpectrum.TsData->at(k) -= HySpectrum.TsData->at(j*seglength) 
00099                                 +(HySpectrum.TsData->at((j+1)*seglength) - HySpectrum.TsData->at(j*seglength))* (k - j*seglength)/seglength;
00100                         
00101                 }*/
00102                 for (int k = j*seglength; k < (j+1)*seglength; ++k)
00103                 {
00104                         exmean += ExSpectrum.TsData->at(k);
00105                         eymean += EySpectrum.TsData->at(k);
00106                         hxmean += HxSpectrum.TsData->at(k);
00107                         hymean += HySpectrum.TsData->at(k);
00108                 }
00109                 exmean /= seglength;
00110                 eymean /= seglength;
00111                 hxmean /= seglength;
00112                 hymean /= seglength;
00113                 for (int k = j*seglength; k < (j+1)*seglength; ++k)
00114                 {
00115                         ExSpectrum.TsData->at(k) -= exmean;
00116                         EySpectrum.TsData->at(k) -= eymean;
00117                         HxSpectrum.TsData->at(k) -= hxmean;
00118                         HySpectrum.TsData->at(k) -= hymean;
00119                         /*window = 0.54 - 0.46 * cos(2 * PI *(k-seglength) /seglength);
00120                         ExSpectrum.TsData->at(k) *= window;
00121                         EySpectrum.TsData->at(k) *= window;
00122                         HxSpectrum.TsData->at(k) *= window;
00123                         HySpectrum.TsData->at(k) *= window;*/
00124                 }
00125                 ExSpectrum.CalcSpectrum(j*seglength,(j+1)*seglength);
00126                 EySpectrum.CalcSpectrum(j*seglength,(j+1)*seglength);
00127                 HxSpectrum.CalcSpectrum(j*seglength,(j+1)*seglength);
00128                 HySpectrum.CalcSpectrum(j*seglength,(j+1)*seglength);
00129                 //cout << "Spectra calculated. " << flush << endl;
00130                 for (int i = 0; i < ExFilterCoeff.size(); ++i)
00131                 {
00132                         ExSpectrum.FreqData->at(i) /= ExFilterCoeff.at(i);
00133                         EySpectrum.FreqData->at(i) /= EyFilterCoeff.at(i);
00134                         HxSpectrum.FreqData->at(i) /= HxFilterCoeff.at(i);
00135                         HySpectrum.FreqData->at(i) /= HyFilterCoeff.at(i);
00136                 }
00137                 //cout << "Filter applied. " << flush << endl;
00138                 ExSpectrum.TsData = &ExSegment;
00139                 EySpectrum.TsData = &EySegment;
00140                 HxSpectrum.TsData = &HxSegment;
00141                 HySpectrum.TsData = &HySegment;
00142                 
00143                 ExSpectrum.CalcTimeSeries();
00144                 EySpectrum.CalcTimeSeries();
00145                 HxSpectrum.CalcTimeSeries();
00146                 HySpectrum.CalcTimeSeries();
00147                 //cout << "TimeSeries calculated. " << flush << endl;
00148                 
00149                 
00150                 /*for (int k = 0; k < seglength; ++k)
00151                 {
00152                         window = 0.54 - 0.46 * cos(2 * PI *(k-seglength) /seglength);
00153                         ExSpectrum.TsData->at(k) /= window;
00154                         EySpectrum.TsData->at(k) /= window;
00155                         HxSpectrum.TsData->at(k) /= window;
00156                         HySpectrum.TsData->at(k) /= window;
00157                         //ExSpectrum.TsData->at(k) += exmean;
00158                         //EySpectrum.TsData->at(k) += eymean;
00159                         //HxSpectrum.TsData->at(k) += hxmean;
00160                         //HySpectrum.TsData->at(k) += hymean;
00161                 }*/
00162                 copy(ExSpectrum.TsData->begin(),ExSpectrum.TsData->end(),MtuData.Data->Ex.data.begin()+j*seglength);
00163                 copy(EySpectrum.TsData->begin(),EySpectrum.TsData->end(),MtuData.Data->Ey.data.begin()+j*seglength);
00164                 copy(HxSpectrum.TsData->begin(),HxSpectrum.TsData->end(),MtuData.Data->Hx.data.begin()+j*seglength);
00165                 copy(HySpectrum.TsData->begin(),HySpectrum.TsData->end(),MtuData.Data->Hy.data.begin()+j*seglength);
00166                 //cout << "Corrections applied. " << flush << endl;
00167         }       
00168         MtuData.Data->Ex.data.resize(nsegs*seglength);
00169         MtuData.Data->Ey.data.resize(nsegs*seglength);
00170         MtuData.Data->Hx.data.resize(nsegs*seglength);
00171         MtuData.Data->Hy.data.resize(nsegs*seglength);
00172         MtuData.WriteAsBirrp((datafilename+"_filt").c_str());
00173 }
Simple Processing/mtufilter.cpp
