IT++ Logo

specmat.cpp

Go to the documentation of this file.
00001 
00031 #include <itpp/base/specmat.h>
00032 #include <itpp/base/math/elem_math.h>
00033 #include <itpp/base/math/log_exp.h>
00034 #include <itpp/base/matfunc.h>
00035 
00036 
00037 namespace itpp
00038 {
00039 
00040 ivec find(const bvec &invector)
00041 {
00042   it_assert(invector.size() > 0, "find(): vector cannot be empty");
00043   ivec temp(invector.size());
00044   int pos = 0;
00045   for (int i = 0;i < invector.size();i++) {
00046     if (invector(i) == bin(1)) {
00047       temp(pos) = i;
00048       pos++;
00049     }
00050   }
00051   temp.set_size(pos, true);
00052   return temp;
00053 }
00054 
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 
00106 
00107 vec impulse(int size)
00108 {
00109   vec t(size);
00110   t.clear();
00111   t[0] = 1.0;
00112   return t;
00113 }
00114 
00115 vec linspace(double from, double to, int points)
00116 {
00117   if (points < 2) {
00118     // This is the "Matlab definition" of linspace
00119     vec output(1);
00120     output(0) = to;
00121     return output;
00122   }
00123   else {
00124     vec output(points);
00125     double step = (to - from) / double(points - 1);
00126     int i;
00127     for (i = 0; i < points; i++)
00128       output(i) = from + i * step;
00129     return output;
00130   }
00131 }
00132 
00133 vec zigzag_space(double t0, double t1, int K)
00134 {
00135   it_assert(K > 0, "zigzag_space:() K must be positive");
00136   ivec N = "0 1";
00137 
00138   int n = 2;
00139   for (int k = 0; k < K; k++) {
00140     ivec Nn = 2 * N;
00141     for (int i = 1; i < length(Nn); i += 2)  {
00142       Nn = concat(Nn, i);
00143       n++;
00144     }
00145     N = Nn;
00146   }
00147 
00148   vec T0 = linspace(t0, t1, n);
00149   vec Tt = zeros(n);
00150   for (int i = 0; i < n; i++) {
00151     Tt(i) = T0(N(i));
00152   }
00153   return Tt;
00154 }
00155 
00156 // Construct a Hadamard-imat of size "size"
00157 imat hadamard(int size)
00158 {
00159   it_assert(size > 0, "hadamard(): size is not a power of 2");
00160   int logsize = ceil_i(::log2(static_cast<double>(size)));
00161   it_assert(pow2i(logsize) == size, "hadamard(): size is not a power of 2");
00162 
00163   imat H(size, size);
00164   H(0, 0) = 1;
00165 
00166   for (int i = 0; i < logsize; ++i) {
00167     int pow2 = 1 << i;
00168     for (int k = 0; k < pow2; ++k) {
00169       for (int l = 0; l < pow2; ++l) {
00170         H(k, l) = H(k, l);
00171         H(k + pow2, l) = H(k, l);
00172         H(k, l + pow2) = H(k, l);
00173         H(k + pow2, l + pow2) = (-1) * H(k, l);
00174       }
00175     }
00176   }
00177   return H;
00178 }
00179 
00180 imat jacobsthal(int p)
00181 {
00182   int quadratic_residue;
00183   imat out(p, p);
00184   int i, j;
00185 
00186   out = -1; // start with all elements equal to "-1"
00187 
00188   // Generate a complete list of quadratic residues
00189   for (i = 0; i < (p - 1) / 2; i++) {
00190     quadratic_residue = ((i + 1) * (i + 1)) % p;
00191     // set this element in all rows (col-row) = quadratic_residue
00192     for (j = 0; j < p; j++) {
00193       out(j, (j + quadratic_residue) % p) = 1;
00194     }
00195   }
00196 
00197   // set diagonal elements to zero
00198   for (i = 0; i < p; i++) {
00199     out(i, i) = 0;
00200   }
00201   return out;
00202 }
00203 
00204 imat conference(int n)
00205 {
00206   it_assert_debug(n % 4 == 2, "conference(int n); wrong size");
00207   int pm = n - 1; // p must be odd prime, not checked
00208   imat out(n, n);
00209 
00210   out.set_submatrix(1, n - 1, 1, n - 1, jacobsthal(pm));
00211   out.set_submatrix(0, 0, 1, n - 1, 1);
00212   out.set_submatrix(1, n - 1, 0, 0, 1);
00213   out(0, 0) = 0;
00214 
00215   return out;
00216 }
00217 
00218 
00219 template <>
00220 const cmat toeplitz(const cvec &c)
00221 {
00222   int s = c.size();
00223   cmat output(s, s);
00224   cvec c_conj = conj(c);
00225   for (int i = 1; i < s; ++i) {
00226     for (int j = 0; j < s - i; ++j) {
00227       output(i + j, j) = c_conj(i);
00228     }
00229   }
00230   // start from j = 0 here, because the main diagonal is not conjugated
00231   for (int j = 0; j < s; ++j) {
00232     for (int i = 0; i < s - j; ++i) {
00233       output(i, i + j) = c(j);
00234     }
00235   }
00236   return output;
00237 }
00238 
00239 
00240 mat rotation_matrix(int dim, int plane1, int plane2, double angle)
00241 {
00242   mat m;
00243   double c = std::cos(angle), s = std::sin(angle);
00244 
00245   it_assert(plane1 >= 0 && plane2 >= 0 &&
00246             plane1 < dim && plane2 < dim && plane1 != plane2,
00247             "Invalid arguments to rotation_matrix()");
00248 
00249   m.set_size(dim, dim, false);
00250   m = 0.0;
00251   for (int i = 0; i < dim; i++)
00252     m(i, i) = 1.0;
00253 
00254   m(plane1, plane1) = c;
00255   m(plane1, plane2) = -s;
00256   m(plane2, plane1) = s;
00257   m(plane2, plane2) = c;
00258 
00259   return m;
00260 }
00261 
00262 void house(const vec &x, vec &v, double &beta)
00263 {
00264   double sigma, mu;
00265   int n = x.size();
00266 
00267   v = x;
00268   if (n == 1) {
00269     v(0) = 1.0;
00270     beta = 0.0;
00271     return;
00272   }
00273   sigma = sum(sqr(x(1, n - 1)));
00274   v(0) = 1.0;
00275   if (sigma == 0.0)
00276     beta = 0.0;
00277   else {
00278     mu = std::sqrt(sqr(x(0)) + sigma);
00279     if (x(0) <= 0.0)
00280       v(0) = x(0) - mu;
00281     else
00282       v(0) = -sigma / (x(0) + mu);
00283     beta = 2 * sqr(v(0)) / (sigma + sqr(v(0)));
00284     v /= v(0);
00285   }
00286 }
00287 
00288 void givens(double a, double b, double &c, double &s)
00289 {
00290   double t;
00291 
00292   if (b == 0) {
00293     c = 1.0;
00294     s = 0.0;
00295   }
00296   else {
00297     if (fabs(b) > fabs(a)) {
00298       t = -a / b;
00299       s = -1.0 / std::sqrt(1 + t * t);
00300       c = s * t;
00301     }
00302     else {
00303       t = -b / a;
00304       c = 1.0 / std::sqrt(1 + t * t);
00305       s = c * t;
00306     }
00307   }
00308 }
00309 
00310 void givens(double a, double b, mat &m)
00311 {
00312   double t, c, s;
00313 
00314   m.set_size(2, 2);
00315 
00316   if (b == 0) {
00317     m(0, 0) = 1.0;
00318     m(1, 1) = 1.0;
00319     m(1, 0) = 0.0;
00320     m(0, 1) = 0.0;
00321   }
00322   else {
00323     if (fabs(b) > fabs(a)) {
00324       t = -a / b;
00325       s = -1.0 / std::sqrt(1 + t * t);
00326       c = s * t;
00327     }
00328     else {
00329       t = -b / a;
00330       c = 1.0 / std::sqrt(1 + t * t);
00331       s = c * t;
00332     }
00333     m(0, 0) = c;
00334     m(1, 1) = c;
00335     m(0, 1) = s;
00336     m(1, 0) = -s;
00337   }
00338 }
00339 
00340 mat givens(double a, double b)
00341 {
00342   mat m(2, 2);
00343   givens(a, b, m);
00344   return m;
00345 }
00346 
00347 
00348 void givens_t(double a, double b, mat &m)
00349 {
00350   double t, c, s;
00351 
00352   m.set_size(2, 2);
00353 
00354   if (b == 0) {
00355     m(0, 0) = 1.0;
00356     m(1, 1) = 1.0;
00357     m(1, 0) = 0.0;
00358     m(0, 1) = 0.0;
00359   }
00360   else {
00361     if (fabs(b) > fabs(a)) {
00362       t = -a / b;
00363       s = -1.0 / std::sqrt(1 + t * t);
00364       c = s * t;
00365     }
00366     else {
00367       t = -b / a;
00368       c = 1.0 / std::sqrt(1 + t * t);
00369       s = c * t;
00370     }
00371     m(0, 0) = c;
00372     m(1, 1) = c;
00373     m(0, 1) = -s;
00374     m(1, 0) = s;
00375   }
00376 }
00377 
00378 mat givens_t(double a, double b)
00379 {
00380   mat m(2, 2);
00381   givens_t(a, b, m);
00382   return m;
00383 }
00384 
00386 template void eye(int, mat &);
00388 template void eye(int, bmat &);
00390 template void eye(int, imat &);
00392 template void eye(int, cmat &);
00393 
00394 } // namespace itpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Wed Jan 20 23:03:04 2010 for IT++ by Doxygen 1.6.2