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