GPLIB++: simple_processing2.cpp Source 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 string version = "$Id: simple_processing2.cpp 1704 2008-06-02 14:50:05Z mmoorkamp $";
00021 int main(int argc, char *argv[])
00022 {
00023         cout << "This is simpleproc2: Simple most estimation of MT transfer function, stacking after Z calculation" << endl<< endl;
00024         cout << " Usage   simpleproc  inputfilename";
00025         cout << " Output will have the same name as MTU Input with '.mtt' appended " << endl<< endl;
00026         cout << " This is Version: " << version << endl << endl;
00027         
00028         TimeSeriesData  MtuData;
00029         MTStation       MeanOutput, MedianOutput;
00030         string filename, basefilename;
00031         const unsigned int      seglength = 120;
00032         const unsigned int navgwindows = 5;
00033         
00034         
00035         if (argc == 2)
00036         {
00037                 filename = argv[1];
00038         }
00039         else
00040         {
00041                 filename = AskFilename("Time Series Filename: ");
00042         }
00043         string datafilename = filename;
00044         size_t dotpos = filename.find('.',0);
00045         if (dotpos != string::npos)
00046                 basefilename = filename.erase(dotpos);
00047         string exfiltername = basefilename + "_1.cts";
00048         string eyfiltername = basefilename + "_2.cts";
00049         string hxfiltername = basefilename + "_3.cts";
00050         string hyfiltername = basefilename + "_4.cts";
00051         MtuData.GetData(datafilename.c_str());
00052         
00053         const unsigned int nfreqs = seglength/2+1;
00054         
00055         
00056         /*if (MtuData.GetData().GetSamplerate() == 1)
00057         {
00058                 double rate;
00059                 cout << "Enter sample rate in Hz: ";
00060                 cin >> rate;
00061                 MtuData.GetData().SetSamplerate(rate);
00062         }*/
00063         
00064         double freqstep = MtuData.GetData().GetSamplerate()/(seglength);
00065         MtuFilter ExFilter(seglength,freqstep),EyFilter(seglength,freqstep),HxFilter(seglength,freqstep),HyFilter(seglength,freqstep);
00066         
00067         ExFilter.GetData(exfiltername);
00068         EyFilter.GetData(eyfiltername);
00069         HxFilter.GetData(hxfiltername);
00070         HyFilter.GetData(hyfiltername);
00071         
00072         gplib::cmat ExTimeFrequency, EyTimeFrequency,HxTimeFrequency,HyTimeFrequency;
00073         
00074         
00075         ExTimeFrequency = TimeFrequency(MtuData.GetData().GetEx().GetData().begin(),MtuData.GetData().GetEx().GetData().end(),seglength,Hamming());
00076         EyTimeFrequency = TimeFrequency(MtuData.GetData().GetEy().GetData().begin(),MtuData.GetData().GetEy().GetData().end(),seglength,Hamming());
00077         HxTimeFrequency = TimeFrequency(MtuData.GetData().GetHx().GetData().begin(),MtuData.GetData().GetHx().GetData().end(),seglength,Hamming());
00078         HyTimeFrequency = TimeFrequency(MtuData.GetData().GetHy().GetData().begin(),MtuData.GetData().GetHy().GetData().end(),seglength,Hamming());
00079         
00080         const unsigned int nsegs = ExTimeFrequency.size1();
00081         
00082         gplib::cmat ExHxCorr(nsegs,nfreqs);
00083         gplib::cmat ExHyCorr(nsegs,nfreqs);
00084         gplib::cmat EyHxCorr(nsegs,nfreqs);
00085         gplib::cmat EyHyCorr(nsegs,nfreqs);
00086         gplib::cmat HxHxCorr(nsegs,nfreqs);
00087         gplib::cmat HxHyCorr(nsegs,nfreqs);
00088         gplib::cmat HyHyCorr(nsegs,nfreqs);
00089         gplib::cmat Hdet(nsegs,nfreqs);
00090         gplib::cmat Zxx(nsegs,nfreqs);
00091         gplib::cmat Zxy(nsegs,nfreqs);
00092         gplib::cmat Zyx(nsegs,nfreqs);
00093         gplib::cmat Zyy(nsegs,nfreqs);
00094         gplib::rmat Weights(nsegs,nfreqs);
00095         ofstream logfile((datafilename+".log").c_str());        
00096         MeanOutput.AssignAll(nfreqs);
00097         MedianOutput.AssignAll(nfreqs);
00098         
00099         NcFile combrescdf((datafilename+".zxyarg.nc").c_str(),NcFile::Replace);
00100         NcDim* xd = combrescdf.add_dim("x",nfreqs);
00101         NcDim* yd = combrescdf.add_dim("y",nsegs);
00102         NcVar* x = combrescdf.add_var("x",ncFloat,xd);
00103         NcVar* y = combrescdf.add_var("y",ncFloat,yd);
00104         NcVar* z = combrescdf.add_var("z",ncFloat,xd,yd);
00105         float *xvals = new float[xd->size()];
00106         float *yvals = new float[yd->size()];
00107         float *zvals = new float[xd->size()*yd->size()];
00108         
00109         for (size_t i = 0; i < nfreqs; ++i)
00110         {
00111                         ofstream zxyfile((datafilename+"_zxy_"+stringify(i*freqstep)).c_str());
00112                         ofstream zxyabsfile((datafilename+"_zxyabs_"+stringify(i*freqstep)).c_str());
00113                         ofstream zxyargfile((datafilename+"_zxyarg_"+stringify(i*freqstep)).c_str());
00114                         ofstream hdetfile((datafilename+"_hdet_"+stringify(i*freqstep)).c_str());
00115                         for (size_t j = 0; j < nsegs; ++j)
00116                         {
00117                                 ExTimeFrequency(j,i) /= ExFilter.GetFilterCoeff().at(i);
00118                                 EyTimeFrequency(j,i) /= EyFilter.GetFilterCoeff().at(i);
00119                                 HxTimeFrequency(j,i) /= HxFilter.GetFilterCoeff().at(i);
00120                                 HyTimeFrequency(j,i) /= HyFilter.GetFilterCoeff().at(i);
00121                                 ExHxCorr(j,i) =(ExTimeFrequency(j,i) * conj(HxTimeFrequency(j,i)));
00122                                 ExHyCorr(j,i) =(ExTimeFrequency(j,i) * conj(HyTimeFrequency(j,i)));
00123                                 EyHxCorr(j,i) =(EyTimeFrequency(j,i) * conj(HxTimeFrequency(j,i)));
00124                                 EyHyCorr(j,i) =(EyTimeFrequency(j,i) * conj(HyTimeFrequency(j,i)));
00125                                 HxHyCorr(j,i) =(HxTimeFrequency(j,i) * conj(HyTimeFrequency(j,i)));
00126                                 HxHxCorr(j,i) =(HxTimeFrequency(j,i) * conj(HxTimeFrequency(j,i)));
00127                                 HyHyCorr(j,i) =(HyTimeFrequency(j,i) * conj(HyTimeFrequency(j,i)));
00128                         }
00129                         for (size_t j = 0; j < nsegs-navgwindows; ++j)
00130                         {
00131                                 complex<double> hxhxavg = accumulate(column(HxHxCorr,i).begin()+j,column(HxHxCorr,i).begin()+j+navgwindows,complex<double>());
00132                                 complex<double> exhxavg = accumulate(column(ExHxCorr,i).begin()+j,column(ExHxCorr,i).begin()+j+navgwindows,complex<double>());
00133                                 complex<double> exhyavg = accumulate(column(ExHyCorr,i).begin()+j,column(ExHyCorr,i).begin()+j+navgwindows,complex<double>());
00134                                 complex<double> eyhxavg = accumulate(column(EyHxCorr,i).begin()+j,column(EyHxCorr,i).begin()+j+navgwindows,complex<double>());
00135                                 complex<double> eyhyavg = accumulate(column(EyHyCorr,i).begin()+j,column(EyHyCorr,i).begin()+j+navgwindows,complex<double>());
00136                                 complex<double> hxhyavg = accumulate(column(HxHyCorr,i).begin()+j,column(HxHyCorr,i).begin()+j+navgwindows,complex<double>());
00137                                 complex<double> hyhyavg = accumulate(column(HyHyCorr,i).begin()+j,column(HyHyCorr,i).begin()+j+navgwindows,complex<double>());
00138                                 Hdet(j,i)=(hxhxavg * hyhyavg - hxhyavg * conj(hxhyavg));
00139                                 hdetfile << setfill(' ') << setw(15) << setprecision(9) << abs(Hdet(j,i).real()) << endl;
00140                                 
00141                                 if (abs(Hdet(j,i)) > 0.0)
00142                                 {
00143                                         
00144                                         Zxx(j,i) =((exhxavg*hyhyavg - exhyavg*conj(hxhyavg))/Hdet(j,i));
00145                                         Zxy(j,i) =((exhyavg*hxhxavg - exhxavg*hxhyavg)/Hdet(j,i));
00146                                         Zyx(j,i) =((eyhxavg*hyhyavg - eyhyavg*conj(hxhyavg))/Hdet(j,i));
00147                                         Zyy(j,i) =((eyhyavg*hxhxavg - eyhxavg*hxhyavg)/Hdet(j,i));
00148                                         zvals[i*nsegs +j] =arg(Zxy(j,i));
00149                                         Weights(j,i) = 1.0;
00150                                         zxyfile << setfill(' ') << setw(15) << setprecision(9) << Zxy(j,i).real() << " ";
00151                                         zxyfile << setfill(' ') << setw(15) << setprecision(9) << Zxy(j,i).imag() << endl; 
00152                                 }
00153                                 else
00154                                 {
00155                                         zvals[i*nsegs +j] = -100;
00156                                         Zxx(j,i) = 0.0;
00157                                         Zxy(j,i) = 0.0;
00158                                         Zyx(j,i) = 0.0;
00159                                         Zyy(j,i) = 0.0;
00160                                         Weights(j,i) = 0.0;
00161                                 }
00162                                 zxyabsfile << setfill(' ') << setw(15) << setprecision(9) << abs(Zxy(j,i)) << endl;
00163                                 zxyargfile << setfill(' ') << setw(15) << setprecision(9) << arg(Zxy(j,i)) << endl;
00164                         }
00165                         
00166                         ofstream weightfile((datafilename+".weights").c_str());
00167                         double currfrequency = i*freqstep;
00168                         xvals[i] = currfrequency;
00169                         complex<double> currzxx,currzxy,currzyx,currzyy;
00170                         const int goodcount = sum(column(Weights,i));
00171                         double normalize = 1./goodcount;
00172                         
00173                         cout << "Freq: " << currfrequency << " Good: " << goodcount << " Total: " << Weights.size1() << endl;
00174                         weightfile << sum(column(Weights,i)) << endl;
00175                         currzxx = prec_inner_prod(column(Zxx,i),column(Weights,i));
00176                         currzxx *= normalize;
00177                         currzxy = prec_inner_prod(column(Zxy,i),column(Weights,i));
00178                         currzxy *= normalize;
00179                         currzyx = prec_inner_prod(column(Zyx,i),column(Weights,i));
00180                         currzyx *= normalize;
00181                         currzyy = prec_inner_prod(column(Zyy,i),column(Weights,i));
00182                         currzyy *= normalize;
00183                         MTTensor CurrZ(currzxx,currzxy,currzyx,currzyy,currfrequency);
00184                         MeanOutput.SetMTData().at(i) = CurrZ;
00185                         zxyfile << endl << endl; // write mean to file as new dataset
00186                         zxyfile << setfill(' ') << setw(15) << setprecision(9) << currzxy.real() << " ";
00187                         zxyfile << setfill(' ') << setw(15) << setprecision(9) << currzxy.imag() << endl;
00188                         zxyfile << endl << endl;  //write median to file as new dataset
00189                         vector<double> zxyreal(goodcount,0.0);
00190                         vector<double> zxyimag(goodcount,0.0);
00191                         int currindex = 0;
00192                         for (int j = 0; j < nsegs; ++j)
00193                         {
00194                                 yvals[j] = j;
00195                                 if (Weights(j,i) > 0.0)
00196                                 {
00197                                         zxyreal.at(currindex) = Zxy(j,i).real();
00198                                         zxyimag.at(currindex) = Zxy(j,i).imag();
00199                                         ++currindex;
00200                                 }
00201                         }
00202                         long double zxyrmedian = Median(zxyreal.begin(),zxyreal.end());
00203                         long double zxyimedian = Median(zxyimag.begin(),zxyimag.end());
00204                         zxyfile << setfill(' ') << setw(15) << setprecision(9) << zxyrmedian << " ";
00205                         zxyfile << setfill(' ') << setw(15) << setprecision(9) << zxyimedian << endl;
00206                         currzxy.real() = zxyrmedian;
00207                         currzxy.imag() = zxyimedian;
00208                         MTTensor NewZ(currzxx,currzxy,currzyx,currzyy,currfrequency);
00209                         MedianOutput.SetMTData().at(i) = NewZ;
00210         }       
00211         x->put(xvals,xd->size());
00212         y->put(yvals,yd->size());
00213         z->put(zvals,z->edges());
00214         MeanOutput.WriteAsMtt((datafilename+".mean").c_str());
00215         MedianOutput.WriteAsMtt((datafilename+"median").c_str());
00216 }