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 #include <limits> 00043 00044 00045 namespace itpp { 00046 00050 typedef signed int QLLR; 00051 00055 typedef Vec<QLLR> QLLRvec; 00056 00060 typedef Mat<QLLR> QLLRmat; 00061 00065 const QLLR QLLR_MAX = (std::numeric_limits<QLLR>::max() >> 2); 00066 // added some margin to make sure the sum of two LLR is still permissible 00067 00117 class LLR_calc_unit { 00118 public: 00120 LLR_calc_unit(); 00121 00126 LLR_calc_unit(short int Dint1, short int Dint2, short int Dint3); 00127 00144 void init_llr_tables(short int Dint1 = 12, short int Dint2 = 300, 00145 short int Dint3 = 7); 00146 // void init_llr_tables(const short int d1=10, const short int d2=300, const short int d3=5); 00147 00149 inline QLLR to_qllr(const double &l); 00150 00152 QLLRvec to_qllr(const vec &l); 00153 00155 QLLRmat to_qllr(const mat &l); 00156 00158 inline double to_double(const QLLR &l) const { return ( ((double) l) / ((double) (1<<Dint1))); }; 00159 00161 vec to_double(const QLLRvec &l); 00162 00164 mat to_double(const QLLRmat &l); 00165 00170 inline QLLR jaclog(QLLR a, QLLR b); 00171 // Note: a version of this function taking "double" values as input 00172 // is deliberately omitted, because this is rather slow. 00173 00180 inline QLLR Boxplus(QLLR a, QLLR b); 00181 00185 QLLR logexp(QLLR x); 00186 00188 ivec get_Dint(); 00189 00191 void operator=(const LLR_calc_unit&); 00192 00194 friend std::ostream &operator<<(std::ostream &os, const LLR_calc_unit &lcu); 00195 00196 private: 00198 ivec construct_logexp_table(); 00199 00201 ivec logexp_table; 00202 00204 short int Dint1, Dint2, Dint3; 00205 00206 }; 00207 00211 std::ostream &operator<<(std::ostream &os, const LLR_calc_unit &lcu); 00212 00213 00214 00215 // ================ IMPLEMENTATIONS OF SOME LIKELIHOOD ALGEBRA FUNCTIONS =========== 00216 00217 inline QLLR LLR_calc_unit::to_qllr(const double &l) 00218 { 00219 it_assert0(l<=to_double(QLLR_MAX),"LLR_calc_unit::to_qllr(): overflow"); 00220 it_assert0(l>=-to_double(QLLR_MAX),"LLR_calc_unit::to_qllr(): overflow"); 00221 return ( (int) std::floor(0.5+(1<<Dint1)*l) ); 00222 }; 00223 00224 inline QLLR LLR_calc_unit::Boxplus(QLLR a, QLLR b) 00225 { 00226 /* These boundary checks are meant to completely eliminate the 00227 possibility of any numerical problem. Perhaps they are not 00228 strictly necessary. */ 00229 if (a>=QLLR_MAX) { 00230 return b; 00231 } 00232 if (b>=QLLR_MAX) { 00233 return a; 00234 } 00235 if (a<=-QLLR_MAX) { 00236 return (-b); 00237 } 00238 if (b<=-QLLR_MAX) { 00239 return (-a); 00240 } 00241 00242 QLLR a_abs = (a>0 ? a : -a); 00243 QLLR b_abs = (b>0 ? b : -b); 00244 QLLR minabs = (a_abs>b_abs ? b_abs : a_abs); 00245 QLLR term1 = (a>0 ? (b>0 ? minabs : -minabs) : (b>0 ? -minabs : minabs)); 00246 QLLR apb = a+b; 00247 QLLR term2 = logexp((apb>0 ? apb : -apb)); 00248 QLLR amb = a-b; 00249 QLLR term3 = logexp((amb>0 ? amb : -amb)); 00250 00251 QLLR result = term1+term2-term3; 00252 if (result>=QLLR_MAX) { 00253 return QLLR_MAX; 00254 } else if (result<=-QLLR_MAX) { 00255 return (-QLLR_MAX); 00256 } else { 00257 return result; 00258 } 00259 } 00260 00261 inline QLLR LLR_calc_unit::logexp(QLLR x) 00262 { 00263 it_assert0(x>=0,"LLR_calc_unit::logexp() is not defined for negative LLR values"); 00264 int ind = x>>Dint3; 00265 if (ind>=Dint2) { 00266 return 0; 00267 } 00268 00269 it_assert0(ind>=0,"LLR_calc_unit::logexp() internal error"); 00270 it_assert0(ind<Dint2,"LLR_calc_unit::logexp() internal error"); 00271 00272 // With interpolation 00273 // int delta=x-(ind<<Dint3); 00274 // return ((delta*logexp_table(ind+1) + ((1<<Dint3)-delta)*logexp_table(ind)) >> Dint3); 00275 00276 // Without interpolation 00277 return logexp_table(ind); 00278 } 00279 00280 inline QLLR LLR_calc_unit::jaclog(QLLR a, QLLR b ) 00281 { 00282 QLLR x,maxab; 00283 00284 if (a>b) { 00285 maxab = a; 00286 x=a-b; 00287 } else { 00288 maxab = b; 00289 x=b-a; 00290 } 00291 00292 if (maxab>=QLLR_MAX) { 00293 return QLLR_MAX; 00294 } else { 00295 return (maxab + logexp(x)); 00296 } 00297 }; 00298 00299 } 00300 00301 #endif
Generated on Fri Jun 8 01:07:12 2007 for IT++ by Doxygen 1.5.2