SpikeSpectrum.cpp

Go to the documentation of this file.

#include "CMtuFormat.h"
#include "CSpikeStats.h"
#include "CTsSpectrum.h"
#include "CDespike.h"
#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include <fftw3.h>
using namespace std;
        
        const int iterations =20;
        //const int blocklength = 2400; 
        CMtuFormat TsData;
        ofstream diag;
        double samplerate;
        const int length = 100;
        const int trailratio = 4;
        
void Filter(const tcompdata &Filter, ttsdata &TimeSeries)
{
        double *timedomain;
        fftw_complex *freqdomain;
        fftw_plan pforward, pinverse;
        int seglength;
        int nsegs;
        double real, imag, mean, coeff;
        
        seglength = Filter.size() * 2;
        nsegs = TimeSeries.size()/seglength;
        timedomain = (double *)fftw_malloc(sizeof(double) * seglength);
        freqdomain = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (seglength/2+1));
        pforward = fftw_plan_dft_r2c_1d(seglength,timedomain,freqdomain,FFTW_MEASURE);
        pinverse = fftw_plan_dft_c2r_1d(seglength,freqdomain,timedomain,FFTW_MEASURE);
        for (int i = 0; i < nsegs-1; ++i)
        {
                mean = 0;
                for (int j = 0; j < seglength; ++j)
                {
                        timedomain[j] = TimeSeries.at(i*seglength +j);
                        mean += timedomain[j];
                }
                for (int j = 0; j < seglength; ++j)
                {
                        timedomain[j] -=mean/seglength;
                        timedomain[j] *= 0.54 - 0.46 * cos(2 * PI *j /seglength);
                }
                fftw_execute(pforward);
                for (int j = 0; j < seglength/2; j++)
                {
                        real = freqdomain[j][0];
                        imag = freqdomain[j][1];
                        coeff = (pow(real,2) + pow(imag,2) - norm(Filter.at(j)))/(pow(real,2) + pow(imag,2));
                        /*coeff = (pow(real - Filter.at(j).real(),2) + pow(imag - Filter.at(j).imag(),2))/
                                (pow(real - Filter.at(j).real(),2) + pow(imag - Filter.at(j).imag(),2)
                                + norm(Filter.at(j)));*/
                                if (coeff < 0)
                                        coeff = 0.0001;
                        freqdomain[j][0] = real * coeff;
                        freqdomain[j][1] = imag * coeff;        
                        if (i == 1)
                                diag << j << " " << coeff << " " << imag << " " << real << " " << Filter.at(j) << endl;
                }
                fftw_execute(pinverse);
                for (int j = 0; j < seglength; ++j)
                {
                        timedomain[j] /= seglength;
                        timedomain[j] /= 0.54 - 0.46 * cos(2 * PI *j /seglength);
                        timedomain[j] += mean/seglength;
                        TimeSeries.at(i*seglength +j)= timedomain[j];
                }       
        }
}

void WriteSpec(string name, const CTsSpectrum &Spectrum)
{
        ofstream spikefreq;
        spikefreq.open(name.c_str());
        int i = 0;
        for (tcompdata::iterator it = Spectrum.FreqData->begin(); it != Spectrum.FreqData->end(); ++it)
        {
                spikefreq.precision(10);
                spikefreq << ++i << " " << (*it).real() << " " << (*it).imag() << endl;
        }
        spikefreq.close();      
}

void WritePower(string name, const CTsSpectrum &Spectrum)
{
        ofstream spikefreq;
        spikefreq.open(name.c_str());
        int i = -1;
        double freqstep = 1./(samplerate * length);
        for (tcompdata::iterator it = Spectrum.FreqData->begin(); it != Spectrum.FreqData->end(); ++it)
        {
                spikefreq.precision(10);
                spikefreq << (++i)*freqstep << " " << norm(*it) << endl;
        }
        spikefreq.close();      
}

void WriteSpike(ttsdata &SpikeForm, string name)
{
        ofstream spike;
        spike.open(name.c_str());
        int i = 0;
        for (ttsdata::iterator it = SpikeForm.begin(); it != SpikeForm.end(); ++it)
        {
                spike.precision(10);
                spike << ++i << " " << *it << endl;
        }
        spike.close();  
}

void ApplyFilter(ttsdata &Component, string name, ttsdata &Time)
{
                double mean = 0;
                CSpikeStats SpikeStats;
                CDespike Despiker;
                CTsSpectrum SpikeSpec;
        
                SpikeStats.Time = &(Time);
                SpikeStats.nValueBins = 1000;
                SpikeStats.nDerivBins = 1000;
                SpikeStats.HeightFraction = 0.01;
                SpikeStats.samplerate = samplerate;
                Despiker.Time = &(Time);
                //Despiker.SpikeLength = 10000;
                //Despiker.TrailPoints = 5000;
                Despiker.SpikeLength = length;
                Despiker.TrailPoints = length/trailratio;
                Despiker.samplerate = samplerate;
                SpikeSpec.TsData = &Despiker.SpikeForm;
                SpikeSpec.MultiCalc = false;
                SpikeSpec.FreqData = new tcompdata;
                SpikeStats.Data = &Component;
                Despiker.Data = &Component;
                SpikeStats.AnalyseSpikes();
                Despiker.HeightThreshold = SpikeStats.HeightThreshold;
                Despiker.StackSpikes();
                //WriteSpike(Despiker.SpikeForm,name+"_spike.dat");
                //for (int i = 0; i < Despiker.SpikeForm.size(); ++i)
                //      mean += Despiker.SpikeForm.at(i);
                for (int i = 0; i < Despiker.SpikeForm.size(); ++i)
                {
                        //Despiker.SpikeForm.at(i) -= mean/Despiker.SpikeForm.size();
                        Despiker.SpikeForm.at(i) *= 0.54 - 0.46 * cos(2 * PI *i 
                                /Despiker.SpikeForm.size());
                }
                SpikeSpec.CalcSpectrum();
                //for (int i = 0; i < SpikeSpec.FreqData->size()/100; ++i)
                //      SpikeSpec.FreqData->at(i) = 0;
                //WriteSpec(name+"_spikefreq.dat",SpikeSpec);   
                //WritePower(name+"_spikepow.dat",SpikeSpec);
                Filter(*SpikeSpec.FreqData,Component);
}

int main()
{
        string infilename,temp;
        ostringstream converter;
        
        double filterfraction, currentfraction;
        
        do
        {
        cout << "Inputfile: ";
        cin >> infilename;
        }
        while (TsData.GetData(infilename.c_str()));
        
        diag.open((infilename+".diag").c_str());
        samplerate = TsData.samplerate;
        for (int i = 0; i < iterations; ++i)
        {
                converter << i;
                temp = infilename + converter.str() + "ey";             
                ApplyFilter(TsData.Ey.data,temp,TsData.t);
                temp = infilename + converter.str() + "ex";             
                ApplyFilter(TsData.Ex.data,temp,TsData.t);
                //temp = infilename + converter.str() + "hx";           
                //ApplyFilter(TsData.Hx.data,temp);
                //temp = infilename + converter.str() + "hy";           
                //ApplyFilter(TsData.Hy.data,temp);
                
        }
        TsData.WriteData(("opt"+infilename).c_str());
        diag.close();
}