00001 #include <iostream>
00002 #include <iomanip>
00003 #include <fstream>
00004 #include "TsSpectrum.h"
00005 #include "TimeSeriesData.h"
00006 #include "statutils.h"
00007 #include "MTStation.h"
00008 #include "types.h"
00009 #include "WFunc.h"
00010 #include "TimeFrequency.h"
00011 #include "VecMat.h"
00012 #include <boost/bind.hpp>
00013 #include <complex>
00014 #include "MtuFilter.h"
00015 #include "convert.h"
00016 #include "netcdfcpp.h"
00017 #include "Util.h"
00018
00019 using namespace std;
00020 using namespace gplib;
00021
00022 string version =
00023 "$Id: simple_processing2.cpp 1816 2009-09-07 11:28:35Z mmoorkamp $";
00024 int main(int argc, char *argv[])
00025 {
00026 cout
00027 << "This is simpleproc2: Simple most estimation of MT transfer function, stacking after Z calculation"
00028 << endl << endl;
00029 cout << " Usage simpleproc inputfilename";
00030 cout
00031 << " Output will have the same name as MTU Input with '.mtt' appended "
00032 << endl << endl;
00033 cout << " This is Version: " << version << endl << endl;
00034
00035 TimeSeriesData MtuData;
00036 MTStation MeanOutput, MedianOutput;
00037 string filename, basefilename;
00038 const unsigned int seglength = 120;
00039 const unsigned int navgwindows = 5;
00040
00041 if (argc == 2)
00042 {
00043 filename = argv[1];
00044 }
00045 else
00046 {
00047 filename = AskFilename("Time Series Filename: ");
00048 }
00049 string datafilename = filename;
00050 size_t dotpos = filename.find('.', 0);
00051 if (dotpos != string::npos)
00052 basefilename = filename.erase(dotpos);
00053 string exfiltername = basefilename + "_1.cts";
00054 string eyfiltername = basefilename + "_2.cts";
00055 string hxfiltername = basefilename + "_3.cts";
00056 string hyfiltername = basefilename + "_4.cts";
00057 MtuData.GetData(datafilename.c_str());
00058
00059 const unsigned int nfreqs = seglength / 2 + 1;
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069 double freqstep = MtuData.GetData().GetSamplerate() / (seglength);
00070 MtuFilter ExFilter(seglength, freqstep), EyFilter(seglength, freqstep),
00071 HxFilter(seglength, freqstep), HyFilter(seglength, freqstep);
00072
00073 ExFilter.GetData(exfiltername);
00074 EyFilter.GetData(eyfiltername);
00075 HxFilter.GetData(hxfiltername);
00076 HyFilter.GetData(hyfiltername);
00077
00078 gplib::cmat ExTimeFrequency, EyTimeFrequency, HxTimeFrequency,
00079 HyTimeFrequency;
00080
00081 ExTimeFrequency = TimeFrequency(
00082 MtuData.GetData().GetEx().GetData().begin(),
00083 MtuData.GetData().GetEx().GetData().end(), seglength, Hamming());
00084 EyTimeFrequency = TimeFrequency(
00085 MtuData.GetData().GetEy().GetData().begin(),
00086 MtuData.GetData().GetEy().GetData().end(), seglength, Hamming());
00087 HxTimeFrequency = TimeFrequency(
00088 MtuData.GetData().GetHx().GetData().begin(),
00089 MtuData.GetData().GetHx().GetData().end(), seglength, Hamming());
00090 HyTimeFrequency = TimeFrequency(
00091 MtuData.GetData().GetHy().GetData().begin(),
00092 MtuData.GetData().GetHy().GetData().end(), seglength, Hamming());
00093
00094 const unsigned int nsegs = ExTimeFrequency.size1();
00095
00096 gplib::cmat ExHxCorr(nsegs, nfreqs);
00097 gplib::cmat ExHyCorr(nsegs, nfreqs);
00098 gplib::cmat EyHxCorr(nsegs, nfreqs);
00099 gplib::cmat EyHyCorr(nsegs, nfreqs);
00100 gplib::cmat HxHxCorr(nsegs, nfreqs);
00101 gplib::cmat HxHyCorr(nsegs, nfreqs);
00102 gplib::cmat HyHyCorr(nsegs, nfreqs);
00103 gplib::cmat Hdet(nsegs, nfreqs);
00104 gplib::cmat Zxx(nsegs, nfreqs);
00105 gplib::cmat Zxy(nsegs, nfreqs);
00106 gplib::cmat Zyx(nsegs, nfreqs);
00107 gplib::cmat Zyy(nsegs, nfreqs);
00108 gplib::rmat Weights(nsegs, nfreqs);
00109 ofstream logfile((datafilename + ".log").c_str());
00110 MeanOutput.AssignAll(nfreqs);
00111 MedianOutput.AssignAll(nfreqs);
00112
00113 NcFile combrescdf((datafilename + ".zxyarg.nc").c_str(), NcFile::Replace);
00114 NcDim* xd = combrescdf.add_dim("x", nfreqs);
00115 NcDim* yd = combrescdf.add_dim("y", nsegs);
00116 NcVar* x = combrescdf.add_var("x", ncFloat, xd);
00117 NcVar* y = combrescdf.add_var("y", ncFloat, yd);
00118 NcVar* z = combrescdf.add_var("z", ncFloat, xd, yd);
00119 float *xvals = new float[xd->size()];
00120 float *yvals = new float[yd->size()];
00121 float *zvals = new float[xd->size() * yd->size()];
00122
00123 for (size_t i = 0; i < nfreqs; ++i)
00124 {
00125 ofstream zxyfile(
00126 (datafilename + "_zxy_" + stringify(i * freqstep)).c_str());
00127 ofstream zxyabsfile((datafilename + "_zxyabs_"
00128 + stringify(i * freqstep)).c_str());
00129 ofstream zxyargfile((datafilename + "_zxyarg_"
00130 + stringify(i * freqstep)).c_str());
00131 ofstream hdetfile(
00132 (datafilename + "_hdet_" + stringify(i * freqstep)).c_str());
00133 for (size_t j = 0; j < nsegs; ++j)
00134 {
00135 ExTimeFrequency(j, i) /= ExFilter.GetFilterCoeff().at(i);
00136 EyTimeFrequency(j, i) /= EyFilter.GetFilterCoeff().at(i);
00137 HxTimeFrequency(j, i) /= HxFilter.GetFilterCoeff().at(i);
00138 HyTimeFrequency(j, i) /= HyFilter.GetFilterCoeff().at(i);
00139 ExHxCorr(j, i) = (ExTimeFrequency(j, i) * conj(
00140 HxTimeFrequency(j, i)));
00141 ExHyCorr(j, i) = (ExTimeFrequency(j, i) * conj(
00142 HyTimeFrequency(j, i)));
00143 EyHxCorr(j, i) = (EyTimeFrequency(j, i) * conj(
00144 HxTimeFrequency(j, i)));
00145 EyHyCorr(j, i) = (EyTimeFrequency(j, i) * conj(
00146 HyTimeFrequency(j, i)));
00147 HxHyCorr(j, i) = (HxTimeFrequency(j, i) * conj(
00148 HyTimeFrequency(j, i)));
00149 HxHxCorr(j, i) = (HxTimeFrequency(j, i) * conj(
00150 HxTimeFrequency(j, i)));
00151 HyHyCorr(j, i) = (HyTimeFrequency(j, i) * conj(
00152 HyTimeFrequency(j, i)));
00153 }
00154 for (size_t j = 0; j < nsegs - navgwindows; ++j)
00155 {
00156 complex<double> hxhxavg = accumulate(column(HxHxCorr, i).begin()
00157 + j, column(HxHxCorr, i).begin() + j + navgwindows, complex<
00158 double> ());
00159 complex<double> exhxavg = accumulate(column(ExHxCorr, i).begin()
00160 + j, column(ExHxCorr, i).begin() + j + navgwindows, complex<
00161 double> ());
00162 complex<double> exhyavg = accumulate(column(ExHyCorr, i).begin()
00163 + j, column(ExHyCorr, i).begin() + j + navgwindows, complex<
00164 double> ());
00165 complex<double> eyhxavg = accumulate(column(EyHxCorr, i).begin()
00166 + j, column(EyHxCorr, i).begin() + j + navgwindows, complex<
00167 double> ());
00168 complex<double> eyhyavg = accumulate(column(EyHyCorr, i).begin()
00169 + j, column(EyHyCorr, i).begin() + j + navgwindows, complex<
00170 double> ());
00171 complex<double> hxhyavg = accumulate(column(HxHyCorr, i).begin()
00172 + j, column(HxHyCorr, i).begin() + j + navgwindows, complex<
00173 double> ());
00174 complex<double> hyhyavg = accumulate(column(HyHyCorr, i).begin()
00175 + j, column(HyHyCorr, i).begin() + j + navgwindows, complex<
00176 double> ());
00177 Hdet(j, i) = (hxhxavg * hyhyavg - hxhyavg * conj(hxhyavg));
00178 hdetfile << setfill(' ') << setw(15) << setprecision(9) << abs(
00179 Hdet(j, i).real()) << endl;
00180
00181 if (abs(Hdet(j, i)) > 0.0)
00182 {
00183
00184 Zxx(j, i) = ((exhxavg * hyhyavg - exhyavg * conj(hxhyavg))
00185 / Hdet(j, i));
00186 Zxy(j, i) = ((exhyavg * hxhxavg - exhxavg * hxhyavg) / Hdet(j,
00187 i));
00188 Zyx(j, i) = ((eyhxavg * hyhyavg - eyhyavg * conj(hxhyavg))
00189 / Hdet(j, i));
00190 Zyy(j, i) = ((eyhyavg * hxhxavg - eyhxavg * hxhyavg) / Hdet(j,
00191 i));
00192 zvals[i * nsegs + j] = arg(Zxy(j, i));
00193 Weights(j, i) = 1.0;
00194 zxyfile << setfill(' ') << setw(15) << setprecision(9) << Zxy(
00195 j, i).real() << " ";
00196 zxyfile << setfill(' ') << setw(15) << setprecision(9) << Zxy(
00197 j, i).imag() << endl;
00198 }
00199 else
00200 {
00201 zvals[i * nsegs + j] = -100;
00202 Zxx(j, i) = 0.0;
00203 Zxy(j, i) = 0.0;
00204 Zyx(j, i) = 0.0;
00205 Zyy(j, i) = 0.0;
00206 Weights(j, i) = 0.0;
00207 }
00208 zxyabsfile << setfill(' ') << setw(15) << setprecision(9) << abs(
00209 Zxy(j, i)) << endl;
00210 zxyargfile << setfill(' ') << setw(15) << setprecision(9) << arg(
00211 Zxy(j, i)) << endl;
00212 }
00213
00214 ofstream weightfile((datafilename + ".weights").c_str());
00215 double currfrequency = i * freqstep;
00216 xvals[i] = currfrequency;
00217 complex<double> currzxx, currzxy, currzyx, currzyy;
00218 const int goodcount = sum(column(Weights, i));
00219 double normalize = 1. / goodcount;
00220
00221 cout << "Freq: " << currfrequency << " Good: " << goodcount
00222 << " Total: " << Weights.size1() << endl;
00223 weightfile << sum(column(Weights, i)) << endl;
00224 currzxx = prec_inner_prod(column(Zxx, i), column(Weights, i));
00225 currzxx *= normalize;
00226 currzxy = prec_inner_prod(column(Zxy, i), column(Weights, i));
00227 currzxy *= normalize;
00228 currzyx = prec_inner_prod(column(Zyx, i), column(Weights, i));
00229 currzyx *= normalize;
00230 currzyy = prec_inner_prod(column(Zyy, i), column(Weights, i));
00231 currzyy *= normalize;
00232 MTTensor CurrZ(currzxx, currzxy, currzyx, currzyy, currfrequency);
00233 MeanOutput.SetMTData().at(i) = CurrZ;
00234 zxyfile << endl << endl;
00235 zxyfile << setfill(' ') << setw(15) << setprecision(9)
00236 << currzxy.real() << " ";
00237 zxyfile << setfill(' ') << setw(15) << setprecision(9)
00238 << currzxy.imag() << endl;
00239 zxyfile << endl << endl;
00240 vector<double> zxyreal(goodcount, 0.0);
00241 vector<double> zxyimag(goodcount, 0.0);
00242 int currindex = 0;
00243 for (int j = 0; j < nsegs; ++j)
00244 {
00245 yvals[j] = j;
00246 if (Weights(j, i) > 0.0)
00247 {
00248 zxyreal.at(currindex) = Zxy(j, i).real();
00249 zxyimag.at(currindex) = Zxy(j, i).imag();
00250 ++currindex;
00251 }
00252 }
00253 long double zxyrmedian = Median(zxyreal.begin(), zxyreal.end());
00254 long double zxyimedian = Median(zxyimag.begin(), zxyimag.end());
00255 zxyfile << setfill(' ') << setw(15) << setprecision(9) << zxyrmedian
00256 << " ";
00257 zxyfile << setfill(' ') << setw(15) << setprecision(9) << zxyimedian
00258 << endl;
00259 currzxy.real() = zxyrmedian;
00260 currzxy.imag() = zxyimedian;
00261 MTTensor NewZ(currzxx, currzxy, currzyx, currzyy, currfrequency);
00262 MedianOutput.SetMTData().at(i) = NewZ;
00263 }
00264 x->put(xvals, xd->size());
00265 y->put(yvals, yd->size());
00266 z->put(zvals, z->edges());
00267 MeanOutput.WriteAsMtt((datafilename + ".mean").c_str());
00268 MedianOutput.WriteAsMtt((datafilename + "median").c_str());
00269 }