OptimumFilter.cpp

Go to the documentation of this file.

#include <CMtuFormat.h>
#include "types.h"
#include <iostream>
#include <fstream>
#include <algorithm>
#include <fftw3.h>

using namespace std;
void Filter(const tcompdata &Filter, ttsdata &TimeSeries)
{
        double *timedomain;
        fftw_complex *freqdomain;
        fftw_plan pforward, pinverse;
        int seglength;
        int nsegs;
        double real, imag, mean;
        
        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.5 * (1 - cos(2 * PI *j /seglength));
                }
                fftw_execute(pforward);
                for (int j = 0; j < seglength/2; j+=2)
                {
                        real = freqdomain[j][0];
                        imag = freqdomain[j][1];
                        real *= (pow(freqdomain[j][0] - Filter.at(j).real(),2) + pow(freqdomain[j][1] - Filter.at(j).imag(),2))/
                                ((pow(freqdomain[j][0] - Filter.at(j).real(),2) + pow(freqdomain[j][1] - Filter.at(j).imag(),2))
                                + norm(Filter.at(j)));
                        imag *= (pow(freqdomain[j][0] - Filter.at(j).real(),2) + pow(freqdomain[j][1] - Filter.at(j).imag(),2))/
                                ((pow(freqdomain[j][0] - Filter.at(j).real(),2) + pow(freqdomain[j][1] - Filter.at(j).imag(),2))
                                + norm(Filter.at(j)));
                        freqdomain[j][0] = real;
                        freqdomain[j][1] = imag;        
                }
                fftw_execute(pinverse);
                for (int j = 0; j < seglength; ++j)
                {
                        TimeSeries.at(i*seglength +j)= timedomain[j]/seglength;
                }       
        }
}

void ReadSpike(string name, tcompdata &Spike)
{
        ifstream spikefile;
        double real,imag,dummy;
        
        spikefile.open(name.c_str());
        while (spikefile.good())
        {
                spikefile >> dummy;
                spikefile >> real;
                spikefile >> imag;
                if (spikefile.good())
                        Spike.push_back(real + I * imag);
        } 
        spikefile.close();
}

int main()
{
        CMtuFormat TsData;
        string infilename;
        tcompdata SpikePowerex, SpikePowerey, SpikePowerhx, SpikePowerhy;
        
        
        cout << "TS File: ";
        cin >> infilename;
        
        TsData.GetData(infilename.c_str());
        
        ReadSpike(infilename+"_spikefreqex.dat", SpikePowerex);
        ReadSpike(infilename+"_spikefreqey.dat", SpikePowerey);
        ReadSpike(infilename+"_spikefreqhx.dat", SpikePowerhx);
        ReadSpike(infilename+"_spikefreqhy.dat", SpikePowerhy);
        Filter(SpikePowerex,TsData.Ex.data);
        Filter(SpikePowerey,TsData.Ey.data);
        Filter(SpikePowerhx,TsData.Hx.data);
        Filter(SpikePowerhy,TsData.Hy.data);

        TsData.WriteData(("opt"+infilename).c_str());
}