CDespike.cpp

Go to the documentation of this file.

#include "CDespike.h"
#include <iostream>
#include <iomanip>
#include <string>
#include <fstream>
#include "tsfuncs.h"

using namespace std;
CDespike::CDespike()
{
        ReachedEnd = false;
        CurrentIndex = 0;
}
CDespike::~CDespike()
{
        Data = NULL;
        Time = NULL;
}

void CDespike::Restart()
{
        ReachedEnd = false;
        CurrentTime = Time->begin();
        CurrentSpike = Data->begin();   
        CurrentIndex = 0;
}

void CDespike::RemoveSpikes(const titdata begin, const titdata end)
{
        CurrentSpike = begin;
        CurrentIndex = begin - Data->begin();
        CurrentTime = Time->begin() + (begin - Data->begin());
        ReachedEnd = false;
        
        vector<double>::iterator datait;
        vector<double>::iterator stackit;
        
        //copy(Data->begin(),Data->end(),back_insert_iterator<vector<double> > (CleanData));
        //Data = &CleanData;
        while (ReachedEnd == false)
        {
                FindNextSpike(begin,end);
                if ( (distance(begin,CurrentSpike) > TrailPoints) && 
                         (distance(CurrentSpike,end) > (SpikeLength - TrailPoints)) )
                {
                        datait = CurrentSpike;
                        advance(datait, -TrailPoints);
                        stackit = SpikeForm.begin();
                        for (int j = 0; j < SpikeLength; ++j)
                        {
                                *(datait) -= *(stackit);
                                //*(datait) += mean;
                                datait++;
                                stackit++;
                        }
                }
        }
}

/*void CDespike::RemoveSpikes()
{
        Restart();
        titdata current;
        titdata next;
        titdata last;
        titdata currenttime;
        titdata lasttime;
        double max = 0;
        int nspikes = 0;
        double prediff,postdiff,timediff;
        titspikeindex index;
        
        current = Data->begin()+1;
        currenttime = Time->begin()+1;
        next = current+1;
        last = current-1;
        lasttime = currenttime -1;
        /*while ( next != Data->end())
        {
                postdiff = *next - *current;
                prediff =  *current - *last;
                timediff = *currenttime - *lasttime;
                if (((fabs(prediff) > HeightThreshold) || (fabs(postdiff) > HeightThreshold))
                        && (timediff <= 2*samplerate) &&( (postdiff * prediff) < 0))
                {       
                        cout << " Before " << *current << " " << *last << " " << *next << " " << postdiff <<" " << prediff;
                        *current = (*next + *last)/2.;
                        cout << " After: " << *current << " " << *last << " " << *next << endl;
                        nspikes++;
                        if (fabs(prediff) > max)
                                max = prediff;
                }*/
/*              index = SpikeIndices.begin();
                while (distance(index,SpikeIndices.end()) > 2)
                {
                        postdiff = Data->at(*(index) +1) - Data->at(*index);
                        prediff = Data->at(*index) - Data->at(*(index) -1);
                        /*Data->at(*index) = (Data->at(*(index) +1) + Data->at(*(index) +2)
                                 + Data->at(*(index) -1) + Data->at(*(index) -2))/4.;*/
/*                      Data->at(*index) = (Data->at(*(index) +1) 
                                 + Data->at(*(index) -1))/2.;
                        index++;
                        nspikes++;
                        if (fabs(prediff) > max)
                                max = prediff;
                }
                /*current++;
                next++;
                last++;
                currenttime++;
                lasttime++;
        }*/
/*      cout << "Max: " << max <<endl;
        cout << "NSpikes: " << nspikes <<endl;
        cout << "NPoints: " << Data->size() <<endl;
}*/
int CDespike::StackSpikes()
{
        return StackSpikes(Data->begin(),Data->end());
}

int CDespike::WriteSpikeForm(string filename)
{
        ofstream outfile;
        
        outfile.open((filename+".spikeform").c_str());
        for (int i = 0; i < SpikeForm.size(); ++i)
        {
                outfile.precision(5);
                outfile << setw(15) << setfill(' ') << i;
                outfile << setw(15) << setfill(' ') << SpikeForm.at(i);
        }
        outfile.close();
        return 0;
}

int CDespike::StackSpikes(const titdata begin, const titdata end)
{
        int NoEstimates = 0;
        titdata StackStart;
        titdata AvgStart;
        double SpikeMean = 0;
        
        CurrentSpike = begin;
        CurrentIndex = distance(Data->begin(),begin);
        CurrentTime = Time->begin() + distance(Data->begin(),begin);
        ReachedEnd = false;
        SpikeIndices.clear();
        SpikeForm.assign(SpikeLength,0);
        SpikeMax = 0;
        
        
        while (ReachedEnd == false)
        {
                FindNextSpike(begin,end);
                SpikeIndices.push_back(CurrentIndex);
                if ( (distance(begin,CurrentSpike) > TrailPoints) && 
                         (distance(CurrentSpike,end) > (SpikeLength - TrailPoints)) )
                {
                        AvgStart = CurrentSpike;
                        advance(AvgStart, -TrailPoints);
                        StackStart = SpikeForm.begin();
                        for (int j = 0; j < SpikeLength; ++j)
                        {
                                *(StackStart) += *(AvgStart);
                                StackStart++;
                                AvgStart++;
                        }
                        NoEstimates++;
                        CurrentSpike += SpikeLength;
                        CurrentIndex += SpikeLength;
                        CurrentTime += SpikeLength;
                        if (CurrentSpike > end)
                                ReachedEnd = true;
                }
        }
        
        for (StackStart = SpikeForm.begin(); StackStart != SpikeForm.end(); ++StackStart)
        {
                *(StackStart) /= NoEstimates;
                SpikeMean += *(StackStart);
                if (*StackStart > SpikeMax)
                        SpikeMax = *StackStart;
        }
        SpikeMean /= SpikeLength;
        for (StackStart = SpikeForm.begin(); StackStart != SpikeForm.end(); ++StackStart)
        {
                *(StackStart) -= SpikeMean;
        }
        cout << "Stacked Estimates :" << NoEstimates << endl;
        return NoEstimates;     
}

void CDespike::FindAllSpikes()
{
        SpikeIndices.clear();
        Restart();
        while (ReachedEnd == false)
        {
                FindNextSpike();
                SpikeIndices.push_back(CurrentIndex);
        }
                        
}

void CDespike::FindNextSpike()
{
        FindNextSpike(Data->begin(),Data->end());
}

void CDespike::FindNextSpike(const titdata begin, const titdata end)
{
        titdata next;
        titdata last;
        double prediff, postdiff;
        
        next = CurrentSpike;
        last = CurrentSpike;
        next++;
        last--;
        
        do
                {
                        CurrentIndex++;
                        CurrentSpike++;
                        CurrentTime++;
                        next++;
                        last++;
                        prediff = *CurrentSpike - *last;
                        postdiff = *next -  *CurrentSpike;
                        
                }
        while ( (next < end)&& !(IsSpikePreRel(prediff,postdiff,*CurrentSpike,HeightThreshold)));
        if (next >= end)
                ReachedEnd = true;               
}

---