miscfunc.cpp

Go to the documentation of this file.

#include "miscfunc.h"
#include "CFatalException.h"
#include <numeric>
#include <algorithm>
#include <fftw3.h>
#include <gsl/gsl_math.h>
#include <fstream>
#include <iomanip>
#include <boost/bind.hpp>
#include <complex>
using namespace std;

void Hilbert(std::vector<double>::iterator Inputbegin,std::vector<double>::iterator Inputend, std::vector<double> &Output)
{
        fftw_complex *timedomain;
        fftw_complex *freqdomain;
        fftw_plan p;
        double swapstore;
        const unsigned int size = distance(Inputbegin,Inputend);
        if (Output.size() != size)
                Output.assign(size,0);
                
        timedomain = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * size);
        freqdomain = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * size);
        p = fftw_plan_dft_1d(size,timedomain,freqdomain,-1,FFTW_ESTIMATE);      
        for (unsigned int i=0; i < size; ++i)
        {
                timedomain[i][0] = *(Inputbegin+i) * ( 0.5 - 0.5 * cos(2 * PI *i /size));
                timedomain[i][1] = 0;
        }
        fftw_execute(p);
        for (unsigned int i=0; i < size; ++i)
        {
                swapstore = freqdomain[i][0];
                if (i < size/2)
                {
                        freqdomain[i][0] = -freqdomain[i][1];
                        freqdomain[i][1] = swapstore;
                }
                else
                {
                        freqdomain[i][0] = freqdomain[i][1];
                        freqdomain[i][1] = -swapstore;
                }
        }
        p = fftw_plan_dft_1d(size,freqdomain,timedomain,1,FFTW_ESTIMATE);
        fftw_execute(p);        
        for (unsigned int i = 0; i < size; ++i)
                Output.at(i) = 1./double(size) * timedomain[i][0];
        fftw_destroy_plan(p);
        fftw_free(timedomain);
        fftw_free(freqdomain);
}



void NormEnvelope(std::vector<double>::iterator Inputbegin,std::vector<double>::iterator Inputend , std::vector<double> &Output)
{
        const int size = distance(Inputbegin,Inputend);
        vector<double> Transformed(size);
        double maximum;
        Output.assign(size,0);
        Hilbert(Inputbegin,Inputend,Transformed);
        for (int i = 0; i < size; ++i)
                Output.at(i) = sqrt(pow(*(Inputbegin+i),2)+pow(Transformed.at(i),2));
        maximum = *max_element(Output.begin(),Output.end());
        for (int i = 0; i < size; ++i)
                Output.at(i) /= maximum;
}


double Cross(ttsdata Master, ttsdata Corr, const int startpoint, const int endpoint)
{
        const int size = endpoint - startpoint;
        double cross = 0;
        double auto1 = 0;
        double auto2 = 0;
        const double mean1 = accumulate(Master.begin()+startpoint,Master.begin()+endpoint,0.)/size;
        const double mean2 = accumulate(Corr.begin()+startpoint,Corr.begin()+endpoint,0.)/size;
        for (int i = startpoint; i < endpoint; ++i)
        {
                cross += (Master.at(i)-mean1)*(Corr.at(i)-mean2);
                auto1 += gsl_pow_2(Master.at(i)-mean1);
                auto2 += gsl_pow_2(Corr.at(i)-mean2);
        }
        if (gsl_fcmp(0.0,auto1,GSL_SQRT_DBL_EPSILON) || gsl_fcmp(0.0,auto2,GSL_SQRT_DBL_EPSILON) )
                return(cross / (sqrt (auto1) * sqrt (auto2)));
        else 
                return 0.0;
}

void index(const trealdata data, vector<int> &result)
{
    int i,j;
    int a;
        const int freqnum = result.size();
    for (i= 0; i < freqnum; ++i)
                result.at(i) = i;
    for (j = 1; j < freqnum; ++j)
    {
                a = result.at(j);
                i = j - 1;
                while (i >= 0 && data.at(result.at(i)) > data.at(a))
                {
                result.at(i+1) = result.at(i);
                --i;
                }
                result.at(i+1) = a;
    }
}

---