GPLIB++
FastICA.h
Go to the documentation of this file.
1 #ifndef FASTICA_H_
2 #define FASTICA_H_
3 #include "VecMat.h"
4 #include <complex>
5 #include <iostream>
6 #include <cmath>
7 #include "PCA.h"
8 
9 namespace ublas = boost::numeric::ublas;
10 
11 namespace gplib
12  {
13  /** \addtogroup sigproc Signal processing methods */
14  /* @{ */
15 
16  void FastICA(rmat &input, rmat &source_estimate, rmat &mixing_matrix)
17  {
18  mixing_matrix = ublas::identity_matrix<double>(input.size1());
19  rmat old_mix(mixing_matrix);
20  rmat new_mix(mixing_matrix);
21  rmat white_mat(mixing_matrix);
22  cmat evec(input.size1(), input.size1());
23  cvec eval(input.size1());
24  PCA(input, evec, eval);
25  white_mat = real(WhiteMat(evec, eval));
26  std::cout << "white_mat: " << white_mat << std::endl;
27  input = prod(white_mat, input);
28  const int maxiterations = 1000;
29  const double min_improv = 0.01;
30  double improvement = 10.0;
31  int niter = 0;
32  while (niter < maxiterations)
33  {
34  //std::cout << std::endl << std::endl << "Iter: " << niter << std::endl;
35  ublas::noalias(source_estimate) = ublas::prod(mixing_matrix, input);
36  rmat pow3mat = prod(trans(input), mixing_matrix);
37 
38  for (size_t i = 0; i < pow3mat.size1(); ++i)
39  {
40  for (size_t j = 0; j < pow3mat.size2(); ++j)
41  pow3mat(i, j) = std::pow(pow3mat(i, j), 3);
42  }
43  new_mix = prod(input, pow3mat);
44  new_mix /= input.size2();
45  new_mix -= 3.0 * mixing_matrix;
46  mixing_matrix = new_mix;
47 
48  rmat wwt = prod(mixing_matrix, trans(mixing_matrix));
49  mixing_matrix /= sqrt(norm_inf(wwt));
50  for (int i = 0; i < maxiterations; ++i)
51  {
52  rmat neww(mixing_matrix);
53  wwt = prod(mixing_matrix, trans(mixing_matrix));
54  neww *= 1.5;
55  mixing_matrix = prod(wwt, mixing_matrix);
56  mixing_matrix *= 0.5;
57  neww -= mixing_matrix;
58  mixing_matrix = neww;
59  //std::cout << "neww: " << neww << std::endl;
60  }
61 
62  ++niter;
63  }
64 
65  noalias(source_estimate) = prod(trans(mixing_matrix), input);
66  mixing_matrix = prod(trans(mixing_matrix), white_mat);
67  }
68  /* @} */
69  }
70 //THIS IS FOR THE FUTURE
71 /*
72  class FastICA
73  {
74  private:
75  boost::numeric::ublas::matrix<double> Weights;
76  public:
77  void Separate(const boost::numeric::vector<double> &In,boost::numeric::vector<double> &Out);
78  FastICA();
79  virtual ~FastICA();
80  };*/
81 
82 #endif /*FASTICA_H_*/
gplib::cmat WhiteMat(gplib::cmat &evectors, gplib::cvec &evalues)
Calculate the Whitening Matrix.
Definition: PCA.h:46
void FastICA(rmat &input, rmat &source_estimate, rmat &mixing_matrix)
Definition: FastICA.h:16
void PCA(const UblasMatrix &observations, gplib::cmat &evectors, gplib::cvec &evalues)
This template function calculates the principal component rotation matrix from a matrix of observatio...
Definition: PCA.h:25