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

Generated on Sat Apr 19 10:43:56 2008 for IT++ by Doxygen 1.5.5