fftw_test.cpp

Go to the documentation of this file.

#include "CTsSpectrum.h"
#include "types.h"
#include "CMtuFormat.h"
#include <fstream>
#include <iostream>
#include <iomanip>
using namespace std;

int main ()
{
        CMtuFormat MtuData;
        ttsdata TimeSeries;
        ttsdata ExHyCross;
        tcompdata FreqDomain;
        tcompdata ExStackedSpec;
        tcompdata HyStackedSpec;
        tcompdata ExHyCorrel;
        ofstream outspec;
        CTsSpectrum ExSpectrum;
        CTsSpectrum HySpectrum;
        CTsSpectrum CorrelSpectrum;
        string infilename;
        const int seglength = 4096;
        int nsegs;
        double mean;
        
        cout << "Infile Name: ";
        cin >> infilename;
        
        MtuData.GetData(infilename); 
        
        ExSpectrum.TsData = &MtuData.Ex.data;
        ExSpectrum.FreqData = &FreqDomain;
        ExSpectrum.MultiCalc = false;
        ExStackedSpec.assign(seglength/2,0);
        nsegs = ExSpectrum.TsData->size()/seglength-1;
        cout << "Nsegs: " << nsegs << endl;
        for (int j = 0; j < nsegs; ++j)
        {
                mean = 0;
                for (int k = j*seglength; k < (j+1)*seglength; ++k)
                        mean += ExSpectrum.TsData->at(k);
                for (int k = j*seglength; k < (j+1)*seglength; ++k)
                {
                        ExSpectrum.TsData->at(k) -= mean;
                        ExSpectrum.TsData->at(k) *= 0.54 - 0.46 * cos(2 * PI *j /seglength);
                }
                ExSpectrum.CalcSpectrum(j*seglength,(j+1)*seglength);
                for (int k = 0; k < seglength/2; ++k)
                        ExStackedSpec.at(k) += FreqDomain[k]; 
        }
        
        HySpectrum.TsData = &MtuData.Hy.data;
        HySpectrum.FreqData = &FreqDomain;
        HySpectrum.MultiCalc = false;
        HyStackedSpec.assign(seglength/2,0);
        nsegs = HySpectrum.TsData->size()/seglength-1;
        for (int j = 0; j < nsegs; ++j)
        {
                mean = 0;
                for (int k = j*seglength; k < (j+1)*seglength; ++k)
                        mean += HySpectrum.TsData->at(k);
                for (int k = j*seglength; k < (j+1)*seglength; ++k)
                        HySpectrum.TsData->at(k) -= mean;
        
                HySpectrum.CalcSpectrum(j*seglength,(j+1)*seglength);
                for (int k = 0; k < seglength/2; ++k)
                        HyStackedSpec.at(k) += FreqDomain[k]; 
        }
        cout << "Spectra calculated" << endl << flush;
        for (int i =0; i < ExStackedSpec.size(); ++i)
                ExHyCorrel.push_back(ExStackedSpec.at(i) * conj(HyStackedSpec.at(i)));
        CorrelSpectrum.TsData = &ExHyCross;
        CorrelSpectrum.FreqData = &ExHyCorrel;
        CorrelSpectrum.MultiCalc = false;
        CorrelSpectrum.CalcTimeSeries();
        outspec.open((infilename+".spec.ex").c_str());
        for (int i = 0; i < ExStackedSpec.size(); ++i)
        { 
                outspec << double(i)/(seglength * MtuData.samplerate) << " " << norm(ExStackedSpec.at(i)) << endl;
        }
        outspec.close();
        outspec.open((infilename+".spec.hy").c_str());
        for (int i = 0; i < HyStackedSpec.size(); ++i)
        { 
                outspec << double(i)/(seglength * MtuData.samplerate) << " " << norm(HyStackedSpec.at(i)) << endl;
        }
        outspec.close();
        outspec.open((infilename+".exhy.corr").c_str());
        for (int i = 0; i < ExHyCorrel.size(); ++i)
        { 
                //outspec << double(i)/(seglength * MtuData.samplerate) << " " << abs(CorrelSpectrum.FreqData->at(i))
                //      /sqrt(abs(HySpectrum.FreqData->at(0)) * abs(ExSpectrum.FreqData->at(0))) << endl;
                outspec << i << " " << setprecision(15) << CorrelSpectrum.TsData->at(i) << endl;
        }
        outspec.close();
}

---