icatest.cpp

Go to the documentation of this file.

#include <iostream>
#include <itpp/base/fastica.h>
#include <itpp/base/mat.h>
#include <itpp/base/sigfun.h>
#include <string>
#include <fstream>
#include "UniformRNG.h"
#include "TimeSeriesData.h"
#include "statutils.h"
using namespace std;

void WriteToFile(const itpp::Vec<double> &orig,const itpp::Vec<double> &in,const  itpp::Vec<double> &out,const string &base)
{
        ofstream origts((base+".orig").c_str());
        ofstream origspec((base+".spec.orig").c_str());
        ofstream ints((base+".in").c_str());
        ofstream inspec((base+".spec.in").c_str());
        ofstream outts((base+".out").c_str());
        ofstream outspec((base+".spec.out").c_str());
        const int size = in.length();
        itpp::Vec<double> insp = itpp::spectrum(in,2048,0);
        itpp::Vec<double> outsp = itpp::spectrum(out,2048,0);
        itpp::Vec<double> origsp = itpp::spectrum(orig,2048,0);
        for (int i = 0; i < size; ++i)
        {
                origts << orig(i) <<  endl;
                origspec << origsp(i) <<  endl;
                ints << in(i) <<  endl;
                inspec << insp(i) <<  endl;
                outts << out(i) <<  endl;
                outspec << outsp(i) <<  endl;
        }
}


int main()
{
        const int nchannels = 4;
        const size_t tslength = 2000;
        itpp::mat Mix(nchannels,nchannels);
        UniformRNG Random;
        
        for (int i = 0; i < nchannels; ++i)
                for (int j = 0; j < nchannels; ++j)
                        Mix(i,j) = Random.GetNumber();
        
        itpp::mat Input(nchannels,tslength);
        itpp::mat Sources(nchannels,tslength);
        for (size_t i = 0; i < tslength; ++i)
        {
                /*Sources(0,i) = sin(i/13.678);
                Sources(1,i) = 2 * Random.GetNumber();
                Sources(2,i) = 300 * (1 - static_cast<double>(i)/tslength);
                Sources(3,i) = sqrt(i) * cos(i/241.646356); */
        }
    for (size_t i = 0; i < tslength; ++i)
    {
        double number = 0.5 - Random.GetNumber();
        Input(0,i) = sin(i/13.678) + 0.1 * Random.GetNumber();
        Input(1,i) = cos(i/13.678) + 0.1 * Random.GetNumber() + 0.5 * number;
        Input(2,i) = 3.3 * number + 0.1 * sin(i/13.678) + 1.2 * cos(i/13.678);
        Input(3,i) = 0.1 * number + 0.01 * Random.GetNumber();
        /*Input(0,i) = Sources.get_col(i) * Mix.get_row(0);
        Input(1,i) = Sources.get_col(i) * Mix.get_row(1);
        Input(2,i) = Sources.get_col(i) * Mix.get_row(2);
        Input(3,i) = Sources.get_col(i) * Mix.get_row(3);*/
        //Input(4,i) = RefData.GetData().Hx.data.at(i);
        //Input(5,i) = RefData.GetData().Hy.data.at(i);
    }
        
        itpp::Fast_ICA fastica(Input);
        fastica.set_nrof_independent_components(nchannels);
        //fastica.set_non_linearity(FICA_NONLIN_TANH);
        //fastica.set_approach(FICA_APPROACH_DEFL);
        fastica.separate();
        itpp::mat Output = fastica.get_independent_components();
        WriteToFile(Sources.get_row(0),Input.get_row(0),Output.get_row(0),"1");
        WriteToFile(Sources.get_row(1),Input.get_row(1),Output.get_row(1),"2");
        WriteToFile(Sources.get_row(2),Input.get_row(2),Output.get_row(2),"3");
        WriteToFile(Sources.get_row(3),Input.get_row(3),Output.get_row(3),"4");
}