[Kst] extragear/graphics/kst/plugins/discretizing_filters

[Nicolas Brisset]

SVN commit 483315 by brisset:

cleanup template stuff

 D             filter.cpp  
 M  +61 -0     filter.h  
 D             polynom.cpp  
 M  +130 -0    polynom.h  


--- trunk/extragear/graphics/kst/plugins/discretizing_filters/filter.h #483314:483315
@@ -39,5 +39,66 @@
   S* x;
 };
 
+template<class S> filter<S>::filter(polynom<S>& Ns, polynom<S>& Ds, double dT)
+:Nz(0), Dz(0)
+{
+  // the greatest degree between N & D
+  n=Ns.GetDegree()>Ds.GetDegree() ? Ns.GetDegree() : Ds.GetDegree();
+  // allocate memory for state vector
+  x = new S[n];
+  // reset it
+  Reset();
+  // as well as the output
+  out=0.0;
+  // computes "z" transfer function
+  polynom<S> potzm1odt(0); potzm1odt[0]=1.0;                    // [2(z-1)/dT]^i
+  polynom<S> dzm1odt(1); dzm1odt[1]=2.0/dT; dzm1odt[0]=-2.0/dT; // 2(z-1)/dT
+  for (int i=0; i<=n; i++) {
+    // (z+1)^(n-i) = 1
+    polynom<S> pozp1(0); pozp1[0]=1.0;
+    // z+1
+    polynom<S> zp1(1); zp1[1]=1.0; zp1[0]=1.0;
+    // (z+1)^(n-i)
+    for (int j=i+1; j<=n; j++) pozp1 = pozp1 * zp1;
+    polynom<S> dNz(0); dNz[0]=Ns[i];
+    dNz = dNz * pozp1 * potzm1odt;
+    Nz = Nz + dNz;
+    polynom<S> dDz(0); dDz[0]=Ds[i];
+    dDz = dDz * pozp1 * potzm1odt;
+    Dz = Dz + dDz;
+    potzm1odt = potzm1odt * dzm1odt;
+  }
+}
+//------------------------------------------------------------------------------
+template<class S> void filter<S>::Reset()
+// reset state vector
+{
+  for (int i=0; i<n; i++) x[i]=0.0;
+}
+//------------------------------------------------------------------------------
+template<class S> void filter<S>::ConnectTo(S& in)
+{
+  // memorizes input address
+  this->in = &in;
+}
+//------------------------------------------------------------------------------
+template<class S> void filter<S>::NextTimeStep()
+{
+  // intermediate output (with new input, but just before clock)
+  S iOut = (x[n-1] + *in * Nz[n])/ Dz[n];
+  // update state vector (x[0] .. x[n-1] (goes from output to input)
+  for (int i=n-1; i>=0; i--)
+    x[i] = (i>0?x[i-1]:0) + Nz[i]*(*in) - Dz[i]*iOut;
+  // output
+  out = (x[n-1] + Nz[n]*(*in)) / Dz[n];
+}
+//------------------------------------------------------------------------------
+template<class S> filter<S>::~filter()
+{
+  // deallocate memory used by state vector
+  delete[] x;
+}
+//------------------------------------------------------------------------------
 
+
 #endif // #ifndef FILTER_H
--- trunk/extragear/graphics/kst/plugins/discretizing_filters/polynom.h #483314:483315
@@ -33,4 +33,134 @@
   int n;                                // degree of polynom (table size is n+1)
 };
 
+// constructor
+template<class S> polynom<S>::polynom(int n)
+{
+  C = new S[n+1];
+  for (int i=0; i<=n; i++) C[i]=0.0;
+  NulC = 0.0;
+  this->n = n;
+};
+//------------------------------------------------------------------------------
+// copy constructor
+template<class S> polynom<S>::polynom(const polynom& P) {
+  n = P.n;
+  C = new S[n+1];
+  for (int i=0; i<=n; i++) C[i]=P.C[i];
+}
+//------------------------------------------------------------------------------
+// destructor
+template<class S> polynom<S>::~polynom()
+{
+  delete[] C;
+};
+//------------------------------------------------------------------------------
+// coefficient acces (read or write)
+template<class S> S& polynom<S>::operator[] (int i)
+{
+  return i<=n ? C[i] : NulC;
+};
+//------------------------------------------------------------------------------
+// value at x
+template<class S> S polynom<S>::operator() (S x)
+{
+  S r=0.0;
+  S xpi=1.0;  // x at the power i
+  for (int i=0; i<=n; i++) {
+    r+=C[i]*xpi;
+    xpi*=x;
+  }
+  return r;
+};
+//------------------------------------------------------------------------------
+// set a polynom equal to another one
+template<class S> polynom<S>& polynom<S>::operator=(const polynom& P)
+{
+  if (n!=P.n) { // reallocate memory if different size
+    delete[] C;
+    n = P.n;
+    C = new S[n+1];
+    this->n = n;
+  }
+  for (int i=0; i<=n; i++) C[i]=P.C[i];
+  return *this;
+};
+//------------------------------------------------------------------------------
+// polynom sum
+template<class S> polynom<S> polynom<S>::operator+(const polynom& P)
+{
+  polynom r(n>P.n?n:P.n);
+  for (int i=0; i<=  n; i++) r.C[i]+=  C[i];
+  for (int i=0; i<=P.n; i++) r.C[i]+=P.C[i];
+  return r;
+};
+//------------------------------------------------------------------------------
+// polynom difference
+template<class S> polynom<S> polynom<S>::operator-(const polynom& P)
+{
+  polynom r(n>P.n?n:P.n);
+  for (int i=0; i<=  n; i++) r.C[i]+=  C[i];
+  for (int i=0; i<=P.n; i++) r.C[i]-=P.C[i];
+  return r;
+};
+//------------------------------------------------------------------------------
+// polynom product
+template<class S> polynom<S> polynom<S>::operator*(const polynom<S>& P)
+{
+  polynom r(n+P.n);
+  for (int i=0; i<=r.n; i++) {
+    r.C[i]=0;
+    for (int j=(0>i-n?0:i-n); j<=(i<P.n?i:P.n); j++) r.C[i]+=C[i-j]*P.C[j];
+  }
+  return r;
+};
+//------------------------------------------------------------------------------
+// scalar multipication
+template<class S> polynom<S> polynom<S>::operator*(const S& x)
+{
+  polynom r(n);
+  for (int i=0; i<=n; i++) r.C[i]=C[i]*x;
+  return r;
+};
+//------------------------------------------------------------------------------
+// scalar division
+template<class S> polynom<S> polynom<S>::operator/(const S& x)
+{
+  polynom<S> r(n);
+  for (int i=0; i<=n; i++) r.C[i]=C[i]/x;
+  return r;
+};
+//------------------------------------------------------------------------------
+// derivative
+template<class S> polynom<S> polynom<S>::Derivative()
+{
+  polynom<S> r(n-1);
+  for (int i=0; i<=r.n; i++) r.C[i]=(i+1)*C[i+1];
+  return r;
+};
+//------------------------------------------------------------------------------
+// integrate
+template<class S> polynom<S> polynom<S>::Integrate()
+{
+  polynom<S> r(n+1);
+  r.C[0]=0;
+  for (int i=1; i<=r.n; i++) r.C[i]=C[i-1]/i;
+  return r;
+};
+//------------------------------------------------------------------------------
+template<class S> int polynom<S>::GetDegree()
+{
+  return n;
+}
+//------------------------------------------------------------------------------
+template<class S> void polynom<S>::SetDegree(int n)
+{
+  if (this->n != n) {
+    delete[] C;
+    C = new double[n+1];
+    this->n = n;
+    for (int i=0; i<=n; i++) C[i]=0.0;
+  }
+}
+
 #endif // ifndef POLYNOM_H