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
Generated on Sun Sep 14 18:52:28 2008 for IT++ by Doxygen 1.5.6