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     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
SourceForge Logo

Generated on Wed Apr 18 11:45:34 2007 for IT++ by Doxygen 1.5.2