simple_processing2.cpp

Go to the documentation of this file.
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     /*if (MtuData.GetData().GetSamplerate() == 1)
00062      {
00063      double rate;
00064      cout << "Enter sample rate in Hz: ";
00065      cin >> rate;
00066      MtuData.GetData().SetSamplerate(rate);
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; // write mean to file as new dataset
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; //write median to file as new dataset
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   }