00001 00029 #include <itpp/comm/exit.h> 00030 #include <itpp/stat/histogram.h> //histogram class for mutual information computation 00031 #include <itpp/base/itcompat.h> 00032 00033 namespace itpp 00034 { 00035 00036 double EXIT::sigma2A;//allocate memory for static member variable 00037 00038 double EXIT::extrinsic_mutual_info(const itpp::vec &obs, const itpp::bvec &cond, const int &N) 00039 { 00040 //initialize histogram 00041 itpp::Histogram<double> hist(itpp::min(obs), itpp::max(obs), N);//common definition interval for both PDFs 00042 00043 //conditional PDF knowing that a bit of 0 was emitted 00044 itpp::ivec idx = itpp::find(cond==itpp::bin(0)); 00045 itpp::vec cond_obs = obs(idx); 00046 hist.reset();//start counting 00047 hist.update(cond_obs); 00048 itpp::vec left_pdf = hist.get_pdf();//the pdf is computed without taking into account the interval length (step) 00049 itpp::ivec left_int = itpp::find(left_pdf!=0);//integration interval for the left PDF 00050 00051 //conditional PDF knowing that a bit of 1 was emitted 00052 idx = itpp::find(cond==itpp::bin(1)); 00053 cond_obs = obs(idx); 00054 hist.reset();//restart counting 00055 hist.update(cond_obs); 00056 itpp::vec right_pdf = hist.get_pdf(); 00057 itpp::ivec right_int = itpp::find(right_pdf!=0);//integration interval for the right PDF 00058 00059 //mutual extrinsic information 00060 itpp::vec left_half = itpp::elem_mult(left_pdf(left_int), itpp::log2(itpp::elem_div(2.0*left_pdf(left_int), left_pdf(left_int)+right_pdf(left_int)))); 00061 double IE = itpp::sum(left_half)-0.5*(left_half(0)+left_half(left_half.length()-1));//numerical integration without taking into account the inteval length (see conditional PDF computation) 00062 itpp::vec right_half = itpp::elem_mult(right_pdf(right_int), itpp::log2(itpp::elem_div(2.0*right_pdf(right_int), left_pdf(right_int)+right_pdf(right_int)))); 00063 IE += itpp::sum(right_half)-0.5*(right_half(0)+right_half(right_half.length()-1));//numerical integration 00064 IE *= 0.5; 00065 00066 return IE; 00067 } 00068 00069 double EXIT::gaussian_fct(const double x) 00070 { 00071 return (1.0/std::sqrt(sigma2A*itpp::m_2pi))*std::exp(-itpp::sqr(x-(sigma2A/2.0))/(2.0*sigma2A))*::log2(1+std::exp(-x)); 00072 }; 00073 00074 }//namespace itpp 00075
Generated on Wed Jul 27 2011 16:27:04 for IT++ by Doxygen 1.7.4