00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
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
00123
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 }