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
Generated on Wed Dec 7 2011 03:38:55 for IT++ by Doxygen 1.7.4