00001 00035 #ifndef LLR_H 00036 #define LLR_H 00037 00038 #include <itpp/base/vec.h> 00039 #include <itpp/base/mat.h> 00040 #include <itpp/base/specmat.h> 00041 #include <itpp/base/matfunc.h> 00042 00043 00044 namespace itpp { 00045 00049 typedef signed int QLLR; 00050 00054 typedef Vec<QLLR> QLLRvec; 00055 00059 typedef Mat<QLLR> QLLRmat; 00060 00064 const QLLR QLLR_MAX=(INT_MAX>>2) ; 00065 // added some margin to make sure the sum of two LLR is still permissible 00066 00116 class LLR_calc_unit { 00117 public: 00119 LLR_calc_unit(); 00120 00125 LLR_calc_unit(short int Dint1, short int Dint2, short int Dint3); 00126 00143 void init_llr_tables(short int Dint1 = 12, short int Dint2 = 300, 00144 short int Dint3 = 7); 00145 // void init_llr_tables(const short int d1=10, const short int d2=300, const short int d3=5); 00146 00148 inline QLLR to_qllr(const double &l); 00149 00151 QLLRvec to_qllr(const vec &l); 00152 00154 QLLRmat to_qllr(const mat &l); 00155 00157 inline double to_double(const QLLR &l) const { return ( ((double) l) / ((double) (1<<Dint1))); }; 00158 00160 vec to_double(const QLLRvec &l); 00161 00163 mat to_double(const QLLRmat &l); 00164 00169 inline QLLR jaclog(QLLR a, QLLR b); 00170 // Note: a version of this function taking "double" values as input 00171 // is deliberately omitted, because this is rather slow. 00172 00179 inline QLLR Boxplus(QLLR a, QLLR b); 00180 00184 QLLR logexp(QLLR x); 00185 00187 ivec get_Dint(); 00188 00190 void operator=(const LLR_calc_unit&); 00191 00193 friend std::ostream &operator<<(std::ostream &os, const LLR_calc_unit &lcu); 00194 00195 private: 00197 ivec construct_logexp_table(); 00198 00200 ivec logexp_table; 00201 00203 short int Dint1, Dint2, Dint3; 00204 00205 }; 00206 00210 std::ostream &operator<<(std::ostream &os, const LLR_calc_unit &lcu); 00211 00212 00213 00214 // ================ IMPLEMENTATIONS OF SOME LIKELIHOOD ALGEBRA FUNCTIONS =========== 00215 00216 inline QLLR LLR_calc_unit::to_qllr(const double &l) 00217 { 00218 it_assert0(l<=to_double(QLLR_MAX),"LLR_calc_unit::to_qllr(): overflow"); 00219 it_assert0(l>=-to_double(QLLR_MAX),"LLR_calc_unit::to_qllr(): overflow"); 00220 return ( (int) std::floor(0.5+(1<<Dint1)*l) ); 00221 }; 00222 00223 inline QLLR LLR_calc_unit::Boxplus(QLLR a, QLLR b) 00224 { 00225 /* These boundary checks are meant to completely eliminate the 00226 possibility of any numerical problem. Perhaps they are not 00227 strictly necessary. */ 00228 if (a>=QLLR_MAX) { 00229 return b; 00230 } 00231 if (b>=QLLR_MAX) { 00232 return a; 00233 } 00234 if (a<=-QLLR_MAX) { 00235 return (-b); 00236 } 00237 if (b<=-QLLR_MAX) { 00238 return (-a); 00239 } 00240 00241 QLLR a_abs = (a>0 ? a : -a); 00242 QLLR b_abs = (b>0 ? b : -b); 00243 QLLR minabs = (a_abs>b_abs ? b_abs : a_abs); 00244 QLLR term1 = (a>0 ? (b>0 ? minabs : -minabs) : (b>0 ? -minabs : minabs)); 00245 QLLR apb = a+b; 00246 QLLR term2 = logexp((apb>0 ? apb : -apb)); 00247 QLLR amb = a-b; 00248 QLLR term3 = logexp((amb>0 ? amb : -amb)); 00249 00250 QLLR result = term1+term2-term3; 00251 if (result>=QLLR_MAX) { 00252 return QLLR_MAX; 00253 } else if (result<=-QLLR_MAX) { 00254 return (-QLLR_MAX); 00255 } else { 00256 return result; 00257 } 00258 } 00259 00260 inline QLLR LLR_calc_unit::logexp(QLLR x) 00261 { 00262 it_assert0(x>=0,"LLR_calc_unit::logexp() is not defined for negative LLR values"); 00263 int ind = x>>Dint3; 00264 if (ind>=Dint2) { 00265 return 0; 00266 } 00267 00268 it_assert0(ind>=0,"LLR_calc_unit::logexp() internal error"); 00269 it_assert0(ind<Dint2,"LLR_calc_unit::logexp() internal error"); 00270 00271 // With interpolation 00272 // int delta=x-(ind<<Dint3); 00273 // return ((delta*logexp_table(ind+1) + ((1<<Dint3)-delta)*logexp_table(ind)) >> Dint3); 00274 00275 // Without interpolation 00276 return logexp_table(ind); 00277 } 00278 00279 inline QLLR LLR_calc_unit::jaclog(QLLR a, QLLR b ) 00280 { 00281 QLLR x,maxab; 00282 00283 if (a>b) { 00284 maxab = a; 00285 x=a-b; 00286 } else { 00287 maxab = b; 00288 x=b-a; 00289 } 00290 00291 if (maxab>=QLLR_MAX) { 00292 return QLLR_MAX; 00293 } else { 00294 return (maxab + logexp(x)); 00295 } 00296 }; 00297 00298 } 00299 00300 #endif
Generated on Thu Apr 19 14:23:58 2007 for IT++ by Doxygen 1.4.6