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
00160
00161
00162
00163
00164
00165
00166
00167
00168
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 }
