Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

chebyshev.h

Go to the documentation of this file.
00001 // 
00002 // author="Tony Kirke"
00003 // Copyright(c) 1993-1998 Tony Kirke
00004 /*
00005  * SPUC - Signal processing using C++ - A DSP library
00006  * 
00007  * This program is free software; you can redistribute it and/or modify
00008  * it under the terms of the GNU General Public License as published by
00009  * the Free Software Foundation; either version 2, or (at your option)
00010  * any later version.
00011  * 
00012  * This program is distributed in the hope that it will be useful,
00013  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015  * GNU General Public License for more details.
00016  * 
00017  * You should have received a copy of the GNU General Public License
00018  * along with this program; if not, write to the Free Software
00019  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
00020 */
00021 #include <spuc.h>
00022 #include <complex.h>
00023 #include <iir_1st.h>
00024 #include <iir_2nd.h>
00025 namespace SPUC {
00033 
00034 
00035 
00036 template <class Numeric> class chebyshev 
00037 {
00038         private:
00039         long order;
00040         long odd;
00041         long n2;
00042         complex<double>* roots;
00043         iir_2nd< Numeric>* iir;
00044         iir_1st< Numeric>* iir_1;
00045         double gain;
00046         double* coeff;
00047         double epi;
00048 
00049         public:
00053         chebyshev(double fcd, long ord=1, double ripple=3.0) {
00054                 gain = 1;
00055                 order = ord;
00056                 epi = pow(10.0,(ripple/10.)) - 1.0;
00057                 epi = pow(epi,(1./(1.0*ord)));
00058                 double wca = tan(0.5*PI*fcd);
00060             n2 = (order+1)/2;
00061                 odd = (order%2);
00062                 roots = new complex<double>[n2];
00063                 if (odd) iir_1 = new iir_1st< Numeric >;
00064                 get_roots(wca,order,n2);
00065                 bilinear(n2);
00066                 iir = new iir_2nd< Numeric >[order/2];
00067                 get_coeff(n2);
00068         }
00070         ~chebyshev() {
00071           delete [] roots;
00072           if (odd) delete iir_1;
00073           delete [] iir;
00074         }
00076         void reset() {
00077           for (int j=odd;j<n2;j++) { 
00078                 iir[j-odd].reset();
00079                 if (odd) iir_1->reset();
00080           }
00081         } 
00083         void print() {
00084                 int j;
00085                 for (j=odd;j<n2;j++) { 
00086                         iir[j-odd].print();
00087                 }
00088                 if (odd) iir_1->print();
00089         } 
00091         Numeric clock(Numeric in) {
00092                 Numeric tmp = in;
00093                 for (int i=odd;i<n2;i++) {
00094                         tmp = iir[i-odd].clock(tmp);
00095                 }
00096                 if (odd) tmp = iir_1->clock(tmp);
00097                 return(gain*tmp);
00098         }
00099         private:
00101         void get_roots(double wp, long n, long n2) {
00102        long l = 0;
00103        if (n%2 == 0) l = 1;                                                 
00104            double arg;
00105            double x = 1/epi;
00106            double asinh = log(x + sqrt(1.0+x*x));
00107            double v0 = asinh/(double(n));
00108            double sm = sinh(v0);
00109            double cm = cosh(v0);
00110            for (int j=0;j<n2;j++) {
00111                  arg = -0.5*PI*l/float(n);
00112                  roots[j] = wp*complex<double>(-sm*cos(arg),cm*sin(arg));
00113                  l += 2;
00114           }
00115         }
00117         void bilinear(long n2) {
00118       double td;
00119           int j;
00120           const double a = 1.;
00121       if (odd) {
00122                 // For s-domain zero (or pole) at infinity we
00123                 // get z-domain zero (or pole) at z = -1
00124                 double tmp = (roots[0].real()-a)/(a+roots[0].real());
00125                 iir_1->set_coeff(-tmp);
00126                 gain *= 0.5*(1+tmp);
00127           }
00128       for (j=odd;j<n2;j++) {                       
00129                 td = a*a - 2*a*roots[j].real() + magsq(roots[j]);
00130                 roots[j] = complex<double>((a*a - magsq(roots[j])),
00131                                         2.0*a*roots[j].imag())/td;
00132           }
00133         }
00135         void get_coeff(long n2) {
00136         for (int j=odd;j<n2;j++) { 
00137                   iir[j-odd].set_coeff(-2*roots[j].real(),magsq(roots[j]));
00138                   gain *= (magsq(roots[j]) - 2*roots[j].real() + 1.0)/(4*magsq(roots[j])); 
00139           }
00140         }
00141 };
00142 } // namespace SPUC 

Generated on Fri Sep 16 11:02:25 2005 for spuc by  doxygen 1.4.4