fftw_test.cpp

Go to the documentation of this file.
00001 #include "CTsSpectrum.h"
00002 #include "types.h"
00003 #include "CMtuFormat.h"
00004 #include <fstream>
00005 #include <iostream>
00006 #include <iomanip>
00007 using namespace std;
00008 
00009 int main ()
00010 {
00011         CMtuFormat MtuData;
00012         ttsdata TimeSeries;
00013         ttsdata ExHyCross;
00014         tcompdata FreqDomain;
00015         tcompdata ExStackedSpec;
00016         tcompdata HyStackedSpec;
00017         tcompdata ExHyCorrel;
00018         ofstream outspec;
00019         CTsSpectrum ExSpectrum;
00020         CTsSpectrum HySpectrum;
00021         CTsSpectrum CorrelSpectrum;
00022         string infilename;
00023         const int seglength = 4096;
00024         int nsegs;
00025         double mean;
00026         
00027         cout << "Infile Name: ";
00028         cin >> infilename;
00029         
00030         MtuData.GetData(infilename); 
00031         
00032         ExSpectrum.TsData = &MtuData.Ex.data;
00033         ExSpectrum.FreqData = &FreqDomain;
00034         ExSpectrum.MultiCalc = false;
00035         ExStackedSpec.assign(seglength/2,0);
00036         nsegs = ExSpectrum.TsData->size()/seglength-1;
00037         cout << "Nsegs: " << nsegs << endl;
00038         for (int j = 0; j < nsegs; ++j)
00039         {
00040                 mean = 0;
00041                 for (int k = j*seglength; k < (j+1)*seglength; ++k)
00042                         mean += ExSpectrum.TsData->at(k);
00043                 for (int k = j*seglength; k < (j+1)*seglength; ++k)
00044                 {
00045                         ExSpectrum.TsData->at(k) -= mean;
00046                         ExSpectrum.TsData->at(k) *= 0.54 - 0.46 * cos(2 * PI *j /seglength);
00047                 }
00048                 ExSpectrum.CalcSpectrum(j*seglength,(j+1)*seglength);
00049                 for (int k = 0; k < seglength/2; ++k)
00050                         ExStackedSpec.at(k) += FreqDomain[k]; 
00051         }
00052         
00053         HySpectrum.TsData = &MtuData.Hy.data;
00054         HySpectrum.FreqData = &FreqDomain;
00055         HySpectrum.MultiCalc = false;
00056         HyStackedSpec.assign(seglength/2,0);
00057         nsegs = HySpectrum.TsData->size()/seglength-1;
00058         for (int j = 0; j < nsegs; ++j)
00059         {
00060                 mean = 0;
00061                 for (int k = j*seglength; k < (j+1)*seglength; ++k)
00062                         mean += HySpectrum.TsData->at(k);
00063                 for (int k = j*seglength; k < (j+1)*seglength; ++k)
00064                         HySpectrum.TsData->at(k) -= mean;
00065         
00066                 HySpectrum.CalcSpectrum(j*seglength,(j+1)*seglength);
00067                 for (int k = 0; k < seglength/2; ++k)
00068                         HyStackedSpec.at(k) += FreqDomain[k]; 
00069         }
00070         cout << "Spectra calculated" << endl << flush;
00071         for (int i =0; i < ExStackedSpec.size(); ++i)
00072                 ExHyCorrel.push_back(ExStackedSpec.at(i) * conj(HyStackedSpec.at(i)));
00073         CorrelSpectrum.TsData = &ExHyCross;
00074         CorrelSpectrum.FreqData = &ExHyCorrel;
00075         CorrelSpectrum.MultiCalc = false;
00076         CorrelSpectrum.CalcTimeSeries();
00077         outspec.open((infilename+".spec.ex").c_str());
00078         for (int i = 0; i < ExStackedSpec.size(); ++i)
00079         { 
00080                 outspec << double(i)/(seglength * MtuData.samplerate) << " " << norm(ExStackedSpec.at(i)) << endl;
00081         }
00082         outspec.close();
00083         outspec.open((infilename+".spec.hy").c_str());
00084         for (int i = 0; i < HyStackedSpec.size(); ++i)
00085         { 
00086                 outspec << double(i)/(seglength * MtuData.samplerate) << " " << norm(HyStackedSpec.at(i)) << endl;
00087         }
00088         outspec.close();
00089         outspec.open((infilename+".exhy.corr").c_str());
00090         for (int i = 0; i < ExHyCorrel.size(); ++i)
00091         { 
00092                 //outspec << double(i)/(seglength * MtuData.samplerate) << " " << abs(CorrelSpectrum.FreqData->at(i))
00093                 //      /sqrt(abs(HySpectrum.FreqData->at(0)) * abs(ExSpectrum.FreqData->at(0))) << endl;
00094                 outspec << i << " " << setprecision(15) << CorrelSpectrum.TsData->at(i) << endl;
00095         }
00096         outspec.close();
00097 }