00001 00030 #ifndef LLR_H 00031 #define LLR_H 00032 00033 #include <limits> 00034 #include <itpp/base/vec.h> 00035 #include <itpp/base/mat.h> 00036 #include <itpp/base/specmat.h> 00037 #include <itpp/base/matfunc.h> 00038 #include <limits> 00039 00040 namespace itpp { 00041 00045 typedef signed int QLLR; 00046 00050 typedef Vec<QLLR> QLLRvec; 00051 00055 typedef Mat<QLLR> QLLRmat; 00056 00060 const QLLR QLLR_MAX = (std::numeric_limits<QLLR>::max() >> 4); 00061 // added some margin to make sure the sum of two LLR is still permissible 00062 00112 class LLR_calc_unit { 00113 public: 00115 LLR_calc_unit(); 00116 00121 LLR_calc_unit(short int Dint1, short int Dint2, short int Dint3); 00122 00150 void init_llr_tables(short int Dint1 = 12, short int Dint2 = 300, 00151 short int Dint3 = 7); 00152 // void init_llr_tables(const short int d1=10, const short int d2=300, 00153 // const short int d3=5); 00154 00156 inline QLLR to_qllr(const double &l) const; 00157 00159 QLLRvec to_qllr(const vec &l) const; 00160 00162 QLLRmat to_qllr(const mat &l) const; 00163 00165 inline double to_double(const QLLR &l) const { return ( ((double) l) / ((double) (1<<Dint1))); }; 00166 00168 vec to_double(const QLLRvec &l) const; 00169 00171 mat to_double(const QLLRmat &l) const; 00172 00177 inline QLLR jaclog(QLLR a, QLLR b) const; 00178 // Note: a version of this function taking "double" values as input 00179 // is deliberately omitted, because this is rather slow. 00180 00187 inline QLLR Boxplus(QLLR a, QLLR b) const; 00188 00192 inline QLLR logexp(QLLR x) const; 00193 00195 ivec get_Dint(); 00196 00198 void operator=(const LLR_calc_unit&); 00199 00201 friend std::ostream &operator<<(std::ostream &os, const LLR_calc_unit &lcu); 00202 00203 private: 00205 ivec construct_logexp_table(); 00206 00208 ivec logexp_table; 00209 00211 short int Dint1, Dint2, Dint3; 00212 00213 }; 00214 00219 std::ostream &operator<<(std::ostream &os, const LLR_calc_unit &lcu); 00220 00221 // ================ IMPLEMENTATIONS OF SOME LIKELIHOOD ALGEBRA FUNCTIONS =========== 00222 00223 inline QLLR LLR_calc_unit::to_qllr(const double &l) const 00224 { 00225 it_assert_debug(l<=to_double(QLLR_MAX),"LLR_calc_unit::to_qllr(): overflow"); 00226 it_assert_debug(l>=-to_double(QLLR_MAX),"LLR_calc_unit::to_qllr(): overflow"); 00227 return ( (int) std::floor(0.5+(1<<Dint1)*l) ); 00228 }; 00229 00230 inline QLLR LLR_calc_unit::Boxplus(QLLR a, QLLR b) const 00231 { 00232 it_assert_debug(a<QLLR_MAX,"LLR_calc_unit::Boxplus(): owerflow"); 00233 it_assert_debug(b<QLLR_MAX,"LLR_calc_unit::Boxplus(): owerflow"); 00234 it_assert_debug(a>-QLLR_MAX,"LLR_calc_unit::Boxplus(): owerflow"); 00235 it_assert_debug(b>-QLLR_MAX,"LLR_calc_unit::Boxplus(): owerflow"); 00236 00237 QLLR a_abs = (a>0 ? a : -a); 00238 QLLR b_abs = (b>0 ? b : -b); 00239 QLLR minabs = (a_abs>b_abs ? b_abs : a_abs); 00240 QLLR term1 = (a>0 ? (b>0 ? minabs : -minabs) : (b>0 ? -minabs : minabs)); 00241 00242 if (Dint2==0) { // logmax approximation - avoid looking into empty table 00243 it_assert_debug(term1<QLLR_MAX,"LLR_calc_unit::Boxplus(): owerflow"); 00244 it_assert_debug(term1>-QLLR_MAX,"LLR_calc_unit::Boxplus(): owerflow"); 00245 return term1; 00246 } 00247 00248 QLLR apb = a+b; 00249 QLLR term2 = logexp((apb>0 ? apb : -apb)); 00250 QLLR amb = a-b; 00251 QLLR term3 = logexp((amb>0 ? amb : -amb)); 00252 00253 QLLR result = term1+term2-term3; 00254 00255 it_assert_debug(result<QLLR_MAX,"LLR_calc_unit::Boxplus(): owerflow"); 00256 it_assert_debug(result>-QLLR_MAX,"LLR_calc_unit::Boxplus(): owerflow"); 00257 return result; 00258 } 00259 00260 inline QLLR LLR_calc_unit::logexp(QLLR x) const 00261 { 00262 it_assert_debug(x>=0,"LLR_calc_unit::logexp() is not defined for negative LLR values"); 00263 int ind = x>>Dint3; 00264 if (ind>=Dint2) { // outside table 00265 return 0; 00266 } 00267 00268 it_assert_debug(ind>=0,"LLR_calc_unit::logexp() internal error"); 00269 it_assert_debug(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 ) const 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 Sun Sep 14 18:52:35 2008 for IT++ by Doxygen 1.5.6