elliptical_proto.cpp
上传用户:jtjnyq9001
上传日期:2014-11-21
资源大小:3974k
文件大小:7k
源码类别:

3G开发

开发平台:

Visual C++

  1. //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2. //
  3. //  File = elliptical_proto.cpp
  4. //
  5. //  Prototype Elliptical Filter
  6. //
  7. #include <math.h>
  8. #include "misdefs.h"
  9. #include "parmfile.h"
  10. #include "elliptical_proto.h"
  11. #include "ipow.h"
  12. #include "filter_types.h"
  13. extern ParmFile ParmInput;
  14. #ifdef _DEBUG
  15.   extern ofstream *DebugFile;
  16. #endif
  17. #define UPPER_SUMMATION_LIMIT 5
  18. //======================================================
  19. //  constructor
  20. EllipticalPrototype::EllipticalPrototype( int filt_order,
  21.                                           double passband_ripple,
  22.                                           double stopband_ripple,
  23.                                           double little_k)
  24.                   :AnalogFilterPrototype( filt_order )
  25. {
  26.   double x, u, u4;
  27.   double term, sum;
  28.   double vv, ww, mu;
  29.   double little_q;
  30.   double big_D;
  31.   int min_order;
  32.   int i, m, r;
  33.   double numer, denom;
  34.   double p_sub_zero;
  35.   double xx, yy;
  36.   double aa, bb, cc;
  37.   std::complex<double> work;
  38.   //-------------------------------------------------
  39.   //  check viability of parameter set
  40.   x = pow( (1.0-little_k*little_k), 0.25);
  41.   u = (1.0 -x)/(2.0*(1.0+x));
  42.   u4 = u*u*u*u;
  43.   little_q = u*(1.+(2*u4*(1.+(7.5*u4*(1.+(10.*u4))))));
  44.   #ifdef _DEBUG
  45.     *DebugFile << "In EllipticalPrototype:" << endl;
  46.     *DebugFile << "   little_q = " << little_q << endl;
  47.   #endif
  48.   big_D = (pow(10.0, stopband_ripple/10.0) - 1.0)/
  49.                     (pow(10.0, passband_ripple/10.0) - 1.0);
  50.   min_order = (int)ceil(log10(16.0*big_D)/
  51.                         log10(1.0/little_q));
  52.   #ifdef _DEBUG
  53.   *DebugFile << "In elliptical filter minimum order is " << min_order << endl;
  54.   #endif
  55.   if(filt_order < min_order){
  56.     *DebugFile << "Fatal Error -- " << endl;
  57.     *DebugFile << "Specified performance cannot be acheived "
  58.               << "with an elliptical filter of specified order."
  59.               << endl;
  60.     *DebugFile << "Filter order must be at least " << min_order
  61.               << endl;
  62.     exit(-1);
  63.     }
  64.   //---------------------------------------------------------
  65.   // if filter is viable, generate prototype
  66.   Num_Poles = Filt_Order;
  67.   Pole_Locs = new std::complex<double>[Num_Poles];
  68.   if(Filt_Order%2)
  69.     Num_Zeros = Filt_Order-1;
  70.   else
  71.     Num_Zeros = Filt_Order;
  72.   Zero_Locs = new std::complex<double>[Num_Zeros];
  73.   //----------------------------------
  74.   numer = pow(10.0, passband_ripple/20.0) + 1.0;
  75.   vv = log(numer/(numer -2.0))/(2.*filt_order);
  76.   #ifdef _DEBUG
  77.     *DebugFile << "   vv = " << vv << endl;
  78.   #endif
  79.   //---------------------------------------
  80.   sum = 0.0;
  81.   for( m=0; m < UPPER_SUMMATION_LIMIT; m++){
  82.     term = ipow(-1.0, m);
  83.     term *= ipow(little_q, m*(m+1));
  84.     term *= sinh((2.0*m + 1)*vv);
  85.     sum += term;
  86.     }
  87.   numer = sum * sqrt(sqrt(little_q));
  88.   sum = 0.0;
  89.   for( m=1; m < UPPER_SUMMATION_LIMIT; m++){
  90.     term = ipow(-1.0, m);
  91.     term *= ipow(little_q, m*m);
  92.     term *= cosh(2.0 * m * vv);
  93.     sum += term;
  94.     }
  95.   p_sub_zero = fabs(numer/(0.5 + sum));
  96.   #ifdef _DEBUG
  97.     *DebugFile << "   p_sub_zero = " << p_sub_zero << endl;
  98.   #endif
  99.   //----------------------------------------
  100.   ww = 1.0 + little_k * p_sub_zero * p_sub_zero;
  101.   ww = sqrt(ww*(1.0 + p_sub_zero * p_sub_zero / little_k));
  102.   #ifdef _DEBUG
  103.     *DebugFile << "   ww = " << ww << endl;
  104.   #endif
  105.   //-----------------------------------------
  106.   r = (filt_order-(filt_order%2))/2;
  107.   B1_Coef = new double[r];
  108.   B0_Coef = new double[r];
  109.   A2_Coef = new double[r];
  110.   A1_Coef = new double[r];
  111.   A0_Coef = new double[r];
  112.   //-------------------------------------------
  113.   H_Zero = 1.0;
  114.   for(i=1; i<=r; i++){
  115.     if(filt_order%2){
  116.       mu = i;
  117.       }
  118.     else{
  119.       mu = i - 0.5;
  120.       }
  121.     //--------------------------------
  122.     sum = 0.0;
  123.     for(m=0; m < UPPER_SUMMATION_LIMIT; m++){
  124.       term = ipow(-1.0, m);
  125.       term *= ipow(little_q, m*(m+1));
  126.       term *= sin( (2*m+1) * PI * mu / filt_order );
  127.       sum += term;
  128.       }
  129.     numer = 2.0 * sum *sqrt(sqrt(little_q));
  130.     //--------------------------------------
  131.     sum = 0.0;
  132.     for(m=1; m < UPPER_SUMMATION_LIMIT; m++){
  133.       term = ipow(-1.0,m);
  134.       term *= ipow(little_q, m*m);
  135.       term *= cos( TWO_PI * m * mu / filt_order );
  136.       sum += term;
  137.       }
  138.     xx = numer/(1.0 + 2.0 * sum);
  139.     //---------------------------------------
  140.     yy = 1.0 - little_k * xx * xx;
  141.     yy = sqrt(yy * (1.0-(xx*xx/little_k)));
  142.     //-------------------------------------
  143.     aa = 1.0/(xx*xx);
  144.     A2_Coef[Num_Biquad_Sects] = 1.0;
  145.     A1_Coef[Num_Biquad_Sects] = 0.0;
  146.     A0_Coef[Num_Biquad_Sects] = aa;
  147.     #ifdef _DEBUG
  148.       *DebugFile << "   i = " << i << ",  aa = " << aa << endl;
  149.     #endif
  150.     //--------------------------------------
  151.     denom = 1.0 + ipow(p_sub_zero * xx, 2);
  152.     bb = 2.0 * p_sub_zero * yy/denom;
  153.     B1_Coef[Num_Biquad_Sects] = bb;
  154.     #ifdef _DEBUG
  155.       *DebugFile << "   i = " << i << ",  bb = " << bb << endl;
  156.     #endif
  157.     //------------------------------------
  158.     denom = ipow(denom, 2);
  159.     numer = ipow(p_sub_zero * yy, 2) + ipow(xx*ww, 2);
  160.     cc = numer/denom;
  161.     B0_Coef[Num_Biquad_Sects] = cc;
  162.     #ifdef _DEBUG
  163.       *DebugFile << "   i = " << i << ",  cc = " << cc << endl;
  164.     #endif
  165.     Num_Biquad_Sects++;
  166.     //---------------------------------------
  167.     //
  168.     H_Zero *= (cc/aa);
  169.     //----------------------------------------
  170.     // compute pair of pole locations
  171.     work = std::sqrt<double>(std::complex<double>(bb*bb - 4.0*cc, 0) );
  172.     Pole_Locs[i-1] = (std::complex<double>(-bb,0)-work)/2.0;
  173.     Pole_Locs[filt_order-i] = (std::complex<double>(-bb,0)+work)/2.0;
  174.     #ifdef _DEBUG
  175.       *DebugFile << "Pole[ " << (i-1) << " ] = ( " << Pole_Locs[i-1].real()
  176.                 << ", " << Pole_Locs[i-1].imag() << ")" << endl;
  177.     #endif
  178.     //--------------------------------------
  179.     // compute pair of zero locations
  180.     //
  181.     Zero_Locs[i-1] = std::complex<double>(0.0, sqrt(aa) );
  182.     Zero_Locs[Num_Zeros-i] = -Zero_Locs[i-1];
  183.     #ifdef _DEBUG
  184.       *DebugFile << "Zero[ " << (i-1) << " ] = ( " << Zero_Locs[i-1].real()
  185.                 << ", " << Zero_Locs[i-1].imag() << ")" << endl;
  186.     #endif
  187.     } //end of loop over i
  188.   //--------------------------------------
  189.   //  Finish up Ho
  190.   if( filt_order%2 ){
  191.     H_Zero *= p_sub_zero;
  192.     #ifdef _DEBUG
  193.       *DebugFile << "In Elliptiocal prototype, p_sub_zero = " << p_sub_zero << endl;
  194.     #endif
  195.     Pole_Locs[(filt_order-1)/2] = std::complex<double>(-p_sub_zero, 0.0);
  196.     Real_Pole = p_sub_zero;
  197.     }
  198.   else{
  199.     H_Zero *= pow(10.0, passband_ripple/(-20.0));
  200.     }
  201.   #ifdef _DEBUG
  202.     *DebugFile << "   H_Zero = " << H_Zero << endl;
  203.   #endif
  204.   return;
  205. }
  206. //==========================================================
  207. EllipticalPrototype::~EllipticalPrototype()
  208. {
  209. }