MaxLloydQuantisierer: src/MaxLloydQuantizer.h Source File

00001 
00008 #ifndef MAXLLOYDQUANTISIERER_H_
00009 #define MAXLLOYDQUANTISIERER_H_
00010 
00011 #include <cmath>
00012 #include <itpp/itbase.h>
00013 #include <itpp/stat/histogram.h>
00014 
00015 using namespace std;
00016 
00017 namespace itpp {
00018 
00023 template<class T>class MaxLloydQuantizer {
00024 public:
00025 
00029         MaxLloydQuantizer() { initialized = trained = false; }
00030 
00035         MaxLloydQuantizer(unsigned int R) { setTargetRate(R); }
00036 
00045         MaxLloydQuantizer(unsigned int R,itpp::Vec<T> trainingdata,double epsilon=1e-3,unsigned int max_iter=1e3,unsigned hist_size=1000) { setTargetRate(R); train(trainingdata,epsilon,max_iter,hist_size); }
00046 
00050         virtual ~MaxLloydQuantizer() {}
00051 
00059         bool train(const itpp::Vec<T> &trainingdata,double epsilon=1e-3,unsigned int max_iter=1e3);
00060 
00061         bool train2(const itpp::Vec<T> &trainingdata,double epsilon=1e-3,unsigned int max_iter=1e3);
00062 
00067         unsigned int getNumberOfQuantizerLevels() const { return N; }
00068 
00073         unsigned int getTargetRate() const { return Rate; }
00074 
00079         bool isTrained() const { return trained; }
00080 
00085         bool isInitialized() const { return initialized; }
00086 
00092         void setTargetRate(unsigned int R)
00093         {
00094                 Rate=R;
00095                 N=1<<R;
00096                 d.set_length(N+1);
00097                 r.set_length(N);
00098                 trained=false;
00099                 initialized=true;
00100         }
00101 
00109         unsigned int quantize(T x);
00110 
00116         T inv_quant(unsigned int level);
00117 
00118 protected:
00119         itpp::Vec<T> d; 
00120         itpp::Vec<T> r; 
00121         unsigned int Rate; 
00122         unsigned int N; 
00123         bool trained; 
00124         bool initialized; 
00125 };
00126 
00127 
00128 template<class T> T MaxLloydQuantizer<T>::inv_quant(unsigned int level)
00129 {
00130 #ifdef DEBUG
00131         if(!initialized) cerr << "MaxLloydQuantizer<T>::inv_quant() called before initializing the quantizer!" << endl;
00132         if(!trained) cerr << "MaxLloydQuantizer<T>::inv_quant() called before training the quantizer!" << endl;
00133         if(level >= N) if(!trained) cerr << "MaxLloydQuantizer<T>::inv_quant(level): level out of range!" << endl;
00134 #endif
00135         return r(level);
00136 }
00137 
00138 
00139 /*template<class T>  bool MaxLloydQuantizer<T>::train2(const itpp::Vec<T> &trainingdata,double epsilon,unsigned int max_iter,unsigned int hist_size)
00140 {
00141 #ifdef DEBUG
00142         if(!initialized) cerr << "MaxLloydQuantizer<T>::train() called before initializing the quantizer!" << endl;
00143 #endif
00144         T minv,maxv;
00145         minv=maxv=trainingdata(0);
00146         for(const T* ptr=trainingdata._data();ptr < trainingdata._data()+trainingdata.length();ptr++)
00147         {
00148                 register T val = *ptr;
00149                 minv = ((val < minv) ? val : minv);
00150                 maxv = ((val > maxv) ? val : maxv);
00151         }
00152 
00153         // berechne nun die Intervallgrenzen eines gleichverteilten Quantisierers
00154         T step=(maxv-minv)/N;
00155         T w=minv;
00156         for(unsigned int k=0;k<=N;k++,w+=step) d(k)=w;
00157 
00158         // nun erstelle die Wahrscheinlichkeitsdichteverteilung
00159         itpp::Histogram<double> hist(minv,maxv,hist_size);
00160         hist.update(trainingdata);
00161         Vec<T> p=hist.get_pdf();
00162 
00163         // nun erfolgt der eigentliche Max-Lloyd-Algorithmus
00164         step = (maxv-minv)/hist_size;
00165         unsigned int iteration=0;
00166         for(;iteration < max_iter;iteration++)
00167         {
00168                 // berechne Rekonstruktionswerte neu
00169                 for(unsigned int i=0;i<N;i++)
00170                 {
00171                         unsigned int llimit = std::floor((d(i)-minv)/step);
00172                         unsigned int ulimit = std::floor((d(i+1)-minv)/step);
00173                         T upperint=0;
00174                         T lowerint=0;
00175                         T x=d(i); //minv+llimit*step;
00176                         for(unsigned int k=llimit;k<ulimit;k++,x+=step)
00177                         {
00178                                 T val=p(k);
00179                                 lowerint += val;
00180                                 upperint += val*x;
00181                         }
00182                         r(i)=upperint/lowerint;
00183                 }
00184                 // berechne Intervallgrenzen neu
00185                 double change=0.0;
00186                 for(unsigned int i=1;i<N;i++)
00187                 {
00188                         T di_neu=(r(i-1)+r(i))/2;
00189                         change += std::abs(double(d(i)-di_neu));
00190                         d(i)=di_neu;
00191                 }
00192                 if(change/N <= epsilon) break;
00193         }
00194 
00195         trained=true;
00196         return (iteration==max_iter);
00197 }*/
00198 
00199 template<class T>  bool MaxLloydQuantizer<T>::train(const itpp::Vec<T> &trainingdata,double epsilon,unsigned int max_iter)
00200 {
00201 #ifdef DEBUG
00202         if(!initialized) cerr << "MaxLloydQuantizer<T>::train() called before initializing the quantizer!" << endl;
00203 #endif
00204         if(trainingdata.length()==0) return false;
00205         T minv,maxv;
00206         minv=maxv=trainingdata(0);
00207         for(const T* ptr=trainingdata._data();ptr < trainingdata._data()+trainingdata.length();ptr++)
00208         {
00209                 register T val = *ptr;
00210                 minv = ((val < minv) ? val : minv);
00211                 maxv = ((val > maxv) ? val : maxv);
00212         }
00213 
00214         // berechne nun die Intervallgrenzen eines gleichverteilten Quantisierers
00215         T step=(maxv-minv)/N;
00216         T w=minv;
00217         for(unsigned int k=0;k<=N;k++,w+=step) d(k)=w;
00218         w=minv+step/2;
00219         for(unsigned int k=0;k<N;k++,w+=step) r(k)=w;
00220 
00221         itpp::Vec<T> r_neu(N);
00222         itpp::Vec<unsigned int> r_anz(N);
00223 
00224         unsigned int iteration=0;
00225         for(;iteration < max_iter;iteration++)
00226         {
00227                 r_neu.clear();
00228                 r_anz.clear();
00229                 // ordne jedem Trainingswert einen Rekonstruktionswert zu
00230                 for(int k=0;k<trainingdata.length();k++)
00231                 {
00232                         T val=trainingdata(k);
00233                         T min_dist=std::abs(r(0)-val);
00234                         unsigned min_index=0;
00235                         for(unsigned int l=1;l<N;l++)
00236                         {
00237                                 T dist=std::abs(T(r(l)-val));
00238                                 if(dist < min_dist)
00239                                 {
00240                                         min_dist=dist;
00241                                         min_index=l;
00242                                 }
00243                         }
00244                         r_neu(min_index)+=val;
00245                         r_anz(min_index)++;
00246                 }
00247                 double change=0.0;
00248                 for(unsigned int k=0;k<N;k++)
00249                 {
00250                         double r_neuw = (r_anz(k) != 0 ? r_neu(k)/r_anz(k) : r(k));
00251                         double tmp=r_neuw-r(k);
00252                         change+=tmp*tmp;
00253                         r(k)=r_neuw;
00254                 }
00255 
00256                 //cout << "Iteration " << iteration << " var: " << change << endl <<  " r=" << r << endl;
00257 
00258                 if(change/N <= epsilon) break;
00259         }
00260         // berechne Intervallgrenzen neu
00261 
00262         for(unsigned int i=1;i<N;i++)
00263         {
00264                 T di_neu=(r(i-1)+r(i))/2;
00265 
00266                 d(i)=di_neu;
00267         }
00268         trained=true;
00269         return (iteration!=max_iter);
00270 }
00271 
00272 template<class T>  bool MaxLloydQuantizer<T>::train2(const itpp::Vec<T> &trainingdata,double epsilon,unsigned int max_iter)
00273 {
00274 #ifdef DEBUG
00275         if(!initialized) cerr << "MaxLloydQuantizer<T>::train() called before initializing the quantizer!" << endl;
00276 #endif
00277         if(trainingdata.length()==0) return false;
00278         T minv,maxv;
00279         minv=maxv=trainingdata(0);
00280         for(const T* ptr=trainingdata._data();ptr < trainingdata._data()+trainingdata.length();ptr++)
00281         {
00282                 register T val = *ptr;
00283                 minv = ((val < minv) ? val : minv);
00284                 maxv = ((val > maxv) ? val : maxv);
00285         }
00286 
00287         // berechne nun die Intervallgrenzen eines gleichverteilten Quantisierers
00288         T step=(maxv-minv)/N;
00289         T w=minv;
00290         for(unsigned int k=0;k<N;k++,w+=step) d(k)=w;
00291         // zur Sicherheit
00292         d(N)=maxv;
00293         w=minv+step/2;
00294         for(unsigned int k=0;k<N;k++,w+=step) r(k)=w;
00295 
00296         itpp::Vec<T> r_neu(N);
00297         itpp::Vec<unsigned int> r_anz(N);
00298 
00299         unsigned int iteration=0;
00300         for(;iteration < max_iter;iteration++)
00301         {
00302                 r_neu.clear();
00303                 r_anz.clear();
00304                 // ordne jedem Trainingswert einen Rekonstruktionswert zu
00305                 for(int k=0;k<trainingdata.length();k++)
00306                 {
00307                         T val=trainingdata(k);
00308                         for(unsigned int l=1;l<=N;l++)
00309                         {
00310                                 if(val <= d(l))
00311                                 {
00312                                         r_neu(l-1)+=val;
00313                                         ++r_anz(l-1);
00314                                         break;
00315                                 }
00316                         }
00317                 }
00318                 double change=0.0;
00319                 for(unsigned int k=0;k<N;k++)
00320                 {
00321                         r(k) = (r_anz(k) != 0 ? r_neu(k)/r_anz(k) : 0.5*(d(k)+d(k+1)));
00322                 }
00323 
00324                 // berechne Intervallgrenzen neu
00325                 for (unsigned int i = 1; i < N; i++) {
00326                         T di_neu = (r(i - 1) + r(i)) / 2;
00327                         change+=std::abs((double) di_neu-d(i));
00328                         d(i) = di_neu;
00329                 }
00330                 cout << "Iteration " << iteration << " var: " << change << endl <<  " r=" << r << endl;
00331 
00332                 if(change/N <= epsilon) break;
00333         }
00334 
00335         trained=true;
00336         return (iteration!=max_iter);
00337 }
00338 
00339 
00340 template<class T> unsigned int MaxLloydQuantizer<T>::quantize(T x)
00341 {
00342 #ifdef DEBUG
00343         if(!initialized) cerr << "MaxLloydQuantizer<T>::quant() called before initializing the quantizer!" << endl;
00344         if(!trained) cerr << "MaxLloydQuantizer<T>::quant() called before training the quantizer!" << endl;
00345 #endif
00346         unsigned int k=0;
00347         while(++k != N) if(x < d(k)) break;
00348         return k-1;
00349 }
00350 
00351 }
00352 #endif /* MAXLLOYDQUANTISIERER_H_ */
00353