00001
00005 #include "bpm/bpm_dsp.h"
00006
00007 static complex_t bessel_pole[] = {
00008 { -1.00000000000e+00, 0.00000000000e+00}, { -1.10160133059e+00, 6.36009824757e-01},
00009 { -1.32267579991e+00, 0.00000000000e+00}, { -1.04740916101e+00, 9.99264436281e-01},
00010 { -1.37006783055e+00, 4.10249717494e-01}, { -9.95208764350e-01, 1.25710573945e+00},
00011 { -1.50231627145e+00, 0.00000000000e+00}, { -1.38087732586e+00, 7.17909587627e-01},
00012 { -9.57676548563e-01, 1.47112432073e+00}, { -1.57149040362e+00, 3.20896374221e-01},
00013 { -1.38185809760e+00, 9.71471890712e-01}, { -9.30656522947e-01, 1.66186326894e+00},
00014 { -1.68436817927e+00, 0.00000000000e+00}, { -1.61203876622e+00, 5.89244506931e-01},
00015 { -1.37890321680e+00, 1.19156677780e+00}, { -9.09867780623e-01, 1.83645135304e+00},
00016 { -1.75740840040e+00, 2.72867575103e-01}, { -1.63693941813e+00, 8.22795625139e-01},
00017 { -1.37384121764e+00, 1.38835657588e+00}, { -8.92869718847e-01, 1.99832584364e+00},
00018 { -1.85660050123e+00, 0.00000000000e+00}, { -1.80717053496e+00, 5.12383730575e-01},
00019 { -1.65239648458e+00, 1.03138956698e+00}, { -1.36758830979e+00, 1.56773371224e+00},
00020 { -8.78399276161e-01, 2.14980052431e+00}, { -1.92761969145e+00, 2.41623471082e-01},
00021 { -1.84219624443e+00, 7.27257597722e-01}, { -1.66181024140e+00, 1.22110021857e+00},
00022 { -1.36069227838e+00, 1.73350574267e+00}, { -8.65756901707e-01, 2.29260483098e+00} };
00023
00024
00025 void _add_splane_pole( filterrep_t *r, complex_t z ) {
00026 if ( c_real( z ) < 0.0 ) r->pole[r->npoles++] = z;
00027 return;
00028 }
00029
00030
00031
00032 filterrep_t* create_splane_representation( filter_t *f ) {
00033
00034 char msg[80];
00035 int p = 0, i = 0;
00036 double theta, rip, eps, y;
00037 filterrep_t *r;
00038
00039
00040 r = (filterrep_t*) calloc( 1, sizeof( filterrep_t) );
00041
00042 if ( ! r ) {
00043 bpm_error( "Cannot allocate memory for s-plane representation.", __FILE__, __LINE__ );
00044 return NULL;
00045 }
00046
00047
00048 r->npoles = 0;
00049
00050
00051 if ( f->options & BESSEL ) {
00052 p = ( f->order * f->order ) / 4;
00053 if ( f->order & 1 ) _add_splane_pole( r, bessel_pole[p++] );
00054 for ( i=0; i<f->order/2; i++ ) {
00055 _add_splane_pole( r, bessel_pole[p] );
00056 _add_splane_pole( r, c_conj(bessel_pole[p]) );
00057 p++;
00058 }
00059 }
00060
00061
00062 if ( f->options & ( BUTTERWORTH | CHEBYSHEV ) ) {
00063 for( i=0; i<2*f->order; i++ ) {
00064 theta = (f->order & 1 ) ? i*PI/f->order : (i+.5)*PI/f->order;
00065 _add_splane_pole( r, c_exp(complex(0.,theta)) );
00066 }
00067 }
00068
00069
00070 if ( f->options & CHEBYSHEV ) {
00071 if ( f->cheb_ripple >= 0.0 ) {
00072 bpm_error( "Chebyshev ripple is must be < 0 dB!", __FILE__, __LINE__ );
00073 return NULL;
00074 }
00075
00076 rip = pow( 10.0, -f->cheb_ripple / 10.0 );
00077 eps = sqrt( rip - 1.0 );
00078 y = asinh( 1.0 / eps ) / (double) f->order;
00079 if ( y <= 0. ) {
00080 sprintf( msg, "Chebyshev ripple coefficient is %f, must be > 0", y );
00081 bpm_error( msg, __FILE__, __LINE__ );
00082 return NULL;
00083 } else {
00084
00085 for( i=0; i<r->npoles; i++ ) {
00086 r->pole[i] = complex( c_real(r->pole[i])*sinh(y), c_imag(r->pole[i])*cosh(y) );
00087 }
00088 }
00089 }
00090
00091 return r;
00092 }