IT++ Logo

specmat.cpp

Go to the documentation of this file.
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
SourceForge Logo

Generated on Sat Apr 19 10:57:51 2008 for IT++ by Doxygen 1.5.5