00001 00033 #include <itpp/base/specmat.h> 00034 #include <itpp/base/elmatfunc.h> 00035 #include <itpp/base/matfunc.h> 00036 #include <itpp/base/stat.h> 00037 00038 00039 namespace itpp { 00040 00041 ivec find(const bvec &invector) 00042 { 00043 it_assert(invector.size()>0,"find(): vector cannot be empty"); 00044 ivec temp(invector.size()); 00045 int pos=0; 00046 for (int i=0;i<invector.size();i++) { 00047 if (invector(i)==bin(1)) { 00048 temp(pos)=i;pos++; 00049 } 00050 } 00051 temp.set_size(pos, true); 00052 return temp; 00053 } 00054 00055 #ifndef DOXYGEN_SHOULD_SKIP_THIS 00056 00057 #define CREATE_SET_FUNS(typef,typem,name,value) \ 00058 typef name(int size) \ 00059 { \ 00060 typef t(size); \ 00061 t = value; \ 00062 return t; \ 00063 } \ 00064 \ 00065 typem name(int rows, int cols) \ 00066 { \ 00067 typem t(rows, cols); \ 00068 t = value; \ 00069 return t; \ 00070 } 00071 00072 #define CREATE_EYE_FUN(type,name,zero,one) \ 00073 type name(int size) { \ 00074 type t(size,size); \ 00075 t = zero; \ 00076 for (int i=0; i<size; i++) \ 00077 t(i,i) = one; \ 00078 return t; \ 00079 } 00080 00081 CREATE_SET_FUNS(vec, mat, ones, 1.0) 00082 CREATE_SET_FUNS(bvec, bmat, ones_b, bin(1)) 00083 CREATE_SET_FUNS(ivec, imat, ones_i, 1) 00084 CREATE_SET_FUNS(cvec, cmat, ones_c, std::complex<double>(1.0)) 00085 00086 CREATE_SET_FUNS(vec, mat, zeros, 0.0) 00087 CREATE_SET_FUNS(bvec, bmat, zeros_b, bin(0)) 00088 CREATE_SET_FUNS(ivec, imat, zeros_i, 0) 00089 CREATE_SET_FUNS(cvec, cmat, zeros_c, std::complex<double>(0.0)) 00090 00091 CREATE_EYE_FUN(mat, eye, 0.0, 1.0) 00092 CREATE_EYE_FUN(bmat, eye_b, bin(0), bin(1)) 00093 CREATE_EYE_FUN(imat, eye_i, 0, 1) 00094 CREATE_EYE_FUN(cmat, eye_c, std::complex<double>(0.0), std::complex<double>(1.0)) 00095 00096 template <class T> 00097 void eye(int size, Mat<T> &m) 00098 { 00099 m.set_size(size, size, false); 00100 m = T(0); 00101 for (int i=size-1; i>=0; i--) 00102 m(i,i) = T(1); 00103 } 00104 00105 #endif //DOXYGEN_SHOULD_SKIP_THIS 00106 00107 vec impulse(int size) { 00108 vec t(size); 00109 t.clear(); 00110 t[0]=1.0; 00111 return t; 00112 } 00113 00114 vec linspace(double from, double to, int points) { 00115 if (points<2){ 00116 // This is the "Matlab definition" of linspace 00117 vec output(1); 00118 output(0)=to; 00119 return output; 00120 } 00121 else{ 00122 vec output(points); 00123 double step = (to - from) / double(points-1); 00124 int i; 00125 for (i=0; i<points; i++) 00126 output(i) = from + i*step; 00127 return output; 00128 } 00129 } 00130 00131 00132 // Construct a Hadamard-imat of size "size" 00133 imat hadamard(int size) { 00134 int i,k,l,pow2,logsize; 00135 imat H(size,size); 00136 logsize = levels2bits(size); 00137 00138 it_assert1(pow2i(logsize)==size,"hadamard size not a power of 2"); 00139 H(0,0)=1;H(0,1)=1;H(1,0)=1;H(1,1)=-1; 00140 00141 for (i=1;i<logsize;i++) { 00142 // pow2=round_i(pow(2,i)); // Unbeliveably slow 00143 pow2 = 1<<i; 00144 for (k=0;k<pow2;k++) { 00145 for (l=0;l<pow2;l++) { 00146 H(k,l)=H(k,l); 00147 H(k+pow2,l)=H(k,l); 00148 H(k,l+pow2)=H(k,l); 00149 H(k+pow2,l+pow2)=(-1)*H(k,l); 00150 } 00151 } 00152 } 00153 return H; 00154 } 00155 00156 imat jacobsthal(int p) 00157 { 00158 int quadratic_residue; 00159 imat out(p,p); 00160 int i, j; 00161 00162 out = -1; // start with all elements equal to "-1" 00163 00164 // Generate a complete list of quadratic residues 00165 for (i=0; i<(p-1)/2; i++) { 00166 quadratic_residue=((i+1)*(i+1))%p; 00167 // set this element in all rows (col-row) = quadratic_residue 00168 for (j=0; j<p; j++) { 00169 out(j, (j+quadratic_residue)%p)=1; 00170 } 00171 } 00172 00173 // set diagonal elements to zero 00174 for (i=0; i<p; i++) { 00175 out(i,i)=0; 00176 } 00177 return out; 00178 } 00179 00180 imat conference(int n) 00181 { 00182 it_assert1(n%4 == 2, "conference(int n); wrong size"); 00183 int pm=n-1; // p must be odd prime, not checked 00184 imat out(n,n); 00185 00186 out.set_submatrix(1,n-1,1,n-1, jacobsthal(pm)); 00187 out.set_submatrix(0,0,1,n-1, 1); 00188 out.set_submatrix(1,n-1,0,0, 1); 00189 out(0,0)=0; 00190 00191 return out; 00192 } 00193 00194 cmat toeplitz(const cvec &c, const cvec &r) { 00195 int size = c.size(); 00196 it_assert(size == r.size(), 00197 "toeplitz(): Incorrect sizes of input vectors."); 00198 cmat output(size, size); 00199 cvec c_conj = conj(c); 00200 c_conj[0] = conj(c_conj[0]); 00201 00202 for (int i = 0; i < size; i++) { 00203 cmat tmp = reshape(c_conj(0, size - 1 - i), size - i, 1); 00204 output.set_submatrix(i, size - 1, i, i, tmp); 00205 } 00206 for (int i = 0; i < size - 1; i++) { 00207 cmat tmp = reshape(r(1, size - 1 - i), 1, size - 1 - i); 00208 output.set_submatrix(i, i, i + 1, size - 1, tmp); 00209 } 00210 00211 return output; 00212 } 00213 00214 cmat toeplitz(const cvec &c) { 00215 int size = c.size(); 00216 cmat output(size, size); 00217 cvec c_conj = conj(c); 00218 c_conj[0] = conj(c_conj[0]); 00219 00220 for (int i = 0; i < size; i++) { 00221 cmat tmp = reshape(c_conj(0, size - 1 - i), size - i, 1); 00222 output.set_submatrix(i, size - 1, i, i, tmp); 00223 } 00224 for (int i = 0; i < size - 1; i++) { 00225 cmat tmp = reshape(c(1, size - 1 - i), 1, size - 1 - i); 00226 output.set_submatrix(i, i, i + 1, size - 1, tmp); 00227 } 00228 00229 return output; 00230 } 00231 00232 mat toeplitz(const vec &c, const vec &r) { 00233 00234 mat output(c.size(), r.size()); 00235 00236 for (int i=0; i<c.size(); i++) { 00237 for(int j=0; j<std::min(r.size(), c.size()-i); j++) 00238 output(i+j,j) = c(i); 00239 } 00240 00241 for (int j=1; j<r.size(); j++) { 00242 for(int i=0; i<std::min(c.size(), r.size()-j); i++) 00243 output(i,i+j) = r(j); 00244 } 00245 00246 return output; 00247 } 00248 00249 mat toeplitz(const vec &c) { 00250 mat output(c.size(), c.size()); 00251 00252 for (int i=0; i<c.size(); i++) { 00253 for(int j=0; j<c.size()-i; j++) 00254 output(i+j,j) = c(i); 00255 } 00256 00257 for (int j=1; j<c.size(); j++) { 00258 for(int i=0; i<c.size()-j; i++) 00259 output(i,i+j) = c(j); 00260 } 00261 00262 return output; 00263 } 00264 00265 mat rotation_matrix(int dim, int plane1, int plane2, double angle) 00266 { 00267 mat m; 00268 double c = std::cos(angle), s = std::sin(angle); 00269 00270 it_assert(plane1>=0 && plane2>=0 && 00271 plane1<dim && plane2<dim && plane1!=plane2, 00272 "Invalid arguments to rotation_matrix()"); 00273 00274 m.set_size(dim, dim, false); 00275 m = 0.0; 00276 for (int i=0; i<dim; i++) 00277 m(i,i) = 1.0; 00278 00279 m(plane1, plane1) = c; 00280 m(plane1, plane2) = -s; 00281 m(plane2, plane1) = s; 00282 m(plane2, plane2) = c; 00283 00284 return m; 00285 } 00286 00287 void house(const vec &x, vec &v, double &beta) 00288 { 00289 double sigma, mu; 00290 int n = x.size(); 00291 00292 v = x; 00293 if (n == 1) { 00294 v(0) = 1.0; 00295 beta = 0.0; 00296 return; 00297 } 00298 sigma = energy(x(1, n-1)); 00299 v(0) = 1.0; 00300 if (sigma == 0.0) 00301 beta = 0.0; 00302 else { 00303 mu = std::sqrt(sqr(x(0)) + sigma); 00304 if (x(0) <= 0.0) 00305 v(0) = x(0) - mu; 00306 else 00307 v(0) = -sigma / (x(0) + mu); 00308 beta = 2 * sqr(v(0)) / (sigma + sqr(v(0))); 00309 v /= v(0); 00310 } 00311 } 00312 00313 void givens(double a, double b, double &c, double &s) 00314 { 00315 double t; 00316 00317 if (b == 0) { 00318 c = 1.0; 00319 s = 0.0; 00320 } 00321 else { 00322 if (fabs(b) > fabs(a)) { 00323 t = -a/b; 00324 s = -1.0 / std::sqrt(1 + t*t); 00325 c = s * t; 00326 } 00327 else { 00328 t = -b/a; 00329 c = 1.0 / std::sqrt(1 + t*t); 00330 s = c * t; 00331 } 00332 } 00333 } 00334 00335 void givens(double a, double b, mat &m) 00336 { 00337 double t, c, s; 00338 00339 m.set_size(2,2); 00340 00341 if (b == 0) { 00342 m(0,0) = 1.0; 00343 m(1,1) = 1.0; 00344 m(1,0) = 0.0; 00345 m(0,1) = 0.0; 00346 } 00347 else { 00348 if (fabs(b) > fabs(a)) { 00349 t = -a/b; 00350 s = -1.0 / std::sqrt(1 + t*t); 00351 c = s * t; 00352 } 00353 else { 00354 t = -b/a; 00355 c = 1.0 / std::sqrt(1 + t*t); 00356 s = c * t; 00357 } 00358 m(0,0) = c; 00359 m(1,1) = c; 00360 m(0,1) = s; 00361 m(1,0) = -s; 00362 } 00363 } 00364 00365 mat givens(double a, double b) 00366 { 00367 mat m(2,2); 00368 givens(a, b, m); 00369 return m; 00370 } 00371 00372 00373 void givens_t(double a, double b, mat &m) 00374 { 00375 double t, c, s; 00376 00377 m.set_size(2,2); 00378 00379 if (b == 0) { 00380 m(0,0) = 1.0; 00381 m(1,1) = 1.0; 00382 m(1,0) = 0.0; 00383 m(0,1) = 0.0; 00384 } 00385 else { 00386 if (fabs(b) > fabs(a)) { 00387 t = -a/b; 00388 s = -1.0 / std::sqrt(1 + t*t); 00389 c = s * t; 00390 } 00391 else { 00392 t = -b/a; 00393 c = 1.0 / std::sqrt(1 + t*t); 00394 s = c * t; 00395 } 00396 m(0,0) = c; 00397 m(1,1) = c; 00398 m(0,1) = -s; 00399 m(1,0) = s; 00400 } 00401 } 00402 00403 mat givens_t(double a, double b) 00404 { 00405 mat m(2,2); 00406 givens_t(a, b, m); 00407 return m; 00408 } 00409 00411 template void eye(int, mat &); 00413 template void eye(int, bmat &); 00415 template void eye(int, imat &); 00417 template void eye(int, cmat &); 00418 00419 } // namespace itpp
Generated on Wed Apr 18 11:23:31 2007 for IT++ by Doxygen 1.5.2