GPLIB++: magnetic_processing.cpp Source File

00001 #include <iostream>
00002 #include "CTsSpectrum.h"
00003 #include "CMtuFormat.h"
00004 #include "CMttData.h"
00005 #include "types.h"
00006 
00007 using namespace std;
00008 int ConstructFilter(tcompdata *FilterCoeff, const string filename, const int seglength, const double freqstep)
00009 {
00010         ifstream filterfile;
00011         ofstream rawout;
00012         ttsdata rawfreqs;
00013         tcompdata rawcoeffs;
00014         double currfreq, currreal, currimag,dummy;
00015         double absresponse,phaseresponse;
00016         
00017         filterfile.open(filename.c_str());
00018         if (!filterfile)
00019         {
00020                 cerr << "Filter-file " << filename << " does not exist ! " << endl;
00021                 return -100;
00022         }
00023         filterfile >> dummy;
00024         while (filterfile.good())
00025         {
00026                 filterfile >> currfreq;
00027                 filterfile.ignore(256,',');
00028                 filterfile >> currreal;
00029                 filterfile.ignore(256,',');
00030                 filterfile >> currimag;
00031                 if (filterfile.good())
00032                 {
00033                         rawfreqs.push_back(currfreq);
00034                         rawcoeffs.push_back(currreal + I * currimag);
00035                 }               
00036         }
00037         filterfile.close();
00038         for (int i = 0; i < seglength/2; ++i)
00039         {
00040                 int j = 0;
00041                 currfreq = i *freqstep;
00042                 while ((j < rawfreqs.size()) && (rawfreqs.at(j) > currfreq))
00043                         ++j;
00044                 --j;
00045                 if (j == rawfreqs.size()-1)
00046                 {
00047                         absresponse = abs(rawcoeffs.at(j)) + (currfreq - rawfreqs.at(j)) * (abs(rawcoeffs.at(j)) - abs(rawcoeffs.at(j-1)))
00048                                 /(rawfreqs.at(j) - rawfreqs.at(j-1));
00049                         phaseresponse = arg(rawcoeffs.at(j)) + (currfreq - rawfreqs.at(j)) * (arg(rawcoeffs.at(j)) - arg(rawcoeffs.at(j-1)))
00050                                 /(rawfreqs.at(j) - rawfreqs.at(j-1));
00051                 }
00052                 else
00053                 {
00054                         absresponse = abs(rawcoeffs.at(j)) + (currfreq - rawfreqs.at(j)) * (abs(rawcoeffs.at(j+1)) - abs(rawcoeffs.at(j)))
00055                                 /(rawfreqs.at(j+1) - rawfreqs.at(j));
00056                         phaseresponse = arg(rawcoeffs.at(j)) + (currfreq - rawfreqs.at(j)) * (arg(rawcoeffs.at(j+1)) - arg(rawcoeffs.at(j)))
00057                                 /(rawfreqs.at(j+1) - rawfreqs.at(j));
00058                 }
00059                 FilterCoeff->push_back(absresponse*(cos(phaseresponse) + I * sin(phaseresponse)));
00060         }
00061         rawout.open((filename+".raw").c_str());
00062         for (int i = 0; i < FilterCoeff->size(); ++i)
00063                 rawout << i * freqstep << " " << FilterCoeff->at(i).real() << " " << FilterCoeff->at(i).imag() << endl; 
00064         rawout.close();
00065         return 0;
00066 }
00067 
00068 int main()
00069 {
00070         CMtuFormat      StationData1;
00071         CMtuFormat      StationData2;
00072         CMttData                MttOutput;
00073         CTsSpectrum HxSpectrum1;
00074         CTsSpectrum HySpectrum1;
00075         CTsSpectrum HxSpectrum2;
00076         CTsSpectrum HySpectrum2;
00077         tcompdata               HxHxCross;
00078         tcompdata               HxHyCross;
00079         tcompdata               HyHxCross;
00080         tcompdata               HyHyCross;
00081         tcompdata       HxHyAuto;
00082         tcompdata               HxHxAuto;
00083         tcompdata               HyHyAuto;
00084         tcompdata       HxFilterCoeff1;
00085         tcompdata       HyFilterCoeff1;
00086         tcompdata       HxFilterCoeff2;
00087         tcompdata       HyFilterCoeff2;
00088         tcompdata               Hdet;
00089         tcompdata               Zxx,Zyy,Zxy,Zyx;
00090         string                  basefilename, hxfiltername1, hyfiltername1, hxfiltername2, hyfiltername2, datafilename1, datafilename2,filename;
00091         int                             seglength = 2400;
00092         int                             nsegs = 0;
00093         int                             dotpos = 0;
00094         double          hxmean1, hymean1, hxmean2, hymean2;
00095         double          freqstep;
00096         
00097         cout << "Filename Station 1: ";
00098         cin >> filename;
00099         datafilename1 = filename;
00100         dotpos = filename.find('.',0);
00101         if (dotpos != string::npos)
00102                 basefilename = filename.erase(dotpos);
00103         hxfiltername1 = basefilename + "_3.cts";
00104         hyfiltername1 = basefilename + "_4.cts";
00105         StationData1.GetData(datafilename1.c_str());
00106         
00107         cout << "Filename Station 2: ";
00108         cin >> filename;
00109         datafilename2 = filename;
00110         dotpos = filename.find('.',0);
00111         if (dotpos != string::npos)
00112                 basefilename = filename.erase(dotpos);
00113         hxfiltername2 = basefilename + "_3.cts";
00114         hyfiltername2 = basefilename + "_4.cts";
00115         StationData2.GetData(datafilename2.c_str());
00116         
00117         HxSpectrum1.TsData = &StationData1.Hx.data;
00118         HySpectrum1.TsData = &StationData1.Hy.data;
00119         HxSpectrum2.TsData = &StationData2.Hx.data;
00120         HySpectrum2.TsData = &StationData2.Hy.data;
00121         
00122         HxSpectrum1.FreqData = new tcompdata;
00123         HySpectrum1.FreqData = new tcompdata;
00124         HxSpectrum2.FreqData = new tcompdata;
00125         HySpectrum2.FreqData = new tcompdata;
00126         
00127         HxSpectrum1.MultiCalc = true;
00128         HySpectrum1.MultiCalc = true;
00129         HxSpectrum2.MultiCalc = true;
00130         HySpectrum2.MultiCalc = true;
00131         
00132         HxHxCross.assign(seglength/2,0);
00133         HxHyCross.assign(seglength/2,0);
00134         HyHxCross.assign(seglength/2,0);
00135         HyHyCross.assign(seglength/2,0);
00136         HxHxAuto.assign(seglength/2,0);
00137         HxHyAuto.assign(seglength/2,0);
00138         HyHyAuto.assign(seglength/2,0);
00139         Hdet.assign(seglength/2,0);
00140         Zxx.assign(seglength/2,0);
00141         Zxy.assign(seglength/2,0);
00142         Zyx.assign(seglength/2,0);
00143         Zyy.assign(seglength/2,0);
00144         
00145         nsegs = HxSpectrum1.TsData->size()/seglength-1;
00146         freqstep = StationData1.samplerate/(seglength);
00147         
00148         ConstructFilter(&HxFilterCoeff1,hxfiltername1,seglength,freqstep);
00149         ConstructFilter(&HyFilterCoeff1,hyfiltername1,seglength,freqstep);
00150         ConstructFilter(&HxFilterCoeff2,hxfiltername2,seglength,freqstep);
00151         ConstructFilter(&HyFilterCoeff2,hyfiltername2,seglength,freqstep);
00152         
00153         for (int j = 0; j < nsegs; ++j)
00154         {
00155                 hxmean1 = 0;
00156                 hymean1 = 0;
00157                 hxmean2 = 0;
00158                 hymean2 = 0;
00159                 /*for (int k = j*seglength; k < (j+1)*seglength; ++k)
00160                 {
00161                         ExSpectrum.TsData->at(k) -= ExSpectrum.TsData->at(j*seglength) 
00162                                 +(ExSpectrum.TsData->at((j+1)*seglength) - ExSpectrum.TsData->at(j*seglength))* (k - j*seglength)/seglength;
00163                         EySpectrum.TsData->at(k) -= EySpectrum.TsData->at(j*seglength) 
00164                                 +(EySpectrum.TsData->at((j+1)*seglength) - EySpectrum.TsData->at(j*seglength))* (k - j*seglength)/seglength;
00165                         HxSpectrum.TsData->at(k) -= HxSpectrum.TsData->at(j*seglength) 
00166                                 +(HxSpectrum.TsData->at((j+1)*seglength) - HxSpectrum.TsData->at(j*seglength))* (k - j*seglength)/seglength;
00167                         HySpectrum.TsData->at(k) -= HySpectrum.TsData->at(j*seglength) 
00168                                 +(HySpectrum.TsData->at((j+1)*seglength) - HySpectrum.TsData->at(j*seglength))* (k - j*seglength)/seglength;
00169                         
00170                 }*/
00171                 for (int k = j*seglength; k < (j+1)*seglength; ++k)
00172                 {
00173                         hxmean1 += HxSpectrum1.TsData->at(k);
00174                         hymean1 += HySpectrum1.TsData->at(k);
00175                         hxmean2 += HxSpectrum2.TsData->at(k);
00176                         hymean2 += HySpectrum2.TsData->at(k);
00177                 }
00178                 for (int k = j*seglength; k < (j+1)*seglength; ++k)
00179                 {
00180                         HxSpectrum1.TsData->at(k) -= hxmean1;
00181                         HySpectrum1.TsData->at(k) -= hymean1;
00182                         HxSpectrum2.TsData->at(k) -= hxmean2;
00183                         HySpectrum2.TsData->at(k) -= hymean2;
00184                         HxSpectrum1.TsData->at(k) *= 0.5 * (1 - cos(2 * PI *(k-seglength) /seglength));
00185                         HySpectrum1.TsData->at(k) *= 0.5 * (1 - cos(2 * PI *(k-seglength) /seglength));
00186                         HxSpectrum2.TsData->at(k) *= 0.5 * (1 - cos(2 * PI *(k-seglength) /seglength));
00187                         HySpectrum2.TsData->at(k) *= 0.5 * (1 - cos(2 * PI *(k-seglength) /seglength));
00188                 }
00189                 HxSpectrum1.CalcSpectrum(j*seglength,(j+1)*seglength);
00190                 HySpectrum1.CalcSpectrum(j*seglength,(j+1)*seglength);
00191                 HxSpectrum2.CalcSpectrum(j*seglength,(j+1)*seglength);
00192                 HySpectrum2.CalcSpectrum(j*seglength,(j+1)*seglength);
00193                 for (int i = 0; i < Zxx.size(); ++i)
00194                 {
00195                         HxSpectrum1.FreqData->at(i) /= HxFilterCoeff1.at(i);
00196                         HySpectrum1.FreqData->at(i) /= HyFilterCoeff1.at(i);
00197                         HxSpectrum2.FreqData->at(i) /= HxFilterCoeff2.at(i);
00198                         HySpectrum2.FreqData->at(i) /= HyFilterCoeff2.at(i);
00199                         HxHxCross.at(i) +=(HxSpectrum1.FreqData->at(i) * conj(HxSpectrum2.FreqData->at(i)));
00200                         HxHyCross.at(i) +=(HxSpectrum1.FreqData->at(i) * conj(HySpectrum2.FreqData->at(i)));
00201                         HyHxCross.at(i) +=(HySpectrum1.FreqData->at(i) * conj(HxSpectrum2.FreqData->at(i)));
00202                         HyHyCross.at(i) +=(HySpectrum1.FreqData->at(i) * conj(HySpectrum2.FreqData->at(i)));
00203                         HxHyAuto.at(i) +=(HxSpectrum1.FreqData->at(i) * conj(HySpectrum1.FreqData->at(i)));
00204                         HxHxAuto.at(i) +=(HxSpectrum1.FreqData->at(i) * conj(HxSpectrum1.FreqData->at(i)));
00205                         HyHyAuto.at(i) +=(HySpectrum1.FreqData->at(i) * conj(HySpectrum1.FreqData->at(i)));
00206                         Hdet.at(i)=(HxHxAuto.at(i) * HyHyAuto.at(i) - HxHyAuto.at(i) * conj(HxHyAuto.at(i)));
00207                         Zxx.at(i) =((HxHxCross.at(i)*HyHyAuto.at(i) - HxHyCross.at(i)*conj(HxHyAuto.at(i)))/Hdet.at(i));
00208                         Zxy.at(i) =((HxHyCross.at(i)*HxHxAuto.at(i) - HxHxCross.at(i)*HxHyAuto.at(i))/Hdet.at(i));
00209                         Zyx.at(i) =((HyHxCross.at(i)*HyHyAuto.at(i) - HyHyCross.at(i)*conj(HxHyAuto.at(i)))/Hdet.at(i));
00210                         Zyy.at(i) =((HyHyCross.at(i)*HxHxAuto.at(i) - HyHxCross.at(i)*HxHyAuto.at(i))/Hdet.at(i));
00211                 }
00212         }       
00213         for (int i = 0; i < Zxx.size(); ++i)
00214         {
00215                 MttOutput.frequency.push_back(i*freqstep);
00216                 MttOutput.DataXX.Z.push_back(Zxx.at(i));
00217                 MttOutput.DataXY.Z.push_back(Zxy.at(i));
00218                 MttOutput.DataYX.Z.push_back(Zyx.at(i));
00219                 MttOutput.DataYY.Z.push_back(Zyy.at(i));
00220         }
00221         MttOutput.frequency.at(0) = 0.00000001;
00222         MttOutput.nfreq = Zxx.size();
00223         MttOutput.DataXX.dZ.assign(Zxx.size(),0);
00224         MttOutput.DataXY.dZ.assign(Zxy.size(),0);
00225         MttOutput.DataYX.dZ.assign(Zyx.size(),0);
00226         MttOutput.DataYY.dZ.assign(Zyy.size(),0);
00227         MttOutput.DataZX.Z.assign(Zyy.size(),0);
00228         MttOutput.DataZY.Z.assign(Zyy.size(),0);
00229         MttOutput.DataZX.dZ.assign(Zyy.size(),0);
00230         MttOutput.DataZY.dZ.assign(Zyy.size(),0);
00231         MttOutput.Rx.assign(Zyy.size(),0);
00232         MttOutput.Ry.assign(Zyy.size(),0);
00233         MttOutput.Rz.assign(Zyy.size(),0);
00234         MttOutput.nu.assign(Zyy.size(),0);
00235         MttOutput.Assign();
00236         MttOutput.WriteData((datafilename1+"m.mtt").c_str());
00237 }