|
GPLIB++
|
00001 #include <iostream> 00002 #include <iomanip> 00003 #include <fstream> 00004 #include <boost/program_options.hpp> 00005 #include "TsSpectrum.h" 00006 #include "TimeSeriesData.h" 00007 #include "statutils.h" 00008 #include "MTStation.h" 00009 #include "types.h" 00010 #include "WFunc.h" 00011 #include "TimeFrequency.h" 00012 #include "VecMat.h" 00013 #include <boost/bind.hpp> 00014 #include <complex> 00015 #include "MtuFilter.h" 00016 #include "convert.h" 00017 #include "netcdfcpp.h" 00018 #include "Util.h" 00019 00020 using namespace std; 00021 using namespace gplib; 00022 namespace po = boost::program_options; 00023 00024 /*! 00025 * \addtogroup UtilProgs Utility Programs 00026 *@{ 00027 * \file simple_processing2.cpp 00028 * Very simple MT processing code based on least squares estimation as described by Sims 1971 similar 00029 * to simple_processing.cpp. Here we calculate the impedance for a small number of sections to create 00030 * a plot of impedance estimates for different times. This can help to identify problems with the data. 00031 * So far no error calculation. 00032 */ 00033 00034 00035 string version = 00036 "$Id: simple_processing2.cpp 1852 2010-05-20 09:14:53Z mmoorkamp $"; 00037 int main(int argc, char *argv[]) 00038 { 00039 cout 00040 << "This is simpleproc2: Simple most estimation of MT transfer function, stacking after Z calculation" 00041 << endl << endl; 00042 cout << " Usage simpleproc inputfilename"; 00043 cout 00044 << " Output will have the same name as MTU Input with '.mtt' appended " 00045 << endl << endl; 00046 cout << " This is Version: " << version << endl << endl; 00047 00048 TimeSeriesData MtuData; 00049 MTStation MeanOutput, MedianOutput; 00050 string infilename, basefilename; 00051 00052 unsigned int navgwindows = 5; 00053 unsigned int seglength = 2400; 00054 bool mtufilter = false; 00055 double rate = 1.0; 00056 po::options_description desc("Allowed options"); 00057 desc.add_options()("help", "produce help message")("seglength", po::value< 00058 unsigned int>(&seglength)->default_value(2400), 00059 "The length of an individual segment for spectral calculations")( 00060 "rate", po::value(&rate)->default_value(1.0), 00061 "Set the sampling rate in Hz")("mtufilter", 00062 po::value<bool>(&mtufilter)->default_value(false), 00063 "Read in filter information for mtu data")("input-file", po::value< 00064 string>(&infilename), "input file"); 00065 00066 po::positional_options_description p; 00067 p.add("input-file", -1); 00068 00069 po::variables_map vm; 00070 po::store( 00071 po::command_line_parser(argc, argv). options(desc).positional(p).run(), 00072 vm); 00073 po::notify(vm); 00074 00075 if (vm.count("help")) 00076 { 00077 cout << desc << "\n"; 00078 return 1; 00079 } 00080 00081 00082 string datafilename = infilename; 00083 size_t dotpos = infilename.find('.', 0); 00084 if (dotpos != string::npos) 00085 basefilename = infilename.erase(dotpos); 00086 string exfiltername = basefilename + "_1.cts"; 00087 string eyfiltername = basefilename + "_2.cts"; 00088 string hxfiltername = basefilename + "_3.cts"; 00089 string hyfiltername = basefilename + "_4.cts"; 00090 MtuData.GetData(datafilename.c_str()); 00091 00092 const unsigned int nfreqs = seglength / 2 + 1; 00093 00094 /*if (MtuData.GetData().GetSamplerate() == 1) 00095 { 00096 double rate; 00097 cout << "Enter sample rate in Hz: "; 00098 cin >> rate; 00099 MtuData.GetData().SetSamplerate(rate); 00100 }*/ 00101 00102 double freqstep = MtuData.GetData().GetSamplerate() / (seglength); 00103 MtuFilter ExFilter(seglength, freqstep), EyFilter(seglength, freqstep), 00104 HxFilter(seglength, freqstep), HyFilter(seglength, freqstep); 00105 00106 if (mtufilter) 00107 { 00108 ExFilter.GetData(exfiltername); 00109 EyFilter.GetData(eyfiltername); 00110 HxFilter.GetData(hxfiltername); 00111 HyFilter.GetData(hyfiltername); 00112 } 00113 gplib::cmat ExTimeFrequency, EyTimeFrequency, HxTimeFrequency, 00114 HyTimeFrequency; 00115 00116 ExTimeFrequency = TimeFrequency( 00117 MtuData.GetData().GetEx().GetData().begin(), 00118 MtuData.GetData().GetEx().GetData().end(), seglength, Hamming()); 00119 EyTimeFrequency = TimeFrequency( 00120 MtuData.GetData().GetEy().GetData().begin(), 00121 MtuData.GetData().GetEy().GetData().end(), seglength, Hamming()); 00122 HxTimeFrequency = TimeFrequency( 00123 MtuData.GetData().GetHx().GetData().begin(), 00124 MtuData.GetData().GetHx().GetData().end(), seglength, Hamming()); 00125 HyTimeFrequency = TimeFrequency( 00126 MtuData.GetData().GetHy().GetData().begin(), 00127 MtuData.GetData().GetHy().GetData().end(), seglength, Hamming()); 00128 00129 const unsigned int nsegs = ExTimeFrequency.size1(); 00130 00131 gplib::cmat ExHxCorr(nsegs, nfreqs); 00132 gplib::cmat ExHyCorr(nsegs, nfreqs); 00133 gplib::cmat EyHxCorr(nsegs, nfreqs); 00134 gplib::cmat EyHyCorr(nsegs, nfreqs); 00135 gplib::cmat HxHxCorr(nsegs, nfreqs); 00136 gplib::cmat HxHyCorr(nsegs, nfreqs); 00137 gplib::cmat HyHyCorr(nsegs, nfreqs); 00138 gplib::cmat Hdet(nsegs, nfreqs); 00139 gplib::cmat Zxx(nsegs, nfreqs); 00140 gplib::cmat Zxy(nsegs, nfreqs); 00141 gplib::cmat Zyx(nsegs, nfreqs); 00142 gplib::cmat Zyy(nsegs, nfreqs); 00143 gplib::rmat Weights(nsegs, nfreqs); 00144 ofstream logfile((datafilename + ".log").c_str()); 00145 MeanOutput.AssignAll(nfreqs); 00146 MedianOutput.AssignAll(nfreqs); 00147 00148 NcFile combrescdf((datafilename + ".zxyarg.nc").c_str(), NcFile::Replace); 00149 NcDim* xd = combrescdf.add_dim("x", nfreqs); 00150 NcDim* yd = combrescdf.add_dim("y", nsegs); 00151 NcVar* x = combrescdf.add_var("x", ncFloat, xd); 00152 NcVar* y = combrescdf.add_var("y", ncFloat, yd); 00153 NcVar* z = combrescdf.add_var("z", ncFloat, xd, yd); 00154 float *xvals = new float[xd->size()]; 00155 float *yvals = new float[yd->size()]; 00156 float *zvals = new float[xd->size() * yd->size()]; 00157 00158 for (size_t i = 0; i < nfreqs; ++i) 00159 { 00160 ofstream zxyfile( 00161 (datafilename + "_zxy_" + stringify(i * freqstep)).c_str()); 00162 ofstream zxyabsfile((datafilename + "_zxyabs_" 00163 + stringify(i * freqstep)).c_str()); 00164 ofstream zxyargfile((datafilename + "_zxyarg_" 00165 + stringify(i * freqstep)).c_str()); 00166 ofstream hdetfile( 00167 (datafilename + "_hdet_" + stringify(i * freqstep)).c_str()); 00168 for (size_t j = 0; j < nsegs; ++j) 00169 { 00170 if (mtufilter) 00171 { 00172 ExTimeFrequency(j, i) /= ExFilter.GetFilterCoeff().at(i); 00173 EyTimeFrequency(j, i) /= EyFilter.GetFilterCoeff().at(i); 00174 HxTimeFrequency(j, i) /= HxFilter.GetFilterCoeff().at(i); 00175 HyTimeFrequency(j, i) /= HyFilter.GetFilterCoeff().at(i); 00176 } 00177 ExHxCorr(j, i) = (ExTimeFrequency(j, i) * conj( 00178 HxTimeFrequency(j, i))); 00179 ExHyCorr(j, i) = (ExTimeFrequency(j, i) * conj( 00180 HyTimeFrequency(j, i))); 00181 EyHxCorr(j, i) = (EyTimeFrequency(j, i) * conj( 00182 HxTimeFrequency(j, i))); 00183 EyHyCorr(j, i) = (EyTimeFrequency(j, i) * conj( 00184 HyTimeFrequency(j, i))); 00185 HxHyCorr(j, i) = (HxTimeFrequency(j, i) * conj( 00186 HyTimeFrequency(j, i))); 00187 HxHxCorr(j, i) = (HxTimeFrequency(j, i) * conj( 00188 HxTimeFrequency(j, i))); 00189 HyHyCorr(j, i) = (HyTimeFrequency(j, i) * conj( 00190 HyTimeFrequency(j, i))); 00191 } 00192 for (size_t j = 0; j < nsegs - navgwindows; ++j) 00193 { 00194 complex<double> hxhxavg = accumulate(column(HxHxCorr, i).begin() 00195 + j, column(HxHxCorr, i).begin() + j + navgwindows, complex< 00196 double> ()); 00197 complex<double> exhxavg = accumulate(column(ExHxCorr, i).begin() 00198 + j, column(ExHxCorr, i).begin() + j + navgwindows, complex< 00199 double> ()); 00200 complex<double> exhyavg = accumulate(column(ExHyCorr, i).begin() 00201 + j, column(ExHyCorr, i).begin() + j + navgwindows, complex< 00202 double> ()); 00203 complex<double> eyhxavg = accumulate(column(EyHxCorr, i).begin() 00204 + j, column(EyHxCorr, i).begin() + j + navgwindows, complex< 00205 double> ()); 00206 complex<double> eyhyavg = accumulate(column(EyHyCorr, i).begin() 00207 + j, column(EyHyCorr, i).begin() + j + navgwindows, complex< 00208 double> ()); 00209 complex<double> hxhyavg = accumulate(column(HxHyCorr, i).begin() 00210 + j, column(HxHyCorr, i).begin() + j + navgwindows, complex< 00211 double> ()); 00212 complex<double> hyhyavg = accumulate(column(HyHyCorr, i).begin() 00213 + j, column(HyHyCorr, i).begin() + j + navgwindows, complex< 00214 double> ()); 00215 Hdet(j, i) = (hxhxavg * hyhyavg - hxhyavg * conj(hxhyavg)); 00216 hdetfile << setfill(' ') << setw(15) << setprecision(9) << abs( 00217 Hdet(j, i).real()) << endl; 00218 00219 if (abs(Hdet(j, i)) > 0.0) 00220 { 00221 00222 Zxx(j, i) = ((exhxavg * hyhyavg - exhyavg * conj(hxhyavg)) 00223 / Hdet(j, i)); 00224 Zxy(j, i) = ((exhyavg * hxhxavg - exhxavg * hxhyavg) / Hdet(j, 00225 i)); 00226 Zyx(j, i) = ((eyhxavg * hyhyavg - eyhyavg * conj(hxhyavg)) 00227 / Hdet(j, i)); 00228 Zyy(j, i) = ((eyhyavg * hxhxavg - eyhxavg * hxhyavg) / Hdet(j, 00229 i)); 00230 zvals[i * nsegs + j] = arg(Zxy(j, i)); 00231 Weights(j, i) = 1.0; 00232 zxyfile << setfill(' ') << setw(15) << setprecision(9) << Zxy( 00233 j, i).real() << " "; 00234 zxyfile << setfill(' ') << setw(15) << setprecision(9) << Zxy( 00235 j, i).imag() << endl; 00236 } 00237 else 00238 { 00239 zvals[i * nsegs + j] = -100; 00240 Zxx(j, i) = 0.0; 00241 Zxy(j, i) = 0.0; 00242 Zyx(j, i) = 0.0; 00243 Zyy(j, i) = 0.0; 00244 Weights(j, i) = 0.0; 00245 } 00246 zxyabsfile << setfill(' ') << setw(15) << setprecision(9) << abs( 00247 Zxy(j, i)) << endl; 00248 zxyargfile << setfill(' ') << setw(15) << setprecision(9) << arg( 00249 Zxy(j, i)) << endl; 00250 } 00251 00252 ofstream weightfile((datafilename + ".weights").c_str()); 00253 double currfrequency = i * freqstep; 00254 xvals[i] = currfrequency; 00255 complex<double> currzxx, currzxy, currzyx, currzyy; 00256 const int goodcount = sum(column(Weights, i)); 00257 double normalize = 1. / goodcount; 00258 00259 cout << "Freq: " << currfrequency << " Good: " << goodcount 00260 << " Total: " << Weights.size1() << endl; 00261 weightfile << sum(column(Weights, i)) << endl; 00262 currzxx = prec_inner_prod(column(Zxx, i), column(Weights, i)); 00263 currzxx *= normalize; 00264 currzxy = prec_inner_prod(column(Zxy, i), column(Weights, i)); 00265 currzxy *= normalize; 00266 currzyx = prec_inner_prod(column(Zyx, i), column(Weights, i)); 00267 currzyx *= normalize; 00268 currzyy = prec_inner_prod(column(Zyy, i), column(Weights, i)); 00269 currzyy *= normalize; 00270 MTTensor CurrZ(currzxx, currzxy, currzyx, currzyy, currfrequency); 00271 MeanOutput.SetMTData().at(i) = CurrZ; 00272 zxyfile << endl << endl; // write mean to file as new dataset 00273 zxyfile << setfill(' ') << setw(15) << setprecision(9) 00274 << currzxy.real() << " "; 00275 zxyfile << setfill(' ') << setw(15) << setprecision(9) 00276 << currzxy.imag() << endl; 00277 zxyfile << endl << endl; //write median to file as new dataset 00278 vector<double> zxyreal(goodcount, 0.0); 00279 vector<double> zxyimag(goodcount, 0.0); 00280 int currindex = 0; 00281 for (unsigned int j = 0; j < nsegs; ++j) 00282 { 00283 yvals[j] = j; 00284 if (Weights(j, i) > 0.0) 00285 { 00286 zxyreal.at(currindex) = Zxy(j, i).real(); 00287 zxyimag.at(currindex) = Zxy(j, i).imag(); 00288 ++currindex; 00289 } 00290 } 00291 long double zxyrmedian = Median(zxyreal.begin(), zxyreal.end()); 00292 long double zxyimedian = Median(zxyimag.begin(), zxyimag.end()); 00293 zxyfile << setfill(' ') << setw(15) << setprecision(9) << zxyrmedian 00294 << " "; 00295 zxyfile << setfill(' ') << setw(15) << setprecision(9) << zxyimedian 00296 << endl; 00297 currzxy = std::complex<double>(zxyrmedian, zxyimedian); 00298 MTTensor NewZ(currzxx, currzxy, currzyx, currzyy, currfrequency); 00299 MedianOutput.SetMTData().at(i) = NewZ; 00300 } 00301 x->put(xvals, xd->size()); 00302 y->put(yvals, yd->size()); 00303 z->put(zvals, z->edges()); 00304 MeanOutput.WriteAsMtt((datafilename + ".mean").c_str()); 00305 MedianOutput.WriteAsMtt((datafilename + "median").c_str()); 00306 } 00307 /*@}*/
1.7.6.1