MtuDespike.cpp

Go to the documentation of this file.

#include "MtuFormat.h"
#include "CDespike.h"
#include "CSpikeStats.h"
#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>
#include <sstream>
#include <algorithm>

const double filterfraction=0.05;
const int sectionlength = 4800;
const int spikelength = 50;
const int trailpoints = 10;
const int nValueBins = 1000;
const int nDerivBins = 1000;
const int MinEstimates = 2;
const int iterations = 1;


using namespace std;


void CleanComponent(TimeSeriesComponent &CurrentComponent, CDespike &Despiker,
        CSpikeStats &SpikeStats, ttsdata &NoiseTs)
{
        titdata sectionstart;
        int estimates = 0;
        int insertpoints = 0;
        titdata spikestart;
        
        SpikeStats.Data = &(CurrentComponent.GetData());
        SpikeStats.AnalyseSpikes();
        Despiker.HeightThreshold = SpikeStats.HeightThreshold;
        Despiker.Data = &(CurrentComponent.GetData());
        for (sectionstart = CurrentComponent.GetData().begin(); 
                        (sectionstart + sectionlength) < CurrentComponent.GetData().end(); sectionstart+=sectionlength)
        {
                estimates = Despiker.StackSpikes(sectionstart,sectionstart+sectionlength);
                if (estimates > MinEstimates)
                {
                        Despiker.RemoveSpikes(sectionstart,sectionstart+sectionlength);
                        for (int i = 0; i < estimates; ++i)
                        {
                                spikestart = NoiseTs.begin() + Despiker.SpikeIndices.at(i) -trailpoints;
                                copy(Despiker.SpikeForm.begin(),Despiker.SpikeForm.end(),spikestart);
                        }
                }
        }
        /*Despiker.Restart();
        while (Despiker.ReachedEnd == false)
        {
                Despiker.FindNextSpike();
                NoiseTs.at(Despiker.CurrentIndex) = *(Despiker.CurrentSpike);
        }*/
        SpikeStats.Data = NULL;
        Despiker.Data = NULL;
}

int main(void)
{
        MtuFormat TestData;
        CDespike Despiker;
        CSpikeStats SpikeStats;
        ttsdata ExNoiseTs;
        ttsdata EyNoiseTs;
        ofstream noisetsfile;
        string infilename, outfilename,iteration;
        
        
        cout << "Inputfile: ";
        cin >> infilename;
        TestData.GetData(infilename);
        
        cout << "Outputfile: ";
        cin >> outfilename;
        
        Despiker.Time = &(TestData.GetTime());
        Despiker.samplerate = TestData.GetSamplerate();
        Despiker.SpikeLength = spikelength;
        Despiker.TrailPoints = trailpoints;
        SpikeStats.Time = &(TestData.GetTime());
        SpikeStats.nValueBins = nValueBins;
        SpikeStats.nDerivBins = nDerivBins;
        SpikeStats.samplerate = TestData.GetSamplerate();
        SpikeStats.HeightFraction = filterfraction;
        ExNoiseTs.assign(TestData.GetEx().GetData().size(),0);
        EyNoiseTs.assign(TestData.GetEy().GetData().size(),0);
        CleanComponent(TestData.GetEx(),Despiker,SpikeStats,ExNoiseTs);
        CleanComponent(TestData.GetEy(),Despiker,SpikeStats,EyNoiseTs);
        SpikeStats.WriteStats(infilename + TestData.GetEx().GetName()); 
        Despiker.WriteSpikeForm(infilename + TestData.GetEx().GetName());
        
        cout << "Writing data ..."  << flush;
        noisetsfile.open((infilename+".noise").c_str());
        for (int i = 0; i < ExNoiseTs.size(); ++i)
        {
                noisetsfile << setw(10) << ExNoiseTs.at(i) << " " << EyNoiseTs.at(i) <<  endl;
        }
        noisetsfile.close();
        TestData.WriteData(outfilename);
        SpikeStats.Time = NULL;
        Despiker.Time = NULL;
        
        cout << " done"  << endl;
}