testlp.cpp

Go to the documentation of this file.

#include <iostream>
#include "WienerInterpolator.h"
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <algorithm>
#include <iterator>
#include "UniformRNG.h"
#include <itpp/itbase.h>
#include <itpp/srccode/lpcfunc.h>
#include <fstream>
#include <numeric>
#include <nr.h>
#include <nrutil.h>
#include <gsl/gsl_math.h>

namespace ublas = boost::numeric::ublas;
int main()
{
        UniformRNG Random;
        
        const int datasize = Random.GetNumber(1000);
        const int filterlength = std::min(datasize/(Random.GetNumber(10)+1)+1,datasize);
        const int inputlength = std::min(datasize,filterlength*(Random.GetNumber(10)+1));
        //const int datasize =4;
        //const int filterlength = 2;
        //const int inputlength = 4;
        float *nrdata;
        float *d;
        nrdata = vector(1,inputlength);
        d = vector(1,filterlength);
        float xms;
        ublas::vector<float> Data(datasize);
        itpp::vec ItppData(inputlength);
        itpp::vec ItppCoeff(filterlength);
        itpp::vec StabCoeff(filterlength);
        
        for (int i = 0; i < datasize; ++i)
        {
                Data(i) = sin(static_cast<double>(i)/1.2)+ (0.5 - Random.GetNumber())*0.1;
                
        }
        WienerInterpolator Wiener(filterlength);
        ublas::vector<float> Input(inputlength);
        ublas::vector<double> Prediction(datasize);
        for (int i = 0; i < inputlength; ++i)
        {
                Input(i) = Data(i);
                ItppData(i) = Data(i);
                Prediction(i) = Data(i);
                nrdata[i+1] = Data(i);
        }
        Wiener.AdaptFilter(Input,Input);
        memcof(&(Input[0])-1,inputlength,filterlength,&xms,d);
        
        const double tolerance = 1e-3;
        std::cout << "Datasize: " << datasize << " Filterlength: " << filterlength << " Inputsize: " << inputlength<< std::endl; 
        for (int i = 0; i < filterlength; ++i)
        {
                if ( gsl_fcmp(Wiener.GetWeights()(i),d[i+1],tolerance) != 0)
                {
                        std::cout << "Index: " << i << " Own: " <<  Wiener.GetWeights()(i) << " NR: " << d[i+1] << std::endl;
                }
        }
        
        /*ItppCoeff = itpp::schurrc(ItppData,filterlength);
        polystab(ItppCoeff,StabCoeff);
        std::ofstream outfile("out");
        for (int i = inputlength - filterlength; i < datasize - filterlength; ++i)
        {
                //vector_range<vector<double> > vr(Data, range(i,i+filterlength));
                //Prediction(i+filterlength) = inner_prod(Wiener.GetWeights(),vr);
                Prediction(i+filterlength) = 0.0;
                for (int j = 0; j < filterlength; ++j)
                {
                        //Prediction(i+filterlength) += ItppCoeff(filterlength-j-1) * Prediction(i+j);
                        Prediction(i+filterlength) += d[filterlength-j] * Prediction(i+j);
                }
        }
         std::cout << "Own: " << std::endl;
        for (int i = 0; i < datasize; ++i)
           outfile << i << " " << Data(i) << " " << Prediction(i) <<  std::endl;
         copy(Wiener.GetWeights().rbegin(),Wiener.GetWeights().rend(),std::ostream_iterator<double>(std::cout," "));
         std::cout << std::endl;
        
        std::cout << "NR: " << std::endl;
         for (int i = 0; i < filterlength; ++i)
                std::cout << d[i+1] << " ";
         std::cout << std::endl;
         
         std::cout << "Itpp: " << std::endl;
         for (int i = 0; i < filterlength; ++i)
                std::cout << ItppCoeff(i) << " ";
          std::cout << std::endl;
          
          /*for (int i = 0; i < filterlength; ++i)
                std::cout << StabCoeff(i) << " ";
          std::cout << std::endl;*/
}

---