IT++ Logo

poly.cpp

Go to the documentation of this file.
00001 
00030 #include <itpp/signal/poly.h>
00031 #include <itpp/base/converters.h>
00032 #include <itpp/base/algebra/eigen.h>
00033 #include <itpp/base/specmat.h>
00034 #include <itpp/base/matfunc.h>
00035 
00036 
00037 namespace itpp
00038 {
00039 
00040 void poly(const vec &r, vec &p)
00041 {
00042   int n = r.size();
00043 
00044   p.set_size(n + 1, false);
00045   p.zeros();
00046   p(0) = 1.0;
00047 
00048   for (int i = 0; i < n; i++)
00049     p.set_subvector(1, i + 1, p(1, i + 1) - r(i)*p(0, i));
00050 }
00051 
00052 void poly(const cvec &r, cvec &p)
00053 {
00054   int n = r.size();
00055 
00056   p.set_size(n + 1, false);
00057   p.zeros();
00058   p(0) = 1.0;
00059 
00060   for (int i = 0; i < n; i++)
00061     p.set_subvector(1, i + 1, p(1, i + 1) - r(i)*p(0, i));
00062 }
00063 
00064 
00065 
00066 void roots(const vec &p, cvec &r)
00067 {
00068   int n = p.size(), m, l;
00069   ivec f = find(p != 0.0);
00070   m = f.size();
00071   vec v = p;
00072   mat A;
00073 
00074   if (m > 0 && n > 1) {
00075     v = v(f(0), f(m - 1));
00076     l = v.size();
00077 
00078     if (l > 1) {
00079 
00080       A = diag(ones(l - 2), -1);
00081       A.set_row(0, -v(1, l - 1) / v(0));
00082       r = eig(A);
00083       cvec d;
00084       cmat V;
00085       eig(A, d , V);
00086 
00087       if (f(m - 1) < n)
00088         r = concat(r, zeros_c(n - f(m - 1) - 1));
00089     }
00090     else {
00091       r.set_size(n - f(m - 1) - 1, false);
00092       r.zeros();
00093     }
00094   }
00095   else
00096     r.set_size(0, false);
00097 }
00098 
00099 void roots(const cvec &p, cvec &r)
00100 {
00101   int n = p.size(), m, l;
00102   ivec f;
00103 
00104   // find all non-zero elements
00105   for (int i = 0; i < n; i++)
00106     if (p(i) != 0.0)
00107       f = concat(f, i);
00108 
00109 
00110   m = f.size();
00111   cvec v = p;
00112   cmat A;
00113 
00114   if (m > 0 && n > 1) {
00115     v = v(f(0), f(m - 1));
00116     l = v.size();
00117 
00118     if (l > 1) {
00119       A = diag(ones_c(l - 2), -1);
00120       A.set_row(0, -v(1, l - 1) / v(0));
00121       r = eig(A);
00122       if (f(m - 1) < n)
00123         r = concat(r, zeros_c(n - f(m - 1) - 1));
00124     }
00125     else {
00126       r.set_size(n - f(m - 1) - 1, false);
00127       r.zeros();
00128     }
00129   }
00130   else
00131     r.set_size(0, false);
00132 }
00133 
00134 
00135 vec polyval(const vec &p, const vec &x)
00136 {
00137   it_error_if(p.size() == 0, "polyval: size of polynomial is zero");
00138   it_error_if(x.size() == 0, "polyval: size of input value vector is zero");
00139 
00140   vec out(x.size());
00141 
00142   out = p(0);
00143 
00144   for (int i = 1; i < p.size(); i++)
00145     out = p(i) + elem_mult(x, out);
00146 
00147   return out;
00148 }
00149 
00150 cvec polyval(const vec &p, const cvec &x)
00151 {
00152   it_error_if(p.size() == 0, "polyval: size of polynomial is zero");
00153   it_error_if(x.size() == 0, "polyval: size of input value vector is zero");
00154 
00155   cvec out(x.size());
00156 
00157   out = p(0);
00158 
00159   for (int i = 1; i < p.size(); i++)
00160     out = std::complex<double>(p(i)) + elem_mult(x, out);
00161 
00162   return out;
00163 }
00164 
00165 cvec polyval(const cvec &p, const vec &x)
00166 {
00167   it_error_if(p.size() == 0, "polyval: size of polynomial is zero");
00168   it_error_if(x.size() == 0, "polyval: size of input value vector is zero");
00169 
00170   cvec out(x.size());
00171 
00172   out = p(0);
00173 
00174   for (int i = 1; i < p.size(); i++)
00175     out = std::complex<double>(p(i)) + elem_mult(to_cvec(x), out);
00176 
00177   return out;
00178 }
00179 
00180 cvec polyval(const cvec &p, const cvec &x)
00181 {
00182   it_error_if(p.size() == 0, "polyval: size of polynomial is zero");
00183   it_error_if(x.size() == 0, "polyval: size of input value vector is zero");
00184 
00185   cvec out(x.size());
00186 
00187   out = p(0);
00188 
00189   for (int i = 1; i < p.size(); i++)
00190     out = p(i) + elem_mult(x, out);
00191 
00192   return out;
00193 }
00194 
00195 
00196 
00197 } // namespace itpp
SourceForge Logo

Generated on Sun Jul 26 08:54:57 2009 for IT++ by Doxygen 1.5.9