00001 00029 #include <itpp/srccode/vqtrain.h> 00030 #include <itpp/base/matfunc.h> 00031 #include <itpp/base/random.h> 00032 #include <itpp/base/sort.h> 00033 #include <itpp/base/specmat.h> 00034 #include <itpp/base/math/min_max.h> 00035 #include <itpp/stat/misc_stat.h> 00036 #include <iostream> 00037 00039 00040 namespace itpp 00041 { 00042 00043 // the cols contains codevectors 00044 double kmeansiter(Array<vec> &DB, mat &codebook) 00045 { 00046 int DIM = DB(0).length(), SIZE = codebook.cols(), T = DB.length(); 00047 vec x, xnum(SIZE); 00048 mat xsum(DIM, SIZE); 00049 int n, MinIndex, i, j, k; 00050 double MinS, S, D, Dold, *xp, *cp; 00051 00052 xsum.clear(); 00053 xnum.clear(); 00054 00055 n = 0; 00056 D = 1E20; 00057 Dold = D; 00058 D = 0; 00059 for (k = 0;k < T;k++) { 00060 x = DB(k); 00061 xp = x._data(); 00062 MinS = 1E20; 00063 MinIndex = 0; 00064 for (i = 0;i < SIZE;i++) { 00065 cp = &codebook(0, i); 00066 S = sqr(xp[0] - cp[0]); 00067 for (j = 1;j < DIM;j++) { 00068 S += sqr(xp[j] - cp[j]); 00069 if (S >= MinS) goto sune; 00070 } 00071 MinS = S; 00072 MinIndex = i; 00073 sune: 00074 i = i; 00075 } 00076 D += MinS; 00077 cp = &xsum(0, MinIndex); 00078 for (j = 0;j < DIM;j++) { 00079 cp[j] += xp[j]; 00080 } 00081 xnum(MinIndex)++; 00082 } 00083 for (i = 0;i < SIZE;i++) { 00084 for (j = 0;j < DIM;j++) { 00085 codebook(j, i) = xsum(j, i) / xnum(i); 00086 } 00087 } 00088 return D; 00089 } 00090 00091 mat kmeans(Array<vec> &DB, int SIZE, int NOITER, bool VERBOSE) 00092 { 00093 int DIM = DB(0).length(), T = DB.length(); 00094 mat codebook(DIM, SIZE); 00095 int n, i, j; 00096 double D, Dold; 00097 ivec ind(SIZE); 00098 00099 for (i = 0;i < SIZE;i++) { 00100 ind(i) = randi(0, T - 1); 00101 j = 0; 00102 while (j < i) { 00103 if (ind(j) == ind(i)) { 00104 ind(i) = randi(0, T - 1); 00105 j = 0; 00106 } 00107 j++; 00108 } 00109 codebook.set_col(i, DB(ind(i))); 00110 } 00111 00112 00113 if (VERBOSE) std::cout << "Training VQ..." << std::endl ; 00114 00115 D = 1E20; 00116 for (n = 0;n < NOITER;n++) { 00117 Dold = D; 00118 D = kmeansiter(DB, codebook); 00119 if (VERBOSE) std::cout << n << ": " << D / T << " "; 00120 if (std::abs((D - Dold) / D) < 1e-4) break; 00121 } 00122 return codebook; 00123 } 00124 00125 mat lbg(Array<vec> &DB, int SIZE, int NOITER, bool VERBOSE) 00126 { 00127 int S = 1, DIM = DB(0).length(), T = DB.length(), i, n; 00128 mat cb; 00129 vec delta = 0.001 * randn(DIM), x; 00130 double D, Dold; 00131 00132 x = zeros(DIM); 00133 for (i = 0;i < T;i++) { 00134 x += DB(i); 00135 } 00136 x /= T; 00137 cb.set_size(DIM, 1); 00138 cb.set_col(0, x); 00139 while (cb.cols() < SIZE) { 00140 S = cb.cols(); 00141 if (2*S <= SIZE) cb.set_size(DIM, 2*S, true); 00142 else cb.set_size(DIM, 2*S, true); 00143 for (i = S;i < cb.cols();i++) { 00144 cb.set_col(i, cb.get_col(i - S) + delta); 00145 } 00146 D = 1E20; 00147 for (n = 0;n < NOITER;n++) { 00148 Dold = D; 00149 D = kmeansiter(DB, cb); 00150 if (VERBOSE) std::cout << n << ": " << D / T << " "; 00151 if (std::abs((D - Dold) / D) < 1e-4) break; 00152 } 00153 } 00154 00155 return cb; 00156 } 00157 00158 mat vqtrain(Array<vec> &DB, int SIZE, int NOITER, double STARTSTEP, bool VERBOSE) 00159 { 00160 int DIM = DB(0).length(); 00161 vec x; 00162 vec codebook(DIM*SIZE); 00163 int n, MinIndex, i, j; 00164 double MinS, S, D, step, *xp, *cp; 00165 00166 for (i = 0;i < SIZE;i++) { 00167 codebook.replace_mid(i*DIM, DB(randi(0, DB.length() - 1))); 00168 } 00169 if (VERBOSE) std::cout << "Training VQ..." << std::endl ; 00170 00171 res: 00172 D = 0; 00173 for (n = 0;n < NOITER;n++) { 00174 step = STARTSTEP * (1.0 - double(n) / NOITER); 00175 if (step < 0) step = 0; 00176 x = DB(randi(0, DB.length() - 1)); // seems unnecessary! Check it up. 00177 xp = x._data(); 00178 00179 MinS = 1E20; 00180 MinIndex = 0; 00181 for (i = 0;i < SIZE;i++) { 00182 cp = &codebook(i * DIM); 00183 S = sqr(xp[0] - cp[0]); 00184 for (j = 1;j < DIM;j++) { 00185 S += sqr(xp[j] - cp[j]); 00186 if (S >= MinS) goto sune; 00187 } 00188 MinS = S; 00189 MinIndex = i; 00190 sune: 00191 i = i; 00192 } 00193 D += MinS; 00194 cp = &codebook(MinIndex * DIM); 00195 for (j = 0;j < DIM;j++) { 00196 cp[j] += step * (xp[j] - cp[j]); 00197 } 00198 if ((n % 20000 == 0) && (n > 1)) { 00199 if (VERBOSE) std::cout << n << ": " << D / 20000 << " "; 00200 D = 0; 00201 } 00202 } 00203 00204 // checking training result 00205 vec dist(SIZE), num(SIZE); 00206 00207 dist.clear(); 00208 num.clear(); 00209 for (n = 0;n < DB.length();n++) { 00210 x = DB(n); 00211 xp = x._data(); 00212 MinS = 1E20; 00213 MinIndex = 0; 00214 for (i = 0;i < SIZE;i++) { 00215 cp = &codebook(i * DIM); 00216 S = sqr(xp[0] - cp[0]); 00217 for (j = 1;j < DIM;j++) { 00218 S += sqr(xp[j] - cp[j]); 00219 if (S >= MinS) goto sune2; 00220 } 00221 MinS = S; 00222 MinIndex = i; 00223 sune2: 00224 i = i; 00225 } 00226 dist(MinIndex) += MinS; 00227 num(MinIndex) += 1; 00228 } 00229 dist = 10 * log10(dist * dist.length() / sum(dist)); 00230 if (VERBOSE) std::cout << std::endl << "Distortion contribution: " << dist << std::endl ; 00231 if (VERBOSE) std::cout << "Num spread: " << num / DB.length()*100 << " %" << std::endl << std::endl ; 00232 if (min(dist) < -30) { 00233 std::cout << "Points without entries! Retraining" << std::endl ; 00234 j = min_index(dist); 00235 i = max_index(num); 00236 codebook.replace_mid(j*DIM, codebook.mid(i*DIM, DIM)); 00237 goto res; 00238 } 00239 00240 mat cb(DIM, SIZE); 00241 for (i = 0;i < SIZE;i++) { 00242 cb.set_col(i, codebook.mid(i*DIM, DIM)); 00243 } 00244 return cb; 00245 } 00246 00247 vec sqtrain(const vec &inDB, int SIZE) 00248 { 00249 vec DB(inDB); 00250 vec Levels, Levels_old; 00251 ivec indexlist(SIZE + 1); 00252 int il, im, ih, i; 00253 int SIZEDB = inDB.length(); 00254 double x; 00255 00256 sort(DB); 00257 Levels = DB(round_i(linspace(0.01 * SIZEDB, 0.99 * SIZEDB, SIZE))); 00258 Levels_old = zeros(SIZE); 00259 00260 while (energy(Levels - Levels_old) > 0.0001) { 00261 Levels_old = Levels; 00262 for (i = 0;i < SIZE - 1;i++) { 00263 x = (Levels(i) + Levels(i + 1)) / 2; 00264 il = 0; 00265 ih = SIZEDB - 1; 00266 while (il < ih - 1) { 00267 im = (il + ih) / 2; 00268 if (x < DB(im)) ih = im; 00269 else il = im; 00270 } 00271 indexlist(i + 1) = il; 00272 } 00273 indexlist(0) = -1; 00274 indexlist(SIZE) = SIZEDB - 1; 00275 for (i = 0;i < SIZE;i++) Levels(i) = mean(DB(indexlist(i) + 1, indexlist(i + 1))); 00276 } 00277 return Levels; 00278 } 00279 00280 ivec bitalloc(const vec &variances, int nobits) 00281 { 00282 ivec bitvec(variances.length()); 00283 bitvec.clear(); 00284 int i, bits = nobits; 00285 vec var = variances; 00286 00287 while (bits > 0) { 00288 i = max_index(var); 00289 var(i) /= 4; 00290 bitvec(i)++; 00291 bits--; 00292 } 00293 return bitvec; 00294 } 00295 00296 } // namespace itpp 00297
Generated on Wed Jul 27 2011 16:27:05 for IT++ by Doxygen 1.7.4