Refactoring qf_cauer and adding tee types

qucs-filter/qf_cauer.cpp

@@ -1,5 +1,5 @@
 /***************************************************************************
-                               qf_cauer.cpp
+                               cauer.cpp
                              ----------------
     begin                : Mon Jan 02 2006
     copyright            : (C) 2006 by Vincent Habchi, F5RCS
@@ -17,644 +17,184 @@
 
 // Elliptic (Cauer) filters, odd order
 
+#include "qf_cauer.h"
+#include "qf_elliptic_functions.h"
 #include <cmath>
-#include <stdlib.h>
 #include <iostream>
+#include <stdio.h>
+#include <stdlib.h>
 #include <string>
 
-#undef _QF_CAUER_DEBUG
+namespace qf {
-//#define _QF_CAUER_DEBUG 1
+cauer::cauer(qf_float amin, qf_float amax, qf_float fc, qf_float fs,
-
+             qf_float imp = 1, qf_float bw = 0, qft type = LOWPASS,
-#include "filter.h"
+             bool is_tee = false)
-#include "qf_poly.h"
+    : filter(CAUER, type, imp, fc, bw, is_tee), a_(NULL) {
-#include "qf_filter.h"
+  if (amin > amax || (amin > 3 || amax < 3) ||
-#include "qf_cauer.h"
+      ((fc > fs && type_ == LOWPASS) || (fc < fs && type_ == HIGHPASS)) ||
-
+      ((type_ == BANDPASS || type_ == BANDSTOP) &&
-qf_cauer::qf_cauer (unsigned n, qf_double_t r, qf_double_t t) :
+       std::abs(fs - (fc * fc) / fs) < bw)) {
-  qf_filter (CAUER, LOWPASS, 1, one_over_pi / 2, 0), rho (r) {
-  ord = n;
-  th = sin (t);
-  a = new qf_double_t[ord + 1];
-  xfer ();
-  values ();
-  synth (LOWPASS);
-}
-
-qf_cauer::qf_cauer (qf_double_t amin, qf_double_t amax, qf_double_t fc,
-		    qf_double_t fs, qf_double_t r = 1, qf_double_t b = 0,
-		    qft type = LOWPASS) :
-  qf_filter (CAUER, type, r, fc, b), a (NULL) {
-  if (amin > amax)
     return;
-
-  if (amin > 3 || amax < 3)
-    return;
-
-  if ((fc > fs && type == LOWPASS) || (fc < fs && type == HIGHPASS))
-    return;
-
-  if ((type == BANDPASS || type == BANDSTOP) &&
-      fabs (fs - (fc * fc) / fs) < bw)
-    return;
-
-  normalize (amin, amax, fs, type);
-  xfer ();
-  values ();
-  synth (type);
-}
-
-qf_cauer::~qf_cauer (void) {
-  if (a != NULL)
-    delete[] a;
-}
-
-static qf_double_t FMAX (qf_double_t x, qf_double_t y) {
-  return ((x > y) ? x : y);
-}
-
-// This is extracted from "Handbook of mathematical functions"
-// Edited by M. Abramowitz & I. A. Stegun
-// U.S. National Bureau of Standards, June '64 / Dec. '72
-
-// K by the arithmetic/geometric mean (AGM) method
-qf_double_t qf_cauer::K (qf_double_t k) {
-  qf_double_t a = 1, b = sqrt (1 - k * k);
-  while (fabs (a - b) > K_ERR) {
-    qf_double_t t = a;
-    a = 0.5 * (a + b);
-    b = sqrt (t * b);
-  }
-  return pi / (2 * a);
-}
-
-// sn (u, m) by descending Landen transforms
-// m = k^2
-qf_double_t qf_cauer::sn (qf_double_t u, qf_double_t m) {
-  if (m < SK_ERR) {
-    // Final approx.
-    return sin (u) - 0.25 * m * cos (u) * (u - 0.5 * sin (2 * u));
-  } else {
-    qf_double_t kp = sqrt (1 - m);
-    qf_double_t smu = (1 - kp) / (1 + kp);
-    qf_double_t v = u / (1 + smu);
-    // Recurse
-    qf_double_t sn1 = sn (v, smu * smu);
-    return (1 + smu) * sn1 / (1 + smu * sn1 * sn1);
-  }
-}
-
-// Computes elliptic jacobi functions
-// Adapted from: Numerical receipes in C, pp. 264 et seq.
-
-// Computes Carlson's elliptic integral of the first kind
-qf_double_t qf_cauer::ellip_RF (qf_double_t x, qf_double_t y, qf_double_t z) {
-  qf_double_t alamb, ave, delx, dely, delz, e2, e3;
-  qf_double_t sqrtx, sqrty, sqrtz, xt, yt, zt;
-
-  // Constants
-  const qf_double_t THIRD = 1.0 / 3.0;
-  const qf_double_t C1 = 1.0 / 24.0;
-  const qf_double_t C2 = 0.1;
-  const qf_double_t C3 = 3.0 / 44.0;
-  const qf_double_t C4 = 1.0 / 14.0;
-
-  xt = x;
-  yt = y;
-  zt = z;
-  do {
-    sqrtx = sqrt (xt);
-    sqrty = sqrt (yt);
-    sqrtz = sqrt (zt);
-    alamb = sqrtx * (sqrty + sqrtz) + sqrty * sqrtz;
-    xt = 0.25 * (xt + alamb);
-    yt = 0.25 * (yt + alamb);
-    zt = 0.25 * (zt + alamb);
-    ave = THIRD * (xt + yt + zt);
-    delx = (ave - xt) / ave;
-    dely = (ave - yt) / ave;
-    delz = (ave - zt) / ave;
-  }
-  while (FMAX (FMAX (fabs (delx), fabs (dely)), fabs (delz)) > K_ERR1);
-
-  e2 = delx * dely - delz * delz;
-  e3 = delx * dely * delz;
-  return (1 + (C1 * e2 - C2 - C3 * e3) * e2 + C4 * e3) / sqrt (ave);
-}
-
-// K(k) = RF(0, 1 - k^2, 1) -> complete elliptic intergral of the 1st kind
-qf_double_t qf_cauer::ellip_K (qf_double_t k) {
-  return ellip_RF (0, 1 - k * k, 1);
-}
-
-// K'(k) = K(sqrt(1 - k^2)), even for small k's
-qf_double_t qf_cauer::Kp (qf_double_t k) {
-  qf_double_t Kp;
-  qf_double_t f1 = 1, f2, w = 1;
-  qf_double_t kb = 1;
-
-  Kp = f2 = 2 * ln2 - log (k);	// K' = ln (4 / k')
-  while (kb > K_ERR2) {
-    kb *= k * k;
-    f1 *= (w / (w + 1));
-    f2 -= 2 / (w * (w + 1));
-    Kp += f1 * f2 * kb;
-    w += 2;
-  }
-  return Kp;
-}
-
-// Compute the Jacobian elliptic functions sn(u,k), cn(u,k) and dn(u,k).
-qf_double_t
-qf_cauer::ellip_sncndn (qf_double_t uu, qf_double_t emmc,
-			qf_double_t& sn, qf_double_t& cn, qf_double_t& dn) {
-  qf_double_t a, b, c, d, emc, u;
-  qf_double_t em[14], en[14];
-  int i, ii, l;
-  bool bo;
-
-  emc = emmc;
-  d = 1 - emc;
-  u = uu;
-
-  if (emc) {
-    bo = (emc < 0);
-    if (bo) {
-      emc /= -1 / d;
-      u *= (d = sqrt (d));
-    }
-    a = 1;
-    dn = 1;
-
-    for (i = 1; i < 14; i++) {
-      l = i;
-      em[i] = a;
-      en[i] = (emc = sqrt (emc));
-      c = 0.5 * (a + emc);
-      if (fabs (a - emc) <= SN_ACC * a)
-	break;
-      emc *= a;
-      a = c;
   }
 
-    u *= c;
+  normalize(amin, amax, fs);
-    sn = sin (u);
+  xfer();
-    cn = cos (u);
+  values();
-
+  synth();
-    if (sn) {
-      a = cn / sn;
-      c *= a;
-      for (ii = l; ii > 0; ii--) {
-	b = em[ii];
-	a *= c;
-	c *= dn;
-	dn = (en[ii] + a) / (b + a);
-	a = c / b;
-      }
-      a = 1 / sqrt (c * c + 1);
-      sn = (sn >= 0 ? a : -a);
-      cn = sn * c;
-    }
-
-    if (bo) {
-      a = dn;
-      dn = cn;
-      cn = a;
-      sn /= d;
-    }
-  } else {
-    cn = 1 / cosh (u);
-    dn = cn;
-    sn = tanh (u);
-  }
-  return sn;
 }
 
-qf_double_t qf_cauer::ellip_sn (qf_double_t x, qf_double_t k) {
+cauer::~cauer(void) {
-  qf_double_t sn, cn, dn;
+  if (a_ != NULL) {
-  return ellip_sncndn (x, 1 - k * k, sn, cn, dn);
+    delete[] a_;
-}
+  }
-
-// Arc sin in degrees
-static qf_double_t ASIND (qf_double_t ang) {
-  return (180 * asin (ang) / pi);
 }
 
 // Normalize the filter parameters to Z = 1 O and w = 1 rad/s
 // and computes order
-void qf_cauer::normalize (qf_double_t amin, qf_double_t amax,
+void cauer::normalize(qf_float amin, qf_float amax, qf_float fs) {
-			  qf_double_t fs, qft type) {
+  qf_float Amax  = pow(10, -amin / 10);
-  qf_double_t Amax = pow (10, -amin / 10);
+  qf_float Amin  = pow(10, -amax / 10);
-  qf_double_t Amin = pow (10, -amax / 10);
+  qf_float Aemax = 1 - Amin;
-  qf_double_t Aemax = 1 - Amin;
+  qf_float Aemin = 1 - Amax;
-  qf_double_t Aemin = 1 - Amax;
+  qf_float sAmin = -10 * log10(Aemax) + amin;
-  qf_double_t sAmin = -10 * log10 (Aemax) + amin;
+  qf_float sAmax = -10 * log10(Aemin) + amax;
-  qf_double_t sAmax = -10 * log10 (Aemin) + amax;
+  qf_float sdiff = sAmax - sAmin;
-  qf_double_t sdiff = sAmax - sAmin;
+  qf_float kA    = pow(10, -sdiff / 20);
+  qf_float KA;
 
-#ifdef _QF_CAUER_DEBUG
+  if (kA < 0.001) {
-  std::cout << "amin + aemin = " << sAmin << " dB\n";
+    KA = Kp(kA) / K(kA);
-  std::cout << "amax + aemax = " << sAmax << " dB\n";
+  } else {
-  std::cout << "D(a) = " << sdiff << " dB\n";
+    KA = K(sqrt(1 - kA * kA)) / K(kA);
-#endif
+  }
 
-  qf_double_t kA = pow (10, -sdiff / 20);
+  rho_ = sqrt(Aemin);
-  qf_double_t KA;
 
-  if (kA < 0.001)
+  switch (type_) {
-    KA = Kp (kA) / K (kA);
-  else
-    KA = K (sqrt (1 - kA * kA)) / K (kA);
-
-  rho = sqrt (Aemin);
-
-  switch (type) {
   case LOWPASS:
-    th = fc / fs;
+    th_ = fc_ / fs;
     break;
   case HIGHPASS:
-    th = fs / fc;
+    th_ = fs / fc_;
     break;
   case BANDPASS:
-    th = bw / fabs (fs - (fc * fc) / fs);
+    th_ = bw_ / std::abs(fs - (fc_ * fc_) / fs);
     break;
-  case BANDSTOP :
+  case BANDSTOP:
-    th = fabs (fs * bw / (fs * fs - fc * fc));
+    th_ = std::abs(fs * bw_ / (fs * fs - fc_ * fc_));
-    bw = fabs (fs - (fc * fc) / fs);
+    bw_ = std::abs(fs - (fc_ * fc_) / fs);
     break;
   }
 
   // Calculate order
-  qf_double_t Kth = K (th) / K (sqrt (1 - th * th));
+  qf_float Kth = K(th_) / K(sqrt(1 - th_ * th_));
-  ord = (unsigned) ceil (Kth * KA);
+  ord_         = (unsigned)ceil(Kth * KA);
-  if ((ord % 2) == 0)
+  if ((ord_ % 2) == 0) {
-    ord++;
+    ord_++;
+  }
 
-#ifdef _QF_CAUER_DEBUG
+  a_ = new qf_float[ord_ + 1];
-  std::cout << "K'/K = " << Kth << ", K1/K'1 = " << KA << '\n';
-  std::cout << "Order = " << ord << '\n';
-#endif
-
-  a = new qf_double_t[ord + 1];
 }
 
 // A Cauer (or elliptic) filter has a symetric D(O)
 // D(O) = F (O) / P (O) = K * O * Prod {(O^2 + a^2(i)) / (a^2(i) * O^2 + 1)}
 // So that it is Chebichev in the passband and in the stopband
-void qf_cauer::xfer (void) {
+void cauer::xfer(void) {
-  int m = (ord - 1) / 2;
+  int m       = (ord_ - 1) / 2;
-  qf_double_t Ws = a[ord] = sqrt (th);
+  qf_float Ws = a_[ord_] = sqrt(th_);
-  qf_double_t k = K (th);
+  qf_float k             = K(th_);
   int u;
 
-#ifdef _QF_CAUER_DEBUG
+  for (unsigned i = 1; i < ord_; i++) {
-  std::cerr << "Computing filter of order " << ord << " with ";
+    qf_float j = (qf_float)i / (qf_float)ord_;
-  std::cerr << "rho = " << rho << " and theta = " << ASIND (th) << "°\n";
+    a_[i]      = Ws * sn(j * k, th_);
-  std::cerr << "k = " << k << '\n';
-#endif
-
-  for (unsigned i = 1; i < ord; i++) {
-    qf_double_t j = (qf_double_t) i / (qf_double_t) ord;
-    a[i] = Ws * sn (j * k, th);
-#ifdef _QF_CAUER_DEBUG
-    std::cerr << "a(" << i << ") = " << a[i] << '\n';
-#endif
   }
 
-#ifdef _QF_CAUER_DEBUG
+  qf_float delta = 1;
-  std::cerr << "Norm. puls. (Ws) = " << Ws << '\n';
+  for (u = 1; u < m + 2; u++) {
-#endif
+    delta *= a_[2 * u - 1] * a_[2 * u - 1];
-
+  }
-  qf_double_t delta = 1;
-  for (u = 1; u < m + 2; u++)
-    delta *= a[2 * u - 1] * a[2 * u - 1];
   delta /= Ws;
-  qf_double_t c = delta * sqrt (1 / (rho * rho) - 1);
+  qf_float c = delta * sqrt(1 / (rho_ * rho_) - 1);
-
-#ifdef _QF_CAUER_DEBUG
-  std::cerr << "D = " << delta << '\n';
-  std::cerr << "c = " << c << '\n';
-#endif
 
   // Computes D
-  F = qf_poly (1, 0, 0, 1);	// F(X) = X
+  F_ = poly(1, 0, 0, 1); // F(X) = X
-  P = qf_poly (c, 0, 0, 0);	// P(X) = c
+  P_ = poly(c, 0, 0, 0); // P(X) = c
 
   for (u = 1; u < m + 1; u++) {
-    qf_poly MF (1, 0, a[2 * u] * a[2 * u], 2);
+    poly MF(1, 0, a_[2 * u] * a_[2 * u], 2);
-    qf_poly MP (a[2 * u] * a[2 * u], 0, 1, 2);
+    poly MP(a_[2 * u] * a_[2 * u], 0, 1, 2);
 
-    MF.disp ("MF");
+    MF.disp("MF");
-    MP.disp ("MP");
+    MP.disp("MP");
-    F *= MF;
+    F_ *= MF;
-    P *= MP;
+    P_ *= MP;
   }
 
-  F.disp ("F");
+  F_.disp("F");
-  P.disp ("P");
+  P_.disp("P");
 
   // E(x) * E(-x) = F(x) * F(-x) + P(x) * P(-x)
-  E = F.hsq () + P.hsq ();
+  E_ = F_.hsq() + P_.hsq();
 
-  E.disp ("E");
+  E_.disp("E");
-  E.hurw ();
+  E_.hurw();
-  E.disp ("E");
+  E_.disp("E");
 
-  BN = E.odd () + F.odd ();
+  BN_ = E_.odd() + F_.odd();
-  BD = E.even () - F.even ();
+  BD_ = E_.even() - F_.even();
 }
 
-void qf_cauer::values (void) {
+void cauer::values(void) {
-
+  n_comp_ = (3 * ord_) / 2;
-  ncomp = (3 * ord) / 2;
-  Comp = (qfc *) malloc (sizeof (qfc) * ncomp);
 
   // For each zero of transmission, we apply the method as in
   // Saal & Ulbrich p. 288 or Zverev pp. 129 et seq.
-  qf_double_t Ws = sqrt (th);
+  qf_float Ws = sqrt(th_);
-  for (unsigned k = 0, l = 2; k < (ncomp - 1); k += 3) {
+  for (unsigned k = 0, l = 2; k < (n_comp_ - 1); k += 3) {
-#ifdef _QF_CAUER_DEBUG
+    extract_pole_pCsLC(1 / a_[l], Ws);
-    std::cerr << "Pole (" << l << ") = " << (1 / (a[l] * Ws)) << "\n";
-#endif
-    extract_pole_pCsLC (1 / a[l], &Comp[k], Ws);
 
     // Zeros mangeling
-    l = ord - l + 1;
+    l = ord_ - l + 1;
-    if (l < (ord + 1) / 2)
+    if (l < (ord_ + 1) / 2) {
       l += 2;
     }
-
+  }
   // Final removal of inifite pole
-  Comp[ncomp - 1].val = Ws * BN.eval (1) / BD.eval (1);
+  subsection final_capa;
+  final_capa.wiring  = Wiring::SHUNT;
+  final_capa.content = Content::CAPA;
+  final_capa.capa_v  = Ws * BN_.eval (1) / BD_.eval (1);
+  proto_subsecs_.push_back(final_capa);
 }
 
-void qf_cauer::synth (qft type) {
+void cauer::synth() {
-  qf_double_t cnrm = 1 / (2 * pi * fc * imp);
+  switch (type_) {
-  qf_double_t lnrm = imp / (2 * pi * fc);
-  unsigned i, node;
-
-  switch (type) {
   case LOWPASS:
-    for (i = 0, node = 1;;) {
+    for (subsection& v : proto_subsecs_) {
-      Comp[i].comp = CAP;	// Parallel capa first
+      v.transform_lp(subsecs_, imp_, fc_);
-      Comp[i].val *= cnrm;	// and last !
-      Comp[i].node1 = node;
-      Comp[i].node2 = 0;
-      i++;
-      if (i == ncomp)
-	break;
-      Comp[i].comp = CAP;	// Resonator cap
-      Comp[i].val *= cnrm;
-      Comp[i].node1 = node;
-      Comp[i].node2 = node + 1;
-      i++;
-      Comp[i].comp = IND;
-      Comp[i].val *= lnrm;
-      Comp[i].node1 = node;
-      Comp[i].node2 = ++node;
-      i++;
     }
     break;
   case HIGHPASS:
-    // First comes a serial cap
+    for (subsection& v : proto_subsecs_) {
-    Comp[0].comp = CAP;
+      v.transform_hp(subsecs_, imp_, fc_);
-    Comp[0].node1 = 1;
-    Comp[0].node2 = 2;
-    Comp[0].val = cnrm / Comp[0].val;
-    for (i = 1, node = 2; i < ncomp;) {
-      // Then a serial resonant circuit
-      Comp[i].comp = CAP;
-      Comp[i].node1 = node;
-      Comp[i].node2 = node + 1;
-      Comp[i].val = cnrm / Comp[i].val;
-      i++;
-      Comp[i].comp = IND;
-      Comp[i].node1 = node + 1;
-      Comp[i].node2 = 0;
-      Comp[i].val = lnrm / Comp[i].val;
-      i++;
-      // Then another serial cap
-      Comp[i].comp = CAP;
-      Comp[i].node1 = node;
-      node += 2;
-      Comp[i].node2 = node;
-      Comp[i].val = cnrm / Comp[i].val;
-      i++;
     }
     break;
-  case BANDPASS: {
-    // We double the number of components
-    ncomp = 3 * ord - 1;
-    qfc * Comp2 = (qfc *) malloc (sizeof (qfc) * ncomp);
-    qf_double_t q = fc / bw;
-
-    for (unsigned i = 0, j = 0, node = 1;;) {
-      Comp2[j].comp = CAP;
-      Comp2[j].node1 = node;
-      Comp2[j].node2 = 0;
-      Comp2[j++].val = Comp[i].val * q * cnrm;
-
-      Comp2[j].comp = IND;
-      Comp2[j].node1 = node;
-      Comp2[j].node2 = 0;
-      Comp2[j++].val = lnrm / (q * Comp[i++].val);
-      if (j == ncomp)
-	break;
-
-      qf_double_t c = Comp[i].val;
-      qf_double_t l = Comp[i + 1].val;
-      qf_double_t iw2 = l * c;
-#ifdef _QF_CAUER_DEBUG
-      std::cout << "O(inf) = " << sqrt (1 / iw2) << '\n';
-#endif
-      qf_double_t b = sqrt (1 + 4 * q * q * iw2);
-
-#ifdef _QF_CAUER_DEBUG
-      std::cout << "b = " << b << '\n';
-#endif
-
-      Comp2[j + 1].comp = IND;
-      Comp2[j + 1].node1 = node;
-      Comp2[j + 1].node2 = node + 1;
-      Comp2[j + 1].val = (b - 1) / (q * c * 2 * b);
-
-      Comp2[j + 3].comp = IND;
-      Comp2[j + 3].node1 = node + 1;
-      Comp2[j + 3].node2 = node + 2;
-      Comp2[j + 3].val = (b + 1) / (q * c * 2 * b);
-
-      Comp2[j].comp = CAP;
-      Comp2[j].node1 = node;
-      Comp2[j].node2 = node + 1;
-      Comp2[j].val = 1 / Comp2[j + 3].val;
-
-      Comp2[j + 2].comp = CAP;
-      Comp2[j + 2].node1 = node + 1;
-      Comp2[j + 2].node2 = node + 2;
-      Comp2[j + 2].val = 1 / Comp2[j + 1].val;
-      Comp2[j].val *= cnrm;
-      Comp2[j + 2].val *= cnrm;
-      Comp2[j + 1].val *= lnrm;
-      Comp2[j + 3].val *= lnrm;
-      j += 4;
-      i += 2;
-      node += 2;
-    }
-
-    delete[] Comp;
-    Comp = Comp2;
-    break;
-  }
-  case BANDSTOP: {
-    ncomp = 3 * ord - 1;
-    qfc * Comp2 = (qfc *) malloc (sizeof (qfc) * ncomp);
-    qf_double_t q = fc / bw;
-
-    for (unsigned i = 0, j = 0, node = 1;;) {
-      Comp2[j].comp = CAP;
-      Comp2[j].node1 = node;
-      Comp2[j].node2 = node + 1;
-      Comp2[j++].val = cnrm * Comp[i].val / q;
-
-      Comp2[j].comp = IND;
-      Comp2[j].node1 = node + 1;
-      Comp2[j].node2 = 0;
-      Comp2[j++].val = lnrm * q / Comp[i++].val;
-      if (j == ncomp)
-	break;
-
-      qf_double_t c = Comp[i].val;
-      qf_double_t l = Comp[i + 1].val;
-      qf_double_t w2 = 1 / (l * c);
-#ifdef _QF_CAUER_DEBUG
-      std::cout << "O(inf) = " << sqrt (w2) << "; q = " << q << '\n';
-#endif
-      qf_double_t b = sqrt (1 + 4 * q * q * w2);
-
-#ifdef _QF_CAUER_DEBUG
-      std::cout << "b = " << b << '\n';
-#endif
-      Comp2[j + 1].comp = IND;
-      Comp2[j + 1].node1 = node;
-      Comp2[j + 1].node2 = node + 2;
-      Comp2[j + 1].val = l * (b - 1) / (2 * b * q);
-
-      Comp2[j + 3].comp = IND;
-      Comp2[j + 3].node1 = node + 2;
-      Comp2[j + 3].node2 = node + 3;
-      Comp2[j + 3].val = l * (b + 1) / (2 * b * q);
-
-      Comp2[j].comp = CAP;
-      Comp2[j].node1 = node;
-      Comp2[j].node2 = node + 2;
-      Comp2[j].val = 1 / Comp2[j + 3].val;
-
-      Comp2[j + 2].comp = CAP;
-      Comp2[j + 2].node1 = node + 2;
-      Comp2[j + 2].node2 = node + 3;
-      Comp2[j + 2].val = 1 / Comp2[j + 1].val;
-      Comp2[j].val *= cnrm;
-      Comp2[j + 2].val *= cnrm;
-      Comp2[j + 1].val *= lnrm;
-      Comp2[j + 3].val *= lnrm;
-      j += 4;
-      i += 2;
-      node += 3;
-    }
-
-    delete[] Comp;
-    Comp = Comp2;
-    break;
-  }
-  }
-}
-
-void qf_cauer::dump () {
-  if (Comp == NULL)
-    return;			// Not initialized
-
-  std::cout << "qf_cauer::dump ()\n";
-
-  switch (type) {
-  case LOWPASS:
-    std::cout << "Lowpass ";
-    break;
-  case HIGHPASS:
-    std::cout << "Highpass ";
-    break;
   case BANDPASS:
-    std::cout << "Bandpass ";
+    for (subsection& v : proto_subsecs_) {
+      v.transform_bp(subsecs_, imp_, fc_, bw_);
+    }
     break;
   case BANDSTOP:
-    std::cout << "Bandstop ";
+    for (subsection& v : proto_subsecs_) {
+      v.transform_bs(subsecs_, imp_, fc_, bw_);
+    }
     break;
   }
-  std::cout << "of order " << ord << ", theta = "
+  if (is_tee_) {
-	    << ASIND (th) << "\x00B0, rho = " << rho << '\n';
+    for (subsection& v : subsecs_) {
-  dump_cout ();
+      v.pi_tee_switch();
-}
-
-void CC (void) {
-  unsigned o;
-  qf_double_t t, r;
-
-  do {
-    std::cout << "Order : ";
-    std::cin >> o;
-    if (o == 0)
-      return;
-    std::cout << "Reflexion (%) : ";
-    std::cin >> r;
-    r /= 100.0;
-    std::cout << "Angle (\x00B0) : ";
-    std::cin >> t;
-    t = pi * t / 180.0;
-    qf_cauer F (o, r, t);
-    F.dump ();
-  }
-  while (true);
-}
-
-#if TEST
-int main (void) {
-  qf_double_t amin, amax, fc, fs, bw, r;
-
-  while (true) {
-    std::cout << "Enter cutoff (Hz) [0 to end]: ";
-    std::cin >> fc;
-    if (fc == 0)
-      return 0;
-    std::cout << "Enter ripple in passband (dB): ";
-    std::cin >> amin;
-    std::cout << "Enter stopband frequency (Hz): ";
-    std::cin >> fs;
-    std::cout << "Enter minimal attenuation in stopband (dB): ";
-    std::cin >> amax;
-    std::cout << "Enter bandwith (0 for high- or low-pass) (Hz): ";
-    std::cin >> bw;
-    std::cout << "Enter reference impedance (Ohm): ";
-    std::cin >> r;
-
-    if (bw == 0) {
-      if (fs > fc) {
-	qf_cauer F (amin, amax, fc, fs, r, 0, LOWPASS);
-	F.dump ();
-      }
-      else {
-	qf_cauer F (amin, amax, fc, fs, r, 0, HIGHPASS);
-	F.dump ();
-      }
-    }
-
-    else {
-      if (amin < amax) {
-	qf_cauer F (amin, amax, fc, fs, r, bw, BANDPASS);
-	F.dump ();
-      }
-      else {
-	qf_cauer F (amax, amin, fc, fs, r, bw, BANDSTOP);
-	F.dump ();
-      }
     }
   }
 }
-#endif /* TEST */
+} // namespace qf

qucs-filter/qf_cauer.h

@@ -1,5 +1,5 @@
 /***************************************************************************
-                               qf_cauer.h
+                               cauer.h
                              --------------
     begin                : Mon Jan 02 2006
     copyright            : (C) 2006 by Vincent Habchi, F5RCS
@@ -18,47 +18,30 @@
 #ifndef _QF_CAUER_H
 #define _QF_CAUER_H
 
-const qf_double_t SN_ACC = 1e-5;	// Accuracy of sn(x) is SN_ACC^2
+#include "qf_filter.h"
-const qf_double_t K_ERR1 = 1e-8;	// Accuracy of K(k)
-const qf_double_t K_ERR2 = 1e-20;	// Accuracy of K(k)
-const qf_double_t K_ERR3 = 1e-6;	// Accuracy of K(k)
-const qf_double_t K_ERR =  4e-16;                // Accuracy of K (k)
-const qf_double_t SK_ERR = sqrt (K_ERR);         // Accuracy of sn (u, k)
-
-class qf_cauer : public qf_filter
-{
- private:
 
+namespace qf {
+class cauer : public filter {
+private:
   // Standard parameters
-  qf_double_t rho;			// Reflection coeff.
+  qf_float rho_; // Reflection coeff.
-  qf_double_t th;			// Modular angle
+  qf_float th_;  // Modular angle
 
   // Zeros of transmission
-  qf_double_t * a;
+  qf_float* a_;
 
 public:
-  qf_cauer (unsigned, qf_double_t, qf_double_t);
+  cauer(qf_float, qf_float, qf_float, qf_float, qf_float, qf_float, qft,
-  qf_cauer (qf_double_t, qf_double_t, qf_double_t, qf_double_t, qf_double_t,
+        bool is_tee);
-	    qf_double_t, qft);
+  virtual ~cauer(void);
-  virtual ~qf_cauer (void);
-
-  // Elliptic functions
-  static qf_double_t K (qf_double_t);
-  static qf_double_t Kp (qf_double_t);
-  static qf_double_t sn (qf_double_t, qf_double_t);
-  static qf_double_t ellip_RF (qf_double_t, qf_double_t, qf_double_t);
-  static qf_double_t ellip_K (qf_double_t);
-  static qf_double_t ellip_sncndn (qf_double_t, qf_double_t,
-				   qf_double_t&, qf_double_t&, qf_double_t&);
-  static qf_double_t ellip_sn (qf_double_t, qf_double_t);
 
   // Computes standard form
-  void normalize (qf_double_t, qf_double_t, qf_double_t, qft);
+  void normalize(qf_float, qf_float, qf_float);
 
-  void xfer (void);		// Computes xfer fctn
+  void xfer(void);      // Computes xfer fctn
-  void values (void);		// Computes norm values
+  void values(void);    // Computes norm values
-  virtual void synth (qft);	// Standard -> Actual form
+  virtual void synth(); // Standard -> Actual form
-  void dump (void);		// Dumps to cout
 };
+} // namespace qf
 
 #endif // _QF_CAUER_H

qucs-filter/qf_elliptic_functions.h

@@ -0,0 +1,215 @@
+#ifndef ELLIPTIC_FUNCTIONS_H
+#define ELLIPTIC_FUNCTIONS_H
+#include "qf_poly.h"
+
+namespace qf {
+const qf_float SN_ACC = 1e-5;        // Accuracy of sn(x) is SN_ACC^2
+const qf_float K_ERR1 = 1e-8;        // Accuracy of K(k)
+const qf_float K_ERR2 = 1e-20;       // Accuracy of K(k)
+const qf_float K_ERR3 = 1e-6;        // Accuracy of K(k)
+const qf_float K_ERR  = 4e-16;       // Accuracy of K (k)clear
+const qf_float SK_ERR = sqrt(K_ERR); // Accuracy of sn (u, k)
+
+static qf_float FMAX(qf_float x, qf_float y) {
+  return ((x > y) ? x : y);
+}
+
+// This is extracted from "Handbook of mathematical functions"
+// Edited by M. Abramowitz & I. A. Stegun
+// U.S. National Bureau of Standards, June '64 / Dec. '72
+// Complete Elliptic integral of the first kind
+// K by the arithmetic/geometric mean (AGM) method
+qf_float K(qf_float k) {
+  qf_float a = 1, b = sqrt(1 - k * k);
+  while (std::abs(a - b) > K_ERR) {
+    qf_float t = a;
+    a          = 0.5 * (a + b);
+    b          = sqrt(t * b);
+  }
+  return pi / (2 * a);
+}
+
+double find_th(double Kth) {
+  double th_low  = 0.0,
+         th_high = 1.0; // Initial guesses for binary search bounds
+  double th_mid, ratio;
+
+  while (th_high - th_low > K_ERR) {
+    th_mid = 0.5 * (th_low + th_high);
+    ratio  = K(th_mid) / K(sqrt(1.0 - th_mid * th_mid));
+
+    if (ratio >
+        Kth) { // If the calculated ratio is too high, need smaller theta
+      th_high = th_mid;
+    } else {
+      th_low = th_mid; // If the calculated ratio is too low, need larger theta
+    }
+  }
+  return th_mid;
+}
+
+double inv_K(double K_result) {
+  double k_low = 0.0, k_high = 1.0; // Initial guesses for binary search bounds
+  double k_mid, K_mid;
+
+  while (k_high - k_low > K_ERR) {
+    k_mid = 0.5 * (k_low + k_high);
+    K_mid = K(k_mid);
+
+    if (K_mid >
+        K_result) { // We need a smaller modulus if K(k_mid) is too large
+      k_high = k_mid;
+    } else {
+      k_low = k_mid; // We need a larger modulus if K(k_mid) is too small
+    }
+  }
+  return k_mid;
+}
+
+// sn (u, m) by descending Landen transforms
+// m = k^2
+qf_float sn(qf_float u, qf_float m) {
+  if (m < SK_ERR) {
+    // Final approx.
+    return sin(u) - 0.25 * m * cos(u) * (u - 0.5 * sin(2 * u));
+  } else {
+    qf_float kp  = sqrt(1 - m);
+    qf_float smu = (1 - kp) / (1 + kp);
+    qf_float v   = u / (1 + smu);
+    // Recurse
+    qf_float sn1 = sn(v, smu * smu);
+    return (1 + smu) * sn1 / (1 + smu * sn1 * sn1);
+  }
+}
+
+// Computes elliptic jacobi functions
+// Adapted from: Numerical receipes in C, pp. 264 et seq.
+
+// Computes Carlson's elliptic integral of the first kind
+qf_float ellip_RF(qf_float x, qf_float y, qf_float z) {
+  qf_float alamb, ave, delx, dely, delz, e2, e3;
+  qf_float sqrtx, sqrty, sqrtz, xt, yt, zt;
+
+  // Constants
+  const qf_float THIRD = 1.0 / 3.0;
+  const qf_float C1    = 1.0 / 24.0;
+  const qf_float C2    = 0.1;
+  const qf_float C3    = 3.0 / 44.0;
+  const qf_float C4    = 1.0 / 14.0;
+
+  xt = x;
+  yt = y;
+  zt = z;
+  do {
+    sqrtx = sqrt(xt);
+    sqrty = sqrt(yt);
+    sqrtz = sqrt(zt);
+    alamb = sqrtx * (sqrty + sqrtz) + sqrty * sqrtz;
+    xt    = 0.25 * (xt + alamb);
+    yt    = 0.25 * (yt + alamb);
+    zt    = 0.25 * (zt + alamb);
+    ave   = THIRD * (xt + yt + zt);
+    delx  = (ave - xt) / ave;
+    dely  = (ave - yt) / ave;
+    delz  = (ave - zt) / ave;
+  } while (FMAX(FMAX(std::abs(delx), std::abs(dely)), std::abs(delz)) > K_ERR1);
+
+  e2 = delx * dely - delz * delz;
+  e3 = delx * dely * delz;
+  return (1 + (C1 * e2 - C2 - C3 * e3) * e2 + C4 * e3) / sqrt(ave);
+}
+
+// K(k) = RF(0, 1 - k^2, 1) -> complete elliptic intergral of the 1st kind
+qf_float ellip_K(qf_float k) {
+  return ellip_RF(0, 1 - k * k, 1);
+}
+
+// K'(k) = K(sqrt(1 - k^2)), even for small k's
+qf_float Kp(qf_float k) {
+  qf_float Kp;
+  qf_float f1 = 1, f2, w = 1;
+  qf_float kb = 1;
+
+  Kp = f2 = 2 * ln2 - log(k); // K' = ln (4 / k')
+  while (kb > K_ERR2) {
+    kb *= k * k;
+    f1 *= (w / (w + 1));
+    f2 -= 2 / (w * (w + 1));
+    Kp += f1 * f2 * kb;
+    w += 2;
+  }
+  return Kp;
+}
+
+// Compute the Jacobian elliptic functions sn(u,k), cn(u,k) and dn(u,k).
+qf_float ellip_sncndn(qf_float uu, qf_float emmc, qf_float& sn, qf_float& cn,
+                      qf_float& dn) {
+  qf_float a, b, c, d, emc, u;
+  qf_float em[14], en[14];
+  int i, ii, l;
+  bool bo;
+
+  emc = emmc;
+  d   = 1 - emc;
+  u   = uu;
+
+  if (emc) {
+    bo = (emc < 0);
+    if (bo) {
+      emc /= -1 / d;
+      u *= (d = sqrt(d));
+    }
+    a  = 1;
+    dn = 1;
+
+    for (i = 1; i < 14; i++) {
+      l     = i;
+      em[i] = a;
+      en[i] = (emc = sqrt(emc));
+      c     = 0.5 * (a + emc);
+      if (std::abs(a - emc) <= SN_ACC * a) {
+        break;
+      }
+      emc *= a;
+      a = c;
+    }
+
+    u *= c;
+    sn = sin(u);
+    cn = cos(u);
+
+    if (sn) {
+      a = cn / sn;
+      c *= a;
+      for (ii = l; ii > 0; ii--) {
+        b = em[ii];
+        a *= c;
+        c *= dn;
+        dn = (en[ii] + a) / (b + a);
+        a  = c / b;
+      }
+      a  = 1 / sqrt(c * c + 1);
+      sn = (sn >= 0 ? a : -a);
+      cn = sn * c;
+    }
+
+    if (bo) {
+      a  = dn;
+      dn = cn;
+      cn = a;
+      sn /= d;
+    }
+  } else {
+    cn = 1 / cosh(u);
+    dn = cn;
+    sn = tanh(u);
+  }
+  return sn;
+}
+
+qf_float ellip_sn(qf_float x, qf_float k) {
+  qf_float sn, cn, dn;
+  return ellip_sncndn(x, 1 - k * k, sn, cn, dn);
+}
+} // namespace qf
+#endif // ELLIPTIC_FUNCTIONS_H

qucs-filter/qf_filter.cpp

@@ -1,5 +1,5 @@
 /***************************************************************************
-                               qf_filter.cpp
+                               filter.cpp
                              -----------------
     begin                : Mon Jan 02 2006
     copyright            : (C) 2006 by Vincent Habchi, F5RCS
@@ -14,432 +14,326 @@
  *   (at your option) any later version.                                   *
  *                                                                         *
  ***************************************************************************/
-
-#ifdef HAVE_CONFIG_H
-# include <config.h>
-#endif
-
 #include <cmath>
 #include <iostream>
 #include <sstream>
-#include <string>
 #include <stdlib.h>
-#include <string.h>
 
-#include "qf_poly.h"
 #include "qf_filter.h"
 
-// Constructor of a filter of nth order
+namespace qf {
 
-// Just initialize a few things
+filter::filter(qfk kind, qft ttype, qf_float imp, qf_float fc = 1,
-qf_filter::qf_filter (int n) :
+               qf_float bw = 0, bool is_tee = false)
-  type (LOWPASS), kind (UNDEF), ord (n), fc (0), bw (0), imp (1),
+    : type_(ttype), kind_(kind), is_tee_(is_tee), fc_(fc), bw_(bw), imp_(imp),
-  ncomp (0), Comp (NULL) {
+      n_comp_(0) {}
-}
-
-// Default constructor
-qf_filter::qf_filter () :
-  type (LOWPASS), kind (UNDEF), ord (0), fc (0), bw (0), imp (1),
-  ncomp (0), Comp (NULL) {
-}
-
-qf_filter::qf_filter (int n, qfk k, qft t) :
-  type (t), kind (k), ord (n), fc (0), bw (0), imp (1),
-  ncomp (0), Comp (NULL) {
-}
-
-qf_filter::qf_filter (qfk k, qft t, qf_double_t r,
-		      qf_double_t f = 1, qf_double_t b = 0) :
-  type (t), kind (k), fc (f), bw (b), imp (r), ncomp (0), Comp (NULL) {
-  fstart = fc - bw / 2;
-  fstop  = fc + bw / 2;
-}
-
-qf_filter::qf_filter (qfk k, qft t) :
-  type (t), kind (k), ord (0), fc (0), bw (0), imp (1),
-  ncomp (0), Comp (NULL) {
-}
 
 // Destructor of a filter
-qf_filter::~qf_filter (void) {
+filter::~filter(void) {
-  if (Comp != NULL)
-    free (Comp);
 }
 
 // Extraction routines
 
 // Extract finite attenuation pole
 // Result is a parallel cap, and a serial resonator (L // C)
-void qf_filter::extract_pole_pCsLC (qf_double_t p, qfc * pComp,
+void filter::extract_pole_pCsLC(qf_float p, qf_float Ws) {
-				    qf_double_t Ws) {
+  qf_float comp[3];
-  BN.disp ("BN");
+  BN_.disp("BN");
-  BD.disp ("BD");
+  BD_.disp("BD");
 
-  qf_double_t pl = -p * p;
+  qf_float pl   = -p * p;
-  qf_double_t bdpl = BD.evalX2 (pl);
+  qf_float bdpl = BD_.evalX2(pl);
 
   // Partial removal of infinite pole (first // cap)
   // c = [B(s)/s] (s^2 = - O^2)
-  qf_double_t c = ((BN << 1).evalX2 (pl)) / bdpl;
+  qf_float c = ((BN_ << 1).evalX2(pl)) / bdpl;
-  pComp->val = c * Ws;
+  comp[0] = c * Ws;
-  qf_poly cS (c, 0, 0, 1);
+  poly cS(c, 0, 0, 1);
-  BN = BN - (BD * cS);		// B = B - cs
+  BN_ = BN_ - (BD_ * cS); // B = B - cs
-  BN.disp ("BN");
+  BN_.disp("BN");
-  BN.div (0, p);
+  BN_.div(0, p);
 
   // Full removal of finite pole
   // c1 = (s B(s) / (s^2 + O^2)) @ s^2 = - O^2
-  BN.disp ("BN");
+  BN_.disp("BN");
-  qf_double_t c1 = (BN >> 1).evalX2 (pl) / bdpl;
+  qf_float c1      = (BN_ >> 1).evalX2(pl) / bdpl;
-  (pComp + 1)->val = c1;
+  comp[1] = c1;
-  (pComp + 2)->val = -Ws / (c1 * pl);
+  comp[2] = -Ws / (c1 * pl);
-  (pComp + 1)->val *= Ws;
+  comp[1] *= Ws;
 
   // 1/B = 1/B - (s/c1) / (s^2 + O^2)
-  BD = BD - (BN >> 1) * (1 / c1);
+  BD_ = BD_ - (BN_ >> 1) * (1 / c1);
-  BD.div (0, p);
+  BD_.div(0, p);
 
-  BN.disp ("BN");
+  BN_.disp("BN");
-  BD.disp ("BD");
+  BD_.disp("BD");
+  subsection fst_subsec;
+  fst_subsec.wiring  = Wiring::SHUNT;
+  fst_subsec.content = Content::CAPA;
+  fst_subsec.capa_v  = comp[0];
+  proto_subsecs_.push_back(fst_subsec);
+  subsection snd_subsec;
+  snd_subsec.wiring  = Wiring::SERIES;
+  snd_subsec.content = Content::PARA_CAPA_INDUC;
+  snd_subsec.capa_v  = comp[1];
+  snd_subsec.indc_v  = comp[2];
+  proto_subsecs_.push_back(snd_subsec);
 }
 
 // User readable value string.
-std::string qf_filter::num2str (qf_double_t num) {
+QString filter::num2str(qf_float num) {
   char c       = 0;
-  qf_double_t cal = fabs(num);
+  qf_float cal = std::abs(num);
-  if(cal > 1e-20) {
+  if (cal > 1e-20) {
     cal = log10(cal) / 3.0;
-    if(cal < -0.2)  cal -= 0.98;
+    if (cal < -0.2) {
+      cal -= 0.98;
+    }
     int expo = int(cal);
 
-    if(expo >= -5) if(expo <= 4)
+    if (expo >= -5) {
-      switch(expo) {
+      if (expo <= 4) {
-        case -5: c = 'f'; break;
+        switch (expo) {
-        case -4: c = 'p'; break;
+        case -5:
-        case -3: c = 'n'; break;
+          c = 'f';
-        case -2: c = 'u'; break;
-        case -1: c = 'm'; break;
-        case  1: c = 'k'; break;
-        case  2: c = 'M'; break;
-        case  3: c = 'G'; break;
-        case  4: c = 'T'; break;
-      }
-
-    if(c) num /= pow(10.0, (qf_double_t)(3*expo));
-  }
-
-  std::stringstream str;
-  str << num;
-  if(c) str << c;
-  return str.str();
-}
-
-// Dumps a filter to std::cout
-void qf_filter::dump_cout (void) {
-  std::string unit;
-
-  for (unsigned i = 0; i < ncomp; i++) {
-    switch (Comp[i].comp) {
-    case CAP:
-      std::cout << "Cap. from node ";
-      unit = "F";
           break;
-    case IND:
+        case -4:
-      std::cout << "Ind. from node ";
+          c = 'p';
-      unit = "H";
           break;
-    case RES:
+        case -3:
-      std::cout << "Res. from node ";
+          c = 'n';
-      unit = "Ohm";
           break;
-    }
+        case -2:
-
+          c = 'u';
-    std::cout << Comp[i].node1 << " to ";
-
-    if (Comp[i].node2 == 0)
-      std::cout << "ground, value = ";
-    else
-      std::cout << "node " << Comp[i].node2 << ", value = ";
-
-    if (Comp[i].val > 1)
-      std::cout << Comp[i].val << " " << unit << '\n';
-    else if (Comp[i].val >= 1e-3)
-      std::cout << Comp[i].val * 1e3 << " m" << unit << '\n';
-    else if (Comp[i].val >= 1e-6)
-      std::cout << Comp[i].val * 1e6 << " u" << unit << '\n';
-    else if (Comp[i].val >= 1e-9)
-      std::cout << Comp[i].val * 1e9 << " n" << unit << '\n';
-    else if (Comp[i].val >= 1e-12)
-      std::cout << Comp[i].val * 1e12 << " p" << unit << '\n';
-    else if (Comp[i].val >= 1e-15)
-      std::cout << Comp[i].val * 1e15 << " f" << unit << '\n';
-    else
-      std::cout << "0 " << unit << '\n';
-  }
-}
-
-std::string qf_filter::to_spice (void) {
-  std::ostringstream res;
-  unsigned i = 0;
-
-  res << "* SPICE netlist\n";
-  switch(type) {
-    case LOWPASS:
-      res << "* low-pass filter " << fc << " Hz cutoff\n";
           break;
-    case HIGHPASS:
+        case -1:
-      res << "* high-pass filter " << fc << " Hz cutoff\n";
+          c = 'm';
           break;
-    case BANDPASS:
+        case 1:
-      res << "* band-pass filter "
+          c = 'k';
-	  << fstart << " Hz ... " << fstop << " Hz\n";
           break;
-    case BANDSTOP:
+        case 2:
-      res << "* band-reject filter "
+          c = 'M';
-	  << fstart << " Hz ... " << fstop << " Hz\n";
           break;
-  }
+        case 3:
-  res << "* PI-type, impedance matching " << imp << " Ohm\n";
+          c = 'G';
-
-  for (unsigned ic = 1, il = 1, ir = 1; i < ncomp; i++) {
-    switch (Comp[i].comp) {
-    case CAP:
-      res << "C" << ic << "\t";
-      res << Comp[i].node1 << "\t" << Comp[i].node2;
-      res << "\t" << Comp[i].val << '\n';
-      // Insert parallel resistor
-      res << "R" << ir << "\t";
-      res << Comp[i].node1 << "\t" << Comp[i].node2;
-      res << "\t10G\n";
-      ic++;
-      ir++;
           break;
-    case IND:
+        case 4:
-      res << "L" << il << "\t";
+          c = 'T';
-      res << Comp[i].node1 << "\t"
-	<< Comp[i].node1 * 100 + Comp[i].node2 << "\t";
-      res << Comp[i].val << "\n";
-      il++;
-      // Insert serial resistor
-      res << "R" << ir << "\t";
-      res << Comp[i].node1 * 100 + Comp[i].node2
-	<< "\t" << Comp[i].node2 << "\t";
-      res << "1u\n";
-      ir++;
-      break;
-    case RES:
-      res << "R" << ir << "\t";
-      ir++;
           break;
         }
       }
-
-  res << "RI\t10000\t1\t" << imp << "\n";
-  res << "RL\t";
-
-  if (Comp[i - 1].node2 == 0)
-    res << Comp[i - 1].node1;
-  else
-    res << Comp[i - 1].node2;
-
-  res << "\t0\t" << imp << "\n";
-  res << "VCC\t10000\t0\tDC\t1\tAC\t1\n";
-  res << ".END\n";
-  return res.str();
-}
-
-void qf_filter::dump_spice (void) {
-  std::cout << to_spice().c_str();
-}
-
-std::string qf_filter::to_qucs (void) {
-  std::ostringstream res;
-  std::ostringstream wir;
-  qf_double_t Value, Value2;
-
-  // create the Qucs schematic
-  res << "<Qucs Schematic " << PACKAGE_VERSION << ">\n";
-
-  int x2, x = 20 + 40;
-  res << "<Components>\n";
-  res << "<Pac P1 1 " << x << " 320 18 -26 0 1 \"1\" 1 \""
-      << imp << " Ohm\" 1 \"0 dBm\" 0 \"1 GHz\" 0>\n";
-  res << "<GND * 1 " << x << " 350 0 0 0 0>\n";
-
-  int yc = 320, yl = 240;
-  unsigned int i, j = 0, k = ord, repser = 0, gndser = 1;
-
-  // connect left source
-  x = 20 + 40;
-  x2 = x + 110;
-  if(type == HIGHPASS) x2 -= 30;
-  wir << "<" << x << " 240 " << x << " 290 \"\" 0 0 0>\n"
-      << "<" << x << " 240 " << x2 << " 240 \"\" 0 0 0>\n";
-
-  x = 170;
-  if (type == BANDPASS || type == BANDSTOP) k = ord * 2;
-  for(i = 0; i < k && ncomp > 0; i++) {
-    // find series components to ground
-    unsigned int a = j;
-    while (Comp[a].node2 != 0) {
-      if (a+1 < ncomp && Comp[a].node2 == Comp[a+1].node1 &&
-	  Comp[a].comp != Comp[a+1].comp)
-	a++;
-      else
-	break;
     }
-    // handle series components to ground
-    if (Comp[a].node2 == 0) {
-      // wire down to shunt components
-      wir << "<" << x << " 240 " << x << " 290 \"\" 0 0 0>\n";
-      wir << "<" << x-35 << " 240 " << x+35 << " 240 \"\" 0 0 0>\n";
-      unsigned int b, c;
-      for (c = 0, b = j; b <= a; b++, c++) {
-	// type of component
-	switch (Comp[b].comp) {
-	case CAP: res << "<C C"; break;
-	case IND: res << "<L L"; break;
-	case RES: res << "<R R"; break;
-	}
-	res << j+1 << " 1 " << x;
-	res << " " << yc+c*60 << " 17 -26 0 1";
-
-	// value of component
-	res << " \"" << num2str(Comp[b].val).c_str();
-	switch (Comp[b].comp) {
-	case CAP: res << "F"; break;
-	case IND: res << "H"; break;
-	case RES: res << "Ohm"; break;
-	}
-	res << "\" 1>\n";
-      }
-      // place ground symbol here
-      res << "<GND * 1 " << x << " " << yc+c*60-30 << " 0 0 0 0>\n";
-      j = a + 1;
-      if (gndser < c) gndser = c;
-      repser = 0;
-    }
-    // handle parallel components in series
-    else {
-      if (repser) x += 35;
-      // connect in series
-      if (i != 0)
-	wir << "<" << x-70 << " 240 " << x-30 << " 240 \"\" 0 0 0>\n";
-      if (i != k - 1)
-	wir << "<" << x+70 << " 240 " << x+30 << " 240 \"\" 0 0 0>\n";
-      unsigned int b = j, c = 0;
-      do {
-	// type of component
-	switch (Comp[b].comp) {
-	case CAP: res << "<C C"; break;
-	case IND: res << "<L L"; break;
-	case RES: res << "<R R"; break;
-	}
-	res << j+1 << " 1 " << x;
-	res << " " << yl-c*35 << " -26 ";
     if (c) {
-	  res << "-44";
+      num /= pow(10.0, (qf_float)(3 * expo));
-	  wir << "<" << x-30 << " " << yl 
-	      << " " << x-30 << " " << yl-c*35 << " \"\" 0 0 0>\n"
-	      << "<" << x+30 << " " << yl
-	      << " " << x+30 << " " << yl-c*35 << " \"\" 0 0 0>\n";
     }
-	else
-	  res << "10";
-	res << " 0 0";
-
-	// value of component
-	res << " \"" << num2str(Comp[b].val).c_str();
-	switch (Comp[b].comp) {
-	case CAP: res << "F"; break;
-	case IND: res << "H"; break;
-	case RES: res << "Ohm"; break;
   }
-	res << "\" 1>\n";
+  QString str = QString::number(num);
-	b++;
+  if (c) {
-	c++;
+    str.append(c);
   }
-      while (b < ncomp && Comp[b-1].node1 == Comp[b].node1 &&
+  return str;
-	     Comp[b-1].node2 == Comp[b].node2);
+}
-      j = b;
+
-      repser++;
+QString filter::to_qucs() {
+  QString compos = "";
+  QString wires  = "";
+  int space      = 100;
+  int x          = 0;
+  // Draw left power source
+
+  compos += QString("<Pac P1 1 %1 290 18 -26 0 1 \"1\" 1 \"%2 Ohm\" 1 \"0 "
+                    "dBm\" 0 \"1 GHz\" 0>\n")
+                .arg(x)
+                .arg(imp_);
+  compos += QString("<GND * 1 %1 320 0 0 0 0>\n").arg(x);
+  wires += QString("<0 200 %1 200 \"\" 0 0 0>\n").arg(space / 2);
+  wires += QString("<%1 200 %1 260 \"\" 0 0 0>\n").arg(x);
+  // Draw filter sections
+  for (auto& subsec : subsecs_) {
+    x += space;
+    if (subsec.wiring == SERIES) {
+      if (subsec.content == INDUC) {
+        wires +=
+            QString("<%1 200 %2 200 0 0 0>\n").arg(x - space / 2).arg(x - 30);
+        compos += QString("<L L1 1 %1 %2 -26 -50 0 0 \"%3H\" 1>\n")
+                      .arg(x)
+                      .arg(200)
+                      .arg(num2str(subsec.indc_v));
+        wires +=
+            QString("<%1 200 %2 200 0 0 0>\n").arg(x + 30).arg(x + space / 2);
+      } else if (subsec.content == CAPA) {
+        wires +=
+            QString("<%1 200 %2 200 0 0 0>\n").arg(x - space / 2).arg(x - 30);
+        compos += QString("<C C1 1 %1 %2 -26 -50 0 0 \"%3F\" 1>\n")
+                      .arg(x)
+                      .arg(200)
+                      .arg(num2str(subsec.capa_v));
+        wires +=
+            QString("<%1 200 %2 200 0 0 0>\n").arg(x + 30).arg(x + space / 2);
+      } else if (subsec.content == PARA_CAPA_INDUC) {
+        wires +=
+            QString("<%1 200 %2 200 0 0 0>\n").arg(x - space / 2).arg(x - 30);
+        compos += QString("<C C1 1 %1 %2 -26 -80 0 0 \"%3F\" 1>\n")
+                      .arg(x)
+                      .arg(200)
+                      .arg(num2str(subsec.capa_v));
+        wires += QString("<%1 %2 %1 %3 0 0 0>\n").arg(x - 30).arg(200).arg(170);
+        wires += QString("<%1 %2 %1 %3 0 0 0>\n").arg(x + 30).arg(200).arg(170);
+        compos += QString("<L L1 1 %1 %2 -26 -90 0 0 \"%3H\" 1>\n")
+                      .arg(x)
+                      .arg(170)
+                      .arg(num2str(subsec.indc_v));
+        wires +=
+            QString("<%1 200 %2 200 0 0 0>\n").arg(x + 30).arg(x + space / 2);
+      } else if (subsec.content == SERIES_CAPA_INDUC) {
+        wires +=
+            QString("<%1 200 %2 200 0 0 0>\n").arg(x - space / 2).arg(x - 30);
+        compos += QString("<C C1 1 %1 %2 -26 -50 0 0 \"%3F\" 1>\n")
+                      .arg(x)
+                      .arg(200)
+                      .arg(num2str(subsec.capa_v));
+        wires +=
+            QString("<%1 200 %2 200 0 0 0>\n").arg(x + 30).arg(x + space / 2);
+        x += space;
+        wires +=
+            QString("<%1 200 %2 200 0 0 0>\n").arg(x - space / 2).arg(x - 30);
+        compos += QString("<L L1 1 %1 %2 -26 -50 0 0 \"%3H\" 1>\n")
+                      .arg(x)
+                      .arg(200)
+                      .arg(num2str(subsec.indc_v));
+        wires +=
+            QString("<%1 200 %2 200 0 0 0>\n").arg(x + 30).arg(x + space / 2);
       }
-    if (j >= ncomp) break;
+    } else if (subsec.wiring == SHUNT) {
-    x += 70;
+      if (subsec.content == INDUC) {
+        wires += QString("<%1 200 %2 200 0 0 0>\n")
+                     .arg(x - space / 2)
+                     .arg(x + space / 2);
+        wires += QString("<%1 %2 %1 %3 0 0 0>\n").arg(x).arg(200).arg(260);
+        compos += QString("<L L1 1 %1 %2 17 -26 0 1 \"%3H\" 1>\n")
+                      .arg(x)
+                      .arg(290)
+                      .arg(num2str(subsec.indc_v));
+        compos += QString("<GND * 1 %1 320 0 0 0 0>\n").arg(x);
+      } else if (subsec.content == CAPA) {
+        wires += QString("<%1 200 %2 200 0 0 0>\n")
+                     .arg(x - space / 2)
+                     .arg(x + space / 2);
+        wires += QString("<%1 %2 %1 %3 0 0 0>\n").arg(x).arg(200).arg(260);
+        compos += QString("<C C1 1 %1 %2 17 -26 0 1 \"%3F\" 1>\n")
+                      .arg(x)
+                      .arg(290)
+                      .arg(num2str(subsec.capa_v));
+        compos += QString("<GND * 1 %1 320 0 0 0 0>\n").arg(x);
+      } else if (subsec.content == SERIES_CAPA_INDUC) {
+        wires += QString("<%1 200 %2 200 0 0 0>\n")
+                     .arg(x - space / 2)
+                     .arg(x + space / 2);
+        compos += QString("<L L1 1 %1 %2 17 -26 0 1 \"%3H\" 1>\n")
+                      .arg(x)
+                      .arg(230)
+                      .arg(num2str(subsec.indc_v));
+        compos += QString("<C C1 1 %1 %2 17 -26 0 1 \"%3F\" 1>\n")
+                      .arg(x)
+                      .arg(290)
+                      .arg(num2str(subsec.capa_v));
+        compos += QString("<GND * 1 %1 320 0 0 0 0>\n").arg(x);
+      } else if (subsec.content == PARA_CAPA_INDUC) {
+        wires += QString("<%1 200 %2 200 0 0 0>\n")
+                     .arg(x - space / 2)
+                     .arg(x + space / 2);
+        wires += QString("<%1 %2 %1 %3 0 0 0>\n").arg(x).arg(200).arg(260);
+        compos += QString("<L L1 1 %1 %2 47 16 0 1 \"%3H\" 1>\n")
+                      .arg(x)
+                      .arg(290)
+                      .arg(num2str(subsec.indc_v));
+        wires += QString("<%1 260 %2 260 0 0 0>\n").arg(x).arg(x + 30);
+        compos += QString("<C C1 1 %1 %2 17 -26 0 1 \"%3F\" 1>\n")
+                      .arg(x + 30)
+                      .arg(290)
+                      .arg(num2str(subsec.capa_v));
+        wires += QString("<%1 320 %2 320 0 0 0>\n").arg(x).arg(x + 30);
+        compos += QString("<GND * 1 %1 320 0 0 0 0>\n").arg(x);
       }
-
-  if (ord & 1)
-    x += 110;
-  else
-    x += 70;
-  res << "<Pac P2 1 "<< x << " 320 18 -26 0 1 \"2\" 1 \""
-      << imp << " Ohm\" 1 \"0 dBm\" 0 \"1 GHz\" 0>\n"
-      << "<GND * 1 " << x << " 350 0 0 0 0>\n";
-
-  yc += 20 + gndser * 60;
-  Value = fc / 10.0;
-  Value2 = 10.0 * fc;
-  res << "<.SP SP1 1 70 " << yc << " 0 50 0 0 \"log\" 1 \""
-      << num2str(Value).c_str() << "Hz\" 1 \"" << num2str(Value2).c_str()
-      << "Hz\" 1 \"200\" 1 \"no\" 0 \"1\" 0 \"2\" 0>\n"
-      << "<Eqn Eqn1 1 260 " << yc+10
-      << " -28 15 0 0 \"dBS21=dB(S[2,1])\" 1 "
-      << "\"dBS11=dB(S[1,1])\" 1 \"yes\" 0>\n"
-      << "</Components>\n";
-
-  res << "<Wires>\n";
-
-  // internal wires
-  res << wir.str().c_str();
-
-  // connect right source
-  x2 = x - 110;
-  if(type == HIGHPASS) x2 += 30;
-  res << "<" << x << " 240 " << x << " 290 \"\" 0 0 0>\n"
-      << "<" << x2 << " 240 " << x << " 240 \"\" 0 0 0>\n";
-
-  res << "</Wires>\n"
-      << "<Diagrams>\n"
-      << "</Diagrams>\n"
-      << "<Paintings>\n";
-
-  res << "<Text 400 " << (yc+10) << " 12 #000000 0 \"";
-
-  switch (kind) {
-  case CAUER: res << "Cauer "; break;
-  case BUTT: res << "Butterworth "; break;
-  case CHEB: res << "Chebichev "; break;
-  case ICHEB: res << "Inverse Chebichev "; break;
-  case BESS: res << "Bessel "; break;
-  case UNDEF: res << "Undefined "; break;
     }
+  }
+  // Draw right power source
+  x += space;
+  wires += QString("<%1 200 %2 200 0 0 0>\n").arg(x - space / 2).arg(x);
+  compos += QString("<Pac P2 1 %1 290 18 -26 0 1 \"2\" 1 \"%2 Ohm\" 1 \"0 "
+                    "dBm\" 0 \"1 GHz\" 0>\n")
+                .arg(x)
+                .arg(imp_);
+  compos += QString("<GND * 1 %1 320 0 0 0 0>\n").arg(x);
+  wires += QString("<%1 200 %1 260 \"\" 0 0 0>\n").arg(x);
 
-  switch(type) {
+  QString s = "";
+  s += "<Qucs Schematic 24.2.1>\n";
+  s += "<Components>\n";
+  s += compos;
+
+  float Value = fc_ / 10.0;
+  float Value2 = 10.0 * fc_;
+  s += "<.SP SP1 1 70 " + QString::number(400)+ " 0 50 0 0 \"log\" 1 \"";
+  s += num2str(Value) + "Hz\" 1 \"" + num2str(Value2);
+  s += "Hz\" 1 \"200\" 1 \"no\" 0 \"1\" 0 \"2\" 0>\n";
+  s += "<Eqn Eqn1 1 260 " + QString::number(410);
+  s += " -28 15 0 0 \"dBS21=dB(S[2,1])\" 1 ";
+  s += "\"dBS11=dB(S[1,1])\" 1 \"yes\" 0>\n";
+  
+  s += "</Components>\n";
+  s += "<Wires>\n";
+  s += wires;
+  s += "</Wires>\n";
+  s += "<Diagrams>\n";
+  s += "</Diagrams>\n";
+  s += "<Paintings>\n";
+  s += "<Text 400 " + QString::number(400) + " 12 #000000 0 \"";
+
+  switch (kind_) {
+  case CAUER:
+    s += "Cauer ";
+    break;
+  case BUTT:
+    s += "Butterworth ";
+    break;
+  case CHEB:
+    s += "Chebichev ";
+    break;
+  case ICHEB:
+    s += "Inverse Chebichev ";
+    break;
+  case BESS:
+    s += "Bessel ";
+    break;
+  case UNDEF:
+    s += "Undefined ";
+    break;
+  }
+  switch (type_) {
   case LOWPASS:
-      res << "low-pass filter\\n"
+    s += "low-pass filter\\n" + num2str(fc_) + "Hz cutoff";
-	  << num2str(fc).c_str() << "Hz cutoff";
     break;
   case HIGHPASS:
-      res << "high-pass filter\\n"
+    s += "high-pass filter\\n" + num2str(fc_) + "Hz cutoff";
-	  << num2str(fc).c_str() << "Hz cutoff";
     break;
   case BANDPASS:
-      res << "band-pass filter\\n"
+    s += "band-pass filter\\n" + num2str(fc_ - bw_ / 2) + "Hz ... " +
-	  << num2str(fstart).c_str() << "Hz ... "
+         num2str(fc_ + bw_ / 2) + "Hz";
-	  << num2str(fstop).c_str() << "Hz";
     break;
   case BANDSTOP:
-      res << "band-reject filter\\n"
+    s += "band-reject filter\\n" + num2str(fc_ - bw_ / 2) + "Hz ... " +
-	  << num2str(fstart).c_str() << "Hz ... "
+         num2str(fc_ + bw_ / 2) + "Hz";
-	  << num2str(fstop).c_str() << "Hz";
     break;
   }
-  res << ", PI-type,\\nimpedance matching " << imp << " Ohm\">\n";
+  s += is_tee_ ? "\\n Tee-type" : "\\n PI-type";
+  s += "\\nimpedance matching " + num2str(imp_) + " Ohm\">\n";
+  s += "</Paintings>\n";
 
-  res << "</Paintings>\n";
+  // s = "wl-copy \'" + s + "\'";
-  return res.str();
+  // std::system(s.toStdString().c_str());
-}
+  // std::cout << s.toStdString();
-
+  return s;
-void qf_filter::dump_qucs (void) {
-  std::cout << to_qucs().c_str();
 }
+} // namespace qf

qucs-filter/qf_filter.h

@@ -1,5 +1,5 @@
 /***************************************************************************
-                               qf_filter.h
+                               filter.h
                              ----------------
     begin                : Mon Jan 02 2006
     copyright            : (C) 2006 by Vincent Habchi, F5RCS
@@ -19,100 +19,70 @@
 #define _QF_FILTER_H
 
 // Header for filter
+#include "qf_poly.h"
+#include "qf_subsection.h"
+#include <QString>
+#include <vector>
 
-enum qf_filter_type
+namespace qf {
-{
+enum filter_type { LOWPASS, HIGHPASS, BANDPASS, BANDSTOP };
-  LOWPASS,
-  HIGHPASS,
-  BANDPASS,
-  BANDSTOP
-};
 
-typedef enum qf_filter_type qft;
+typedef enum filter_type qft;
 
-enum qf_filter_kind
+enum filter_kind { BUTT, CHEB, ICHEB, BESS, CAUER, UNDEF };
-{
-  BUTT,
-  CHEB,
-  ICHEB,
-  BESS,
-  CAUER,
-  UNDEF
-};
 
-typedef enum qf_filter_kind qfk;
+typedef enum filter_kind qfk;
 
-enum qf_ctype
+enum ctype { CAP, IND, RES };
-{
-  CAP,
-  IND,
-  RES
-};
 
-typedef enum qf_ctype qfct;
+typedef enum ctype qfct;
 
-struct qf_comp
+struct comp {
-{
   qfct comp;
-  qf_double_t val;
+  qf_float val;
-  unsigned node1;
-  unsigned node2;
 };
 
-typedef struct qf_comp qfc;
+typedef struct comp qfc;
 
 // Generic filter class
 
-class qf_filter
+class filter {
-{
+protected:
- protected:
+  const qft type_; // Lowpass, highpass...
-  const qft type;		// Lowpass, highpass...
+  const qfk kind_; // Butterworth, Chebichev...
-  const qfk kind;		// Butterworth, Chebichev...
+  unsigned ord_;   // Order of filter
-  unsigned ord;			// Order of filter
+  bool is_tee_;    // Tee or Pi
 
-  const qf_double_t fc;		// Cutoff / Center
+  const qf_float fc_;  // Cutoff for lp or hp / Center for bp and bs
-  qf_double_t bw;	        // Bandwidth
+  qf_float bw_;        // Bandwidth
-  const qf_double_t imp;	// Terminating impedance
+  const qf_float imp_; // Terminating impedance
-  qf_double_t fstart;           // Start frequency
+  std::vector<subsection> proto_subsecs_;
-  qf_double_t fstop;            // Stop frequency
+  std::vector<subsection> subsecs_;
 
   // Polynomial description
 
-  qf_poly E;			// H(w) = E(w)/P(w)
+  poly E_; // H(w) = E(w)/P(w)
-  qf_poly F;			// D(w) = F(w)/P(w)
+  poly F_; // D(w) = F(w)/P(w)
-  qf_poly P;
+  poly P_;
 
-  qf_poly BN;			// B(w) susceptance of filter
+  poly BN_; // B(w) susceptance of filter
-  qf_poly BD;			// B(w) susceptance of filter
+  poly BD_; // B(w) susceptance of filter
 
-  unsigned ncomp;		// Number of components
+  unsigned n_comp_; // Number of components
-  qfc * Comp;			// Table of components
 
 public:
-  qf_filter (void);		// Default init
+  filter(qfk, qft, qf_float, qf_float, qf_float, bool);
-  qf_filter (int);		// Init
+  virtual ~filter(); // Exit
-  qf_filter (qfk, qft);
-  qf_filter (int, qfk, qft);
-  qf_filter (qfk, qft, qf_double_t, qf_double_t, qf_double_t);
-  virtual ~qf_filter (void);	// Exit
 
   // Common routines to perform extraction of poles and zeros
-
   // This one extracts a finite pole of transmission
-  void extract_pole_pCsLC (qf_double_t, qfc *, qf_double_t);
+  void extract_pole_pCsLC(qf_float, qf_float);
-
+  int order() { return ord_; }
-  int order (void) { return ord; }
+  virtual void synth() = 0; // Synthesize filter
-
+  QString to_qucs();
-  virtual void synth (qft) = 0;	// Synthesize filter
-
-  std::string to_qucs (void);	// Outputs Qucs
-  std::string to_spice (void);	// Outputs SPICE
-  void dump_qucs (void);        // Outputs Qucs to std::cout
-  void dump_spice (void);	// Outputs SPICE to std::cout
-  void dump_cout (void);	// Outputs to std::cout
 
 private:
-  std::string num2str (qf_double_t);
+  QString num2str(qf_float);
 };
-
+} // namespace qf
 #endif // _QF_FILTER_H

qucs-filter/qf_math.h

@@ -0,0 +1,10 @@
+#ifndef QF_MATH
+#define QF_MATH
+
+namespace qf {
+typedef double qf_float;
+static const double pi = 3.1415926535897932384626433832795029; /* pi       */
+static const double one_over_pi = 0.3183098861837906715377675267450287;                       /* 1/pi     */
+static const double ln2 = 0.6931471805599453094172321214581766; /* log_e(2) */
+} // namespace qf
+#endif

qucs-filter/qf_matrix.h

@@ -1,5 +1,5 @@
 /***************************************************************************
-                               qf_matrix.h
+                               matrix.h
                              ----------------
     begin                : Mon Jan 02 2006
     copyright            : (C) 2006 by Stefan Jahn
@@ -18,33 +18,33 @@
 #ifndef _QF_MATRIX_H
 #define _QF_MATRIX_H
 
-class qf_matrix
+#include "qf_math.h"
-{
+#include <stdio.h>
- public:
+#include <stdlib.h>
+
+namespace qf {
+
+class matrix {
+public:
   // constructor
-  qf_matrix (unsigned int d) {
+  matrix(unsigned int d) {
-    data = (qf_double_t *) calloc (d * d, sizeof (qf_double_t));
+    data = (qf_float*)calloc(d * d, sizeof(qf_float));
     n    = d;
   }
 
   // destructor
-  ~qf_matrix () {
+  ~matrix() { free(data); }
-    free (data);
-  }
 
   // accessor operators
-  qf_double_t operator () (int r, int c) const {
+  qf_float operator()(int r, int c) const { return data[r * n + c]; }
-    return data[r * n + c];
+  qf_float& operator()(int r, int c) { return data[r * n + c]; }
-  }
-  qf_double_t & operator () (int r, int c) {
-    return data[r * n + c];
-  }
 
   // size of matrix
   unsigned int n;
 
- private:
+private:
-  qf_double_t * data;
+  qf_float* data;
 };
+} // namespace qf
 
 #endif // _QF_MATRIX_H

qucs-filter/qf_poly.h

@@ -1,5 +1,5 @@
 /***************************************************************************
-                               qf_poly.h
+                               poly.h
                              ----------------
     begin                : Mon Jan 02 2006
     copyright            : (C) 2006 by Vincent Habchi, F5RCS
@@ -20,99 +20,98 @@
 
 /* Headers for R[X] arithmetic */
 
-#define qf_double_t long double
+#include "qf_math.h"
-
 #include "qf_matrix.h"
+#include <cmath>
 
+namespace qf {
 // A polynom can be described either by a product of monoms equal to
 // (x - r[i]) where r[i] is the ith root and a constant factor, or by
 // the classical series of coefficient a[0]...a[n], or both.
-
+enum poly_rep {
-enum qf_poly_rep
-{
   NONE  = 0, // Not initialized
   ROOTS = 1, // P(X) = k * prod (x - r[i])
   COEFF = 2, // P(X) = sum (a[i] * x^i)
   BOTH  = 3  // Both have been computed
 };
 
-typedef enum qf_poly_rep qpr;
+typedef enum poly_rep qpr;
 
-class qf_poly
+class poly {
-{
+private:
- private:
   qpr rep;       // Type of representation
   unsigned d;    // Current degree
-  qf_double_t krts;		// Constant k
+  qf_float krts; // Constant k
-  qf_double_t * p;		// Table of coefficients
+  qf_float* p;   // Table of coefficients
-  qf_double_t * rts;		// Table of complex roots
+  qf_float* rts; // Table of complex roots
 
   // Functions used by solve
-  void qf_bcm (qf_matrix &);
+  void bcm(matrix&);
-  int  qf_qrc (qf_matrix &, qf_double_t *);
+  int qrc(matrix&, qf_float*);
-  void qf_scm (qf_matrix &);
+  void scm(matrix&);
 
- public:
+public:
-  qf_poly ();
+  poly();
-  qf_poly (unsigned);		// id with d°
+  poly(unsigned);                               // id with d°
-  qf_poly (qf_double_t, qf_double_t, qf_double_t, unsigned); // Up to d°=2    
+  poly(qf_float, qf_float, qf_float, unsigned); // Up to d°=2
-  qf_poly (int, const qf_double_t[]);	// Id, with inst.
+  poly(int, const qf_float[]);                  // Id, with inst.
-  qf_poly (int, qf_double_t, const qf_double_t[]);
+  poly(int, qf_float, const qf_float[]);
-  qf_poly (const qf_poly &);	// Copy
+  poly(const poly&); // Copy
-  ~qf_poly ();
+  ~poly();
 
   // access operators
-  qf_poly & operator = (const qf_poly &); // P = Q
+  poly& operator=(const poly&); // P = Q
-  qf_double_t & operator [] (int i);	  // Access to element
+  qf_float& operator[](int i);  // Access to element
 
   // arithmetic operators
-  qf_poly operator - (void);	// Unary -
+  poly operator-(void); // Unary -
 
-  friend qf_poly operator + (qf_poly, qf_poly);
+  friend poly operator+(poly, poly);
-  friend qf_poly operator - (qf_poly, qf_poly);
+  friend poly operator-(poly, poly);
-  friend qf_poly operator * (qf_poly, qf_poly);
+  friend poly operator*(poly, poly);
-  friend qf_poly operator * (qf_poly, const qf_double_t);
+  friend poly operator*(poly, const qf_float);
 
-  qf_poly operator += (qf_poly);
+  poly operator+=(poly);
-  qf_poly operator -= (qf_poly);
+  poly operator-=(poly);
-  qf_poly operator *= (qf_poly);	// P(X) = P(X)*Q(X)
+  poly operator*=(poly); // P(X) = P(X)*Q(X)
-  qf_poly operator *= (const qf_double_t);
+  poly operator*=(const qf_float);
 
-  qf_poly operator << (unsigned);	// Basic div by X^n
+  poly operator<<(unsigned); // Basic div by X^n
-  qf_poly operator >> (unsigned);	// Multiply by X^n
+  poly operator>>(unsigned); // Multiply by X^n
 
-  bool operator == (qf_poly);	// Test
+  bool operator==(poly); // Test
-  bool operator != (qf_poly);	// Test
+  bool operator!=(poly); // Test
-  bool is_null (void);
+  bool is_null(void);
 
-  unsigned deg (void);		// Degree of poly
+  unsigned deg(void);        // Degree of poly
-  void spl (void);		// Simplify
+  void spl(void);            // Simplify
-  qf_poly odd (void);		// Odd part
+  poly odd(void);            // Odd part
-  qf_poly even (void);		// Even part
+  poly even(void);           // Even part
-  qf_poly mnx (void);		// P(X) -> P(-X)
+  poly mnx(void);            // P(X) -> P(-X)
-  qf_poly hsq (void);		// P(X)*P(-X)
+  poly hsq(void);            // P(X)*P(-X)
-  qf_poly sqr (void);		// Q(X) = P(X^2)
+  poly sqr(void);            // Q(X) = P(X^2)
-  qf_double_t eval (qf_double_t);	// P(X = a) 
+  qf_float eval(qf_float);   // P(X = a)
-  qf_double_t evalX2 (qf_double_t);	// P(X^2 = a) 
+  qf_float evalX2(qf_float); // P(X^2 = a)
 
-  void to_roots (void);		// Solves
+  void to_roots(void);          // Solves
-  qf_double_t k (void);		// Return krts factor
+  qf_float k(void);             // Return krts factor
-  void to_coeff (void);		// Calculate normal form
+  void to_coeff(void);          // Calculate normal form
-  void div (qf_double_t, qf_double_t);	// Simple division
+  void div(qf_float, qf_float); // Simple division
-  void hurw (void);		// "Hurwitzes" polynom
+  void hurw(void);              // "Hurwitzes" polynom
 
-  void disp (const char *);	// Prints P(X)
+  void disp(const char*); // Prints P(X)
-  void disp_c (void);
+  void disp_c(void);
-  void disp_r (void);
+  void disp_r(void);
 
-  friend void smpf (qf_poly &, qf_poly &);	// Simplify
+  friend void smpf(poly&, poly&); // Simplify
 };
 
 // For solve, we need some gibber
 
 // Save complex value elements
-#define SET_COMPLEX_PACKED(zp,n,r,i) \
+#define SET_COMPLEX_PACKED(zp, n, r, i)                                        \
-	*((zp)+2*(n)+0)=(r); *((zp)+2*(n)+1)=(i);
+  *((zp) + 2 * (n) + 0) = (r);                                                 \
+  *((zp) + 2 * (n) + 1) = (i);
 
 // Some constants
 
@@ -123,15 +122,16 @@ class qf_poly
 #define ROOT_PREC 1e-9
 #define ROOT_TOL 1e-7
 
-inline qf_double_t ROUND_ROOT (qf_double_t k) {
+inline qf_float ROUND_ROOT(qf_float k) {
-  if (k > 0)
+  if (k > 0) {
-    return floor (k / ROOT_PREC) * ROOT_PREC;
+    return std::floor(k / ROOT_PREC) * ROOT_PREC;
-  else
+  } else {
-    return ceil (k / ROOT_PREC) * ROOT_PREC;
+    return std::ceil(k / ROOT_PREC) * ROOT_PREC;
+  }
 }
 
 #define RADIX 2
-#define RADIX2		(RADIX*RADIX)
+#define RADIX2 (RADIX * RADIX)
 #define MAX_ITERATIONS 60
-
+} // namespace qf
 #endif // _QF_POLY_H

qucs-filter/qf_subsection.h

@@ -0,0 +1,174 @@
+#ifndef QF_SUBSECTION
+#define QF_SUBSECTION
+
+#include "qf_math.h"
+#include <iostream>
+#include <map>
+#include <string>
+#include <vector>
+
+namespace qf {
+
+enum Wiring { SHUNT, SERIES };
+
+enum Content { CAPA, INDUC, PARA_CAPA_INDUC, SERIES_CAPA_INDUC };
+
+struct subsection {
+  Wiring wiring;
+  Content content;
+  qf_float capa_v = 0;
+  qf_float indc_v = 0;
+  qf_float imp_   = 1;
+
+  void transform_lp(std::vector<subsection>& out_vec, qf_float imp,
+                    qf_float fc) {
+    qf_float cnrm = 1 / (2 * pi * fc * imp);
+    qf_float lnrm = imp / (2 * pi * fc);
+    subsection subsec;
+    subsec.wiring  = wiring;
+    subsec.content = content;
+    subsec.capa_v  = capa_v * cnrm;
+    subsec.indc_v  = indc_v * lnrm;
+    subsec.imp_    = imp;
+    out_vec.push_back(subsec);
+  }
+
+  void transform_hp(std::vector<subsection>& out_vec, qf_float imp,
+                    qf_float fc) {
+    subsection subsec;
+    qf_float lnrm = imp / (2 * pi * fc);
+    qf_float cnrm = 1 / (2 * pi * fc * imp);
+    subsec.indc_v = cnrm / capa_v;
+    subsec.capa_v = lnrm / indc_v;
+    subsec.wiring = wiring;
+    subsec.imp_   = imp;
+    if (content == CAPA) {
+      subsec.content = INDUC;
+    } else if (content == INDUC) {
+      subsec.content = CAPA;
+    } else {
+      subsec.content = content;
+    }
+
+    out_vec.push_back(subsec);
+  }
+
+  void transform_bp(std::vector<subsection>& out_vec, qf_float imp, qf_float fc,
+                    qf_float bw) {
+    qf_float q    = fc / bw;
+    qf_float cnrm = 1 / (2 * pi * fc * imp);
+    qf_float lnrm = imp / (2 * pi * fc);
+    if (content == CAPA) {
+      subsection subsec;
+      subsec.wiring  = wiring;
+      subsec.content = PARA_CAPA_INDUC;
+      subsec.indc_v  = lnrm / (q * capa_v);
+      subsec.capa_v  = cnrm * (q * capa_v);
+      subsec.imp_    = imp;
+      out_vec.push_back(subsec);
+    } else if (content == INDUC) {
+      subsection subsec;
+      subsec.wiring  = wiring;
+      subsec.content = SERIES_CAPA_INDUC;
+      subsec.capa_v  = lnrm / (q * indc_v);
+      subsec.indc_v  = cnrm / (indc_v * q);
+      subsec.imp_    = imp;
+      out_vec.push_back(subsec);
+    } else if (content == PARA_CAPA_INDUC) {
+      subsection fst_subsec;
+      subsection sec_subsec;
+      fst_subsec.wiring  = wiring;
+      sec_subsec.wiring  = wiring;
+      fst_subsec.imp_    = imp;
+      sec_subsec.imp_    = imp;
+      fst_subsec.content = PARA_CAPA_INDUC;
+      sec_subsec.content = PARA_CAPA_INDUC;
+      qf_float iw2       = indc_v * capa_v;
+      qf_float b         = sqrt(1 + 4 * q * q * iw2);
+      fst_subsec.capa_v  = cnrm * (q * capa_v * 2 * b) / (b + 1);
+      fst_subsec.indc_v  = lnrm * (b - 1) / (q * capa_v * 2 * b);
+      sec_subsec.capa_v  = cnrm * (q * capa_v * 2 * b) / (b - 1);
+      sec_subsec.indc_v  = lnrm * (b + 1) / (q * capa_v * 2 * b);
+      out_vec.push_back(fst_subsec);
+      out_vec.push_back(sec_subsec);
+    }
+  }
+
+  void transform_bs(std::vector<subsection>& out_vec, qf_float imp, qf_float fc,
+                    qf_float bw) {
+    qf_float q    = fc / bw;
+    qf_float cnrm = 1 / (2 * pi * fc * imp);
+    qf_float lnrm = imp / (2 * pi * fc);
+    if (content == CAPA) {
+      subsection subsec;
+      subsec.wiring  = wiring;
+      subsec.content = SERIES_CAPA_INDUC;
+      subsec.capa_v  = capa_v * cnrm / q;
+      subsec.indc_v  = q * lnrm / capa_v;
+      subsec.imp_    = imp;
+      out_vec.push_back(subsec);
+    } else if (content == INDUC) {
+      subsection subsec;
+      subsec.wiring  = wiring;
+      subsec.content = PARA_CAPA_INDUC;
+      subsec.indc_v  = indc_v * lnrm / q;
+      subsec.capa_v  = q * cnrm / indc_v;
+      subsec.imp_    = imp;
+      out_vec.push_back(subsec);
+    } else if (content == PARA_CAPA_INDUC) {
+      subsection fst_subsec;
+      subsection sec_subsec;
+      fst_subsec.wiring  = wiring;
+      sec_subsec.wiring  = wiring;
+      fst_subsec.imp_    = imp;
+      sec_subsec.imp_    = imp;
+      fst_subsec.content = PARA_CAPA_INDUC;
+      sec_subsec.content = PARA_CAPA_INDUC;
+      qf_float w2        = 1 / (indc_v * capa_v);
+      qf_float b         = sqrt(1 + 4 * q * q * w2);
+      fst_subsec.capa_v  = cnrm * (2 * b * q) / (indc_v * (b + 1));
+      fst_subsec.indc_v  = lnrm * indc_v * (b - 1) / (2 * b * q);
+      sec_subsec.capa_v  = cnrm * (2 * b * q) / (indc_v * (b - 1));
+      sec_subsec.indc_v  = lnrm * indc_v * (b + 1) / (2 * b * q);
+      out_vec.push_back(fst_subsec);
+      out_vec.push_back(sec_subsec);
+    }
+  }
+
+  void pi_tee_switch() {
+    wiring = wiring == Wiring::SHUNT ? Wiring::SERIES : Wiring::SHUNT;
+    std::map<Content, Content> content_map = {
+        {CAPA, INDUC},
+        {INDUC, CAPA},
+        {PARA_CAPA_INDUC, SERIES_CAPA_INDUC},
+        {SERIES_CAPA_INDUC, PARA_CAPA_INDUC}};
+    content = content_map[content];
+    std::swap(capa_v, indc_v);
+    capa_v /= imp_ * imp_;
+    indc_v *= imp_ * imp_;
+  }
+
+  void dump() {
+    std::string desc = "";
+    desc += wiring == Wiring::SHUNT ? "Shunt " : "Series ";
+    switch (content) {
+    case Content::CAPA:
+      desc += "Capa ";
+      break;
+    case Content::INDUC:
+      desc += "Induc ";
+      break;
+    case Content::SERIES_CAPA_INDUC:
+      desc += "Series Capa Induc ";
+      break;
+    case Content::PARA_CAPA_INDUC:
+      desc += "Shunt Capa Induc ";
+    default:
+      break;
+    }
+    std::cout << desc << capa_v << " " << indc_v << "\n";
+  }
+};
+} // namespace qf
+
+#endif

qucs-filter/qucsfilter.cpp

@@ -406,7 +406,7 @@ QString * QucsFilter::calculateFilter(struct tFilter * Filter)
         s = LC_Filter::createSchematic(Filter, false);
     }
     else  {
-      qf_cauer * F = NULL;
+      qf::cauer * F = NULL;
       double amin, amax, fc, fs, bw, r;
       fc = Filter->Frequency;
       amin = Filter->Ripple;
@@ -415,24 +415,25 @@ QString * QucsFilter::calculateFilter(struct tFilter * Filter)
       amax = Filter->Attenuation;
       bw = Filter->Frequency2 - fc;
 
+      bool is_tee = ComboRealize->currentIndex() == 1;
       switch (Filter->Class) {
       case CLASS_LOWPASS:
-        F = new qf_cauer (amin, amax, fc, fs, r, 0, LOWPASS);
+        F = new qf::cauer (amin, amax, fc, fs, r, 0, qf::LOWPASS, is_tee);
         break;
       case CLASS_HIGHPASS:
-        F = new qf_cauer (amin, amax, fc, fs, r, 0, HIGHPASS);
+        F = new qf::cauer (amin, amax, fc, fs, r, 0, qf::HIGHPASS, is_tee);
         break;  
       case CLASS_BANDPASS:
-        F = new qf_cauer (amin, amax, fc + bw / 2, fs, r, bw, BANDPASS);
+        F = new qf::cauer (amin, amax, fc + bw / 2, fs, r, bw, qf::BANDPASS, is_tee);
         break;
       case CLASS_BANDSTOP:
-        F = new qf_cauer (amin, amax, fc + bw / 2, fs, r, bw, BANDSTOP);
+        F = new qf::cauer (amin, amax, fc + bw / 2, fs, r, bw, qf::BANDSTOP, is_tee);
         break;
       }
       if (F) {
         //F->dump();
         EditOrder->setText(QString::number(F->order()));
-        s = new QString(F->to_qucs().c_str());
+        s = new QString(F->to_qucs());
         delete F;
       }
       else {