IT++ Logo Newcom Logo

specmat.cpp

Go to the documentation of this file.
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     it_assert(size > 0, "hadamard(): size is not a power of 2"); 
00135     int logsize = ceil_i(log2(static_cast<double>(size)));
00136     it_assert(pow2i(logsize) == size, "hadamard(): size is not a power of 2");
00137 
00138     imat H(size, size); H(0,0) = 1; 
00139         
00140     for (int i = 0; i < logsize; ++i) {
00141       int pow2 = 1 << i;
00142       for (int k = 0; k < pow2; ++k) {
00143         for (int l = 0; l < pow2; ++l) {
00144           H(k, l) = H(k, l);
00145           H(k+pow2, l) = H(k, l);
00146           H(k, l+pow2) = H(k, l);
00147           H(k+pow2, l+pow2) = (-1) * H(k, l);
00148         }
00149       }
00150     }
00151     return H;
00152   }
00153 
00154   imat jacobsthal(int p)
00155   {
00156     int quadratic_residue;
00157     imat out(p,p);
00158     int i, j;
00159 
00160     out = -1; // start with all elements equal to "-1"
00161   
00162     // Generate a complete list of quadratic residues
00163     for (i=0; i<(p-1)/2; i++) {
00164       quadratic_residue=((i+1)*(i+1))%p;
00165       // set this element in all rows (col-row) = quadratic_residue
00166       for (j=0; j<p; j++) { 
00167         out(j, (j+quadratic_residue)%p)=1;
00168       }
00169     }
00170 
00171     // set diagonal elements to zero
00172     for (i=0; i<p; i++) {
00173       out(i,i)=0;
00174     }
00175     return out;
00176   }
00177 
00178   imat conference(int n)
00179   {
00180     it_assert1(n%4 == 2, "conference(int n); wrong size");
00181     int pm=n-1; // p must be odd prime, not checked
00182     imat out(n,n);
00183 
00184     out.set_submatrix(1,n-1,1,n-1, jacobsthal(pm));
00185     out.set_submatrix(0,0,1,n-1, 1);
00186     out.set_submatrix(1,n-1,0,0, 1);
00187     out(0,0)=0;
00188 
00189     return out;
00190   }
00191 
00192   cmat toeplitz(const cvec &c, const cvec &r) {
00193     int size = c.size();
00194     it_assert(size == r.size(), 
00195               "toeplitz(): Incorrect sizes of input vectors.");
00196     cmat output(size, size);
00197     cvec c_conj = conj(c);
00198     c_conj[0] = conj(c_conj[0]);
00199 
00200     for (int i = 0; i < size; i++) {
00201       cmat tmp = reshape(c_conj(0, size - 1 - i), size - i, 1);
00202       output.set_submatrix(i, size - 1, i, i, tmp);
00203     }
00204     for (int i = 0; i < size - 1; i++) {
00205       cmat tmp = reshape(r(1, size - 1 - i), 1, size - 1 - i);
00206       output.set_submatrix(i, i, i + 1, size - 1, tmp);
00207     }
00208 
00209     return output;
00210   }
00211 
00212   cmat toeplitz(const cvec &c) {
00213     int size = c.size();
00214     cmat output(size, size);
00215     cvec c_conj = conj(c);
00216     c_conj[0] = conj(c_conj[0]);
00217 
00218     for (int i = 0; i < size; i++) {
00219       cmat tmp = reshape(c_conj(0, size - 1 - i), size - i, 1);
00220       output.set_submatrix(i, size - 1, i, i, tmp);
00221     }
00222     for (int i = 0; i < size - 1; i++) {
00223       cmat tmp = reshape(c(1, size - 1 - i), 1, size - 1 - i);
00224       output.set_submatrix(i, i, i + 1, size - 1, tmp);
00225     }
00226 
00227     return output;
00228   }
00229 
00230   mat toeplitz(const vec &c, const vec &r) {
00231 
00232     mat output(c.size(), r.size());
00233 
00234     for (int i=0; i<c.size(); i++) {
00235       for(int j=0; j<std::min(r.size(), c.size()-i); j++)
00236         output(i+j,j) = c(i);
00237     }
00238 
00239     for (int j=1; j<r.size(); j++) {
00240       for(int i=0; i<std::min(c.size(), r.size()-j); i++)
00241         output(i,i+j) = r(j);
00242     }
00243 
00244     return output;
00245   }
00246 
00247   mat toeplitz(const vec &c) {
00248     mat output(c.size(), c.size());
00249 
00250     for (int i=0; i<c.size(); i++) {
00251       for(int j=0; j<c.size()-i; j++)
00252         output(i+j,j) = c(i);
00253     }
00254 
00255     for (int j=1; j<c.size(); j++) {
00256       for(int i=0; i<c.size()-j; i++)
00257         output(i,i+j) = c(j);
00258     }
00259 
00260     return output;
00261   }
00262 
00263   mat rotation_matrix(int dim, int plane1, int plane2, double angle)
00264   {
00265     mat m;
00266     double c = std::cos(angle), s = std::sin(angle);
00267   
00268     it_assert(plane1>=0 && plane2>=0 &&
00269               plane1<dim && plane2<dim && plane1!=plane2,
00270               "Invalid arguments to rotation_matrix()");
00271 
00272     m.set_size(dim, dim, false);
00273     m = 0.0;
00274     for (int i=0; i<dim; i++)
00275       m(i,i) = 1.0;
00276 
00277     m(plane1, plane1) = c;
00278     m(plane1, plane2) = -s;
00279     m(plane2, plane1) = s;
00280     m(plane2, plane2) = c;
00281 
00282     return m;
00283   }
00284 
00285   void house(const vec &x, vec &v, double &beta)
00286   {
00287     double sigma, mu;
00288     int n = x.size();
00289 
00290     v = x;
00291     if (n == 1) {
00292       v(0) = 1.0;
00293       beta = 0.0;
00294       return;
00295     }
00296     sigma = energy(x(1, n-1));
00297     v(0) = 1.0;
00298     if (sigma == 0.0)
00299       beta = 0.0;
00300     else {
00301       mu = std::sqrt(sqr(x(0)) + sigma);
00302       if (x(0) <= 0.0)
00303         v(0) = x(0) - mu;
00304       else
00305         v(0) = -sigma / (x(0) + mu);
00306       beta = 2 * sqr(v(0)) / (sigma + sqr(v(0)));
00307       v /= v(0);
00308     }
00309   }
00310 
00311   void givens(double a, double b, double &c, double &s)
00312   {
00313     double t;
00314     
00315     if (b == 0) {
00316       c = 1.0;
00317       s = 0.0;
00318     }
00319     else {
00320       if (fabs(b) > fabs(a)) {
00321         t = -a/b;
00322         s = -1.0 / std::sqrt(1 + t*t);
00323         c = s * t;
00324       }
00325       else {
00326         t = -b/a;
00327         c = 1.0 / std::sqrt(1 + t*t);
00328         s = c * t;
00329       }
00330     }
00331   }
00332 
00333   void givens(double a, double b, mat &m)
00334   {
00335     double t, c, s;
00336 
00337     m.set_size(2,2);
00338     
00339     if (b == 0) {
00340       m(0,0) = 1.0;
00341       m(1,1) = 1.0;
00342       m(1,0) = 0.0;
00343       m(0,1) = 0.0;
00344     }
00345     else {
00346       if (fabs(b) > fabs(a)) {
00347         t = -a/b;
00348         s = -1.0 / std::sqrt(1 + t*t);
00349         c = s * t;
00350       }
00351       else {
00352         t = -b/a;
00353         c = 1.0 / std::sqrt(1 + t*t);
00354         s = c * t;
00355       }
00356       m(0,0) = c;
00357       m(1,1) = c;
00358       m(0,1) = s;
00359       m(1,0) = -s;
00360     }
00361   }
00362 
00363   mat givens(double a, double b)
00364   {
00365     mat m(2,2);
00366     givens(a, b, m);
00367     return m;
00368   }
00369 
00370 
00371   void givens_t(double a, double b, mat &m)
00372   {
00373     double t, c, s;
00374 
00375     m.set_size(2,2);
00376     
00377     if (b == 0) {
00378       m(0,0) = 1.0;
00379       m(1,1) = 1.0;
00380       m(1,0) = 0.0;
00381       m(0,1) = 0.0;
00382     }
00383     else {
00384       if (fabs(b) > fabs(a)) {
00385         t = -a/b;
00386         s = -1.0 / std::sqrt(1 + t*t);
00387         c = s * t;
00388       }
00389       else {
00390         t = -b/a;
00391         c = 1.0 / std::sqrt(1 + t*t);
00392         s = c * t;
00393       }
00394       m(0,0) = c;
00395       m(1,1) = c;
00396       m(0,1) = -s;
00397       m(1,0) = s;
00398     }
00399   }
00400 
00401   mat givens_t(double a, double b)
00402   {
00403     mat m(2,2);
00404     givens_t(a, b, m);
00405     return m;
00406   }
00407 
00409   template void eye(int, mat &);
00411   template void eye(int, bmat &);
00413   template void eye(int, imat &);
00415   template void eye(int, cmat &);
00416 
00417 } // namespace itpp
SourceForge Logo

Generated on Sat Aug 25 23:40:25 2007 for IT++ by Doxygen 1.5.2