From 4f1ee28d687d583601d41ff58e1618b381d2675f Mon Sep 17 00:00:00 2001 From: Katja Date: Sun, 6 Nov 2011 14:41:44 +0000 Subject: made creb compliant with double precision - changed float to t_float - adapted subnormal detection svn path=/trunk/externals/creb/; revision=15706 --- modules++/DSPIfilters.h | 206 ++++++++++++++++++++++++------------------------ 1 file changed, 103 insertions(+), 103 deletions(-) (limited to 'modules++/DSPIfilters.h') diff --git a/modules++/DSPIfilters.h b/modules++/DSPIfilters.h index 09268de..4fa53ea 100644 --- a/modules++/DSPIfilters.h +++ b/modules++/DSPIfilters.h @@ -33,9 +33,9 @@ class DSPIfilterOrtho { inline DSPIfilterOrtho(){resetState();resetCoef();resetSCoef();} inline ~DSPIfilterOrtho(){} - inline void resetState(){d1A = d1B = d2A = d2B = 0.0f;} - inline void resetCoef(){ai = ar = c0 = c1 = c2 = 0.0f;} - inline void resetSCoef(){s_ai = s_ar = s_c0 = s_c1 = s_c2 = 0.0f;} + inline void resetState(){d1A = d1B = d2A = d2B = 0.0;} + inline void resetCoef(){ai = ar = c0 = c1 = c2 = 0.0;} + inline void resetSCoef(){s_ai = s_ar = s_c0 = s_c1 = s_c2 = 0.0;} /* * Biquad filters remarks @@ -56,15 +56,15 @@ class DSPIfilterOrtho { */ // make sure freq and Q are positive and within bounds - inline void checkBounds(float &freq, float &Q) + inline void checkBounds(t_float &freq, t_float &Q) { freq = fabs(freq); Q = fabs(Q); - float epsilon = .0001f; // stability guard - float fmin = 0.0f + epsilon; - float fmax = 0.5f - epsilon; - float Qmin = 1.1f; + t_float epsilon = .0001; // stability guard + t_float fmin = 0.0 + epsilon; + t_float fmax = 0.5 - epsilon; + t_float Qmin = 1.1; if (freq < fmin) freq = fmin; if (freq > fmax) freq = fmax; @@ -73,7 +73,7 @@ class DSPIfilterOrtho { } - inline void setAP(float freq, float Q) // allpass + inline void setAP(t_float freq, t_float Q) // allpass { // prototype: H(s) = (1 - 2s/Qw0 + (s/w0)^2) / (1 + 2s/Qw0 + (s/w0)^2) @@ -83,46 +83,46 @@ class DSPIfilterOrtho { // prewarp for bilin transfo freq = bilin_prewarp(freq); - float zeta = 1.0f/Q; + t_float zeta = 1.0/Q; - DSPIcomplex p = bilin_stoz(DSPIcomplex(-zeta, (1.0f-zeta*zeta))*freq); - DSPIcomplex z = 1.0f / p; + DSPIcomplex p = bilin_stoz(DSPIcomplex(-zeta, (1.0-zeta*zeta))*freq); + DSPIcomplex z = 1.0 / p; setPoleZeroNormalized(p, z, DSPIcomplex(1,0)); } - inline void setLP(float freq, float Q) // low pass + inline void setLP(t_float freq, t_float Q) // low pass { // prototype: H(s) = 1 / (1 + 2s/Qw0 + (s/w0)^2) // the bilinear transform has 2 zeros at NY checkBounds(freq, Q); freq = bilin_prewarp(freq); - float zeta = 1/Q; + t_float zeta = 1/Q; - DSPIcomplex p = bilin_stoz(DSPIcomplex(-zeta, (1.0f-zeta*zeta))*freq); + DSPIcomplex p = bilin_stoz(DSPIcomplex(-zeta, (1.0-zeta*zeta))*freq); setPoleZeroNormalized(p, DSPIcomplex(-1, 0), DSPIcomplex(1, 0)); } - inline void setHP(float freq, float Q) // hi pass + inline void setHP(t_float freq, t_float Q) // hi pass { // prototype: H(s) = (s/w0)^2 / (1 + 2s/Qw0 + (s/w0)^2) // the bilinear transform has 2 zeros at DC checkBounds(freq, Q); freq = bilin_prewarp(freq); - float zeta = 1/Q; + t_float zeta = 1/Q; - DSPIcomplex p = bilin_stoz(DSPIcomplex(-zeta, (1.0f-zeta*zeta))*freq); + DSPIcomplex p = bilin_stoz(DSPIcomplex(-zeta, (1.0-zeta*zeta))*freq); setPoleZeroNormalized(p, DSPIcomplex(1, 0), DSPIcomplex(-1, 0)); } - inline void setBP(float freq, float Q) // band pass (1-allpass) + inline void setBP(t_float freq, t_float Q) // band pass (1-allpass) { // prototype: 1/2 * (1 - H_allpass(z)) setAP(freq, Q); - float h = -0.5f; + t_float h = -0.5; c0 *= h; c1 *= h; c2 *= h; @@ -130,46 +130,46 @@ class DSPIfilterOrtho { } - inline void setBR(float freq, float Q) // band reject + inline void setBR(t_float freq, t_float Q) // band reject { // prototype: H(s) = (1 - (s/w0)^2) / (1 + 2s/Qw0 + (s/w0)^2) checkBounds(freq, Q); // pole phasor - DSPIcomplex z = DSPIcomplex(2.0f * M_PI * freq); + DSPIcomplex z = DSPIcomplex(2.0 * M_PI * freq); // prewarp for bilin transfo freq = bilin_prewarp(freq); - float zeta = 1/Q; + t_float zeta = 1/Q; - DSPIcomplex p = bilin_stoz(DSPIcomplex(-zeta, (1.0f-zeta*zeta))*freq); + DSPIcomplex p = bilin_stoz(DSPIcomplex(-zeta, (1.0-zeta*zeta))*freq); setPoleZeroNormalized(p, z, DSPIcomplex(1,0)); } - inline void setHS(float freq, float gain) // low shelf + inline void setHS(t_float freq, t_float gain) // low shelf { // hi shelf = LP - g(LP-1) - float Q = M_SQRT2; + t_float Q = M_SQRT2; setLP(freq, Q); - c0 -= gain * (c0 - 1.0f); + c0 -= gain * (c0 - 1.0); c1 -= gain * (c1); c2 -= gain * (c2); } - inline void setLS(float freq, float gain) // low shelf + inline void setLS(t_float freq, t_float gain) // low shelf { // hi shelf = HP - g(HP-1) - float Q = M_SQRT2; + t_float Q = M_SQRT2; setHP(freq, Q); - c0 -= gain * (c0 - 1.0f); + c0 -= gain * (c0 - 1.0); c1 -= gain * (c1); c2 -= gain * (c2); } - inline void setEQ(float freq, float Q, float gain)// param EQ + inline void setEQ(t_float freq, t_float Q, t_float gain)// param EQ { // EQ = (1+A)/2 + (1-A)/2 AP - float a0 = 0.5f * (1.0f + gain); - float a1 = 0.5f * (1.0f - gain); + t_float a0 = 0.5 * (1.0 + gain); + t_float a1 = 0.5 * (1.0 - gain); setAP(freq, Q); c0 *= a1; c1 *= a1; @@ -186,8 +186,8 @@ class DSPIfilterOrtho { ar = a.r(); ai = a.i(); - c0 = 1.0f; - c1 = 2.0f * (a.r() - b.r()); + c0 = 1.0; + c1 = 2.0 * (a.r() - b.r()); c2 = (a.norm2() - b.norm2() - c1 * a.r()) / a.i(); } @@ -201,7 +201,7 @@ class DSPIfilterOrtho { { setPoleZero(a, b); DSPIcomplex invComplexGain = ((c-a)*(c-a.conj()))/((c-b)*(c-b.conj())); - float invGain = invComplexGain.norm(); + t_float invGain = invComplexGain.norm(); c0 *= invGain; c1 *= invGain; c2 *= invGain; @@ -212,12 +212,12 @@ class DSPIfilterOrtho { // one channel bang inline void Bang ( - float &input, - float &output + t_float &input, + t_float &output ) { - float d1t = ar * d1A + ai * d2A + input; - float d2t = ar * d2A - ai * d1A; + t_float d1t = ar * d1A + ai * d2A + input; + t_float d2t = ar * d2A - ai * d1A; output = c0 * input + c1 * d1A + c2 * d2A; d1A = d1t; d2A = d2t; @@ -227,13 +227,13 @@ class DSPIfilterOrtho { // a default s could be s = (1 - (.1)^(1/n)) inline void BangSmooth ( - float &input, // input ref - float &output, // output ref - float s // smooth pole + t_float &input, // input ref + t_float &output, // output ref + t_float s // smooth pole ) { - float d1t = s_ar * d1A + s_ai * d2A + input; - float d2t = s_ar * d2A - s_ai * d1A; + t_float d1t = s_ar * d1A + s_ai * d2A + input; + t_float d2t = s_ar * d2A - s_ai * d1A; s_ar += s * (ar - s_ar); s_ai += s * (ai - s_ai); output = s_c0 * input + s_c1 * d1A + s_c2 * d2A; @@ -247,16 +247,16 @@ class DSPIfilterOrtho { // two channel bang inline void Bang2 ( - float &input1, - float &input2, - float &output1, - float &output2 + t_float &input1, + t_float &input2, + t_float &output1, + t_float &output2 ) { - float d1tA = ar * d1A + ai * d2A + input1; - float d1tB = ar * d1B + ai * d2B + input2; - float d2tA = ar * d2A - ai * d1A; - float d2tB = ar * d2B - ai * d1B; + t_float d1tA = ar * d1A + ai * d2A + input1; + t_float d1tB = ar * d1B + ai * d2B + input2; + t_float d2tA = ar * d2A - ai * d1A; + t_float d2tB = ar * d2B - ai * d1B; output1 = c0 * input1 + d1A * c1 + d2A * c2; output2 = c0 * input2 + d1B * c1 + d2B * c2; d1A = d1tA; @@ -268,17 +268,17 @@ class DSPIfilterOrtho { // two channel bang smooth inline void Bang2Smooth ( - float &input1, - float &input2, - float &output1, - float &output2, - float s + t_float &input1, + t_float &input2, + t_float &output1, + t_float &output2, + t_float s ) { - float d1tA = s_ar * d1A + s_ai * d2A + input1; - float d1tB = s_ar * d1B + s_ai * d2B + input2; - float d2tA = s_ar * d2A - s_ai * d1A; - float d2tB = s_ar * d2B - s_ai * d1B; + t_float d1tA = s_ar * d1A + s_ai * d2A + input1; + t_float d1tB = s_ar * d1B + s_ai * d2B + input2; + t_float d2tA = s_ar * d2A - s_ai * d1A; + t_float d2tB = s_ar * d2B - s_ai * d1B; s_ar += s * (ar - s_ar); s_ai += s * (ai - s_ai); output1 = s_c0 * input1 + d1A * s_c1 + d2A * s_c2; @@ -296,7 +296,7 @@ class DSPIfilterOrtho { { // state data - float zero = 0.0f; + t_float zero = 0.0; d1A = DSPI_IS_DENORMAL(d1A) ? zero : d1A; d2A = DSPI_IS_DENORMAL(d2A) ? zero : d2A; @@ -310,11 +310,11 @@ class DSPIfilterOrtho { // smooth data - float dai = ai - s_ai; - float dar = ar - s_ar; - float dc0 = c0 - s_c0; - float dc1 = c1 - s_c1; - float dc2 = c2 - s_c2; + t_float dai = ai - s_ai; + t_float dar = ar - s_ar; + t_float dc0 = c0 - s_c0; + t_float dc1 = c1 - s_c1; + t_float dc2 = c2 - s_c2; s_ai = DSPI_IS_DENORMAL(dai) ? ai : s_ai; @@ -332,25 +332,25 @@ class DSPIfilterOrtho { private: // state data - float d1A; - float d2A; + t_float d1A; + t_float d2A; - float d1B; - float d2B; + t_float d1B; + t_float d2B; // pole data - float ai; - float s_ai; - float ar; - float s_ar; + t_float ai; + t_float s_ai; + t_float ar; + t_float s_ar; // zero data - float c0; - float s_c0; - float c1; - float s_c1; - float c2; - float s_c2; + t_float c0; + t_float s_c0; + t_float c1; + t_float s_c1; + t_float c2; + t_float s_c2; }; @@ -369,22 +369,22 @@ class DSPIfilterSeries{ biquad = new DSPIfilterOrtho[numberOfSections]; } - inline void setButterHP(float freq) + inline void setButterHP(t_float freq) { /* This member function computes the poles for a highpass butterworth filter. * The filter is transformed to the digital domain using a bilinear transform. * Every biquad section is normalized at NY. */ - float epsilon = .0001f; // stability guard - float min = 0.0f + epsilon; - float max = 0.5f - epsilon; + t_float epsilon = .0001; // stability guard + t_float min = 0.0 + epsilon; + t_float max = 0.5 - epsilon; if (freq < min) freq = min; if (freq > max) freq = max; // prewarp cutoff frequency - float omega = bilin_prewarp(freq); + t_float omega = bilin_prewarp(freq); DSPIcomplex NY(-1,0); //normalize at NY DSPIcomplex DC(1,0); //all zeros will be at DC @@ -403,7 +403,7 @@ class DSPIfilterSeries{ } - inline void setButterLP(float freq) + inline void setButterLP(t_float freq) { /* This member function computes the poles for a lowpass butterworth filter. * The filter is transformed to the digital domain using a bilinear transform. @@ -414,16 +414,16 @@ class DSPIfilterSeries{ */ - float epsilon = .0001f; // stability guard - float min = 0.0f + epsilon; - float max = 0.5f - epsilon; + t_float epsilon = .0001; // stability guard + t_float min = 0.0 + epsilon; + t_float max = 0.5 - epsilon; if (freq < min) freq = min; if (freq > max) freq = max; // prewarp cutoff frequency - float omega = bilin_prewarp(freq); + t_float omega = bilin_prewarp(freq); DSPIcomplex DC(1,0); //normalize at DC DSPIcomplex NY(-1,0); //all zeros will be at NY @@ -444,19 +444,19 @@ class DSPIfilterSeries{ for (int i=0; i