From 494a07a361fe4ee0e54f77468a976b1a77818770 Mon Sep 17 00:00:00 2001 From: Tom Schouten Date: Fri, 12 Sep 2003 22:26:57 +0000 Subject: creb 0.9.0 svn path=/trunk/externals/creb/; revision=956 --- include/dspi/DSPIcomplex.h | 2 +- include/extlib_util.h | 12 ++- include/filters.h | 226 +++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 238 insertions(+), 2 deletions(-) create mode 100644 include/filters.h (limited to 'include') diff --git a/include/dspi/DSPIcomplex.h b/include/dspi/DSPIcomplex.h index c0e2d63..ad3e041 100644 --- a/include/dspi/DSPIcomplex.h +++ b/include/dspi/DSPIcomplex.h @@ -22,7 +22,7 @@ #define DSPIcomplex_h #include -#include +#include class DSPIcomplex { diff --git a/include/extlib_util.h b/include/extlib_util.h index 6e4d94b..3f92ab0 100644 --- a/include/extlib_util.h +++ b/include/extlib_util.h @@ -30,5 +30,15 @@ /* convert milliseconds to 1-p, with p a real pole */ float milliseconds_2_one_minus_realpole(float time); + +typedef union +{ + unsigned int i; + float f; +} t_flint; + /* check if floating point number is denormal */ -#define IS_DENORMAL(f) (((*(unsigned int *)&(f))&0x7f800000) == 0) + +//#define IS_DENORMAL(f) (((*(unsigned int *)&(f))&0x7f800000) == 0) + +#define IS_DENORMAL(f) (((((t_flint)(f)).i) & 0x7f800000) == 0) diff --git a/include/filters.h b/include/filters.h new file mode 100644 index 0000000..e0d1c49 --- /dev/null +++ b/include/filters.h @@ -0,0 +1,226 @@ +/* this file contains a 37th attempt to write a general purpose iir filter toolset */ + +/* defined as inline functions with call by reference to enable compiler ref/deref optim */ + +/* the typedef */ +#ifndef T +#define T float +#endif + + +/* the prototype for a word */ +#define P static inline void +#define PP static void + + +/* the 'reference' arguments */ +#define A *a +#define B *b +#define C *c +#define D *d +#define X *x +#define Y *y +#define S *s + + +/* the opcodes */ + +/* add */ +P cadd (T X, T Y, T A, T B, T C, T D) { X = A + C; Y = B + D;} +P cadd2 (T A, T B, T C, T D) { A += C; B += D;} +P vcadd (T X, T A, T C) { cadd(x,x+1,a,a+1,c,c+1); } +P vcadd2 (T A, T C) { cadd2(a,a+1,c,c+1); } + + +/* mul */ +P cmul_r (T X, T A, T B, T C, T D) { X = A * C - B * D;} +P cmul_i (T Y, T A, T B, T C, T D) { Y = A * D + B * C;} +P cmul (T X, T Y, T A, T B, T C, T D) +{ + cmul_r (x, a, b, c, d); + cmul_i (y, a, b, c, d); +} +P cmul2 (T A, T B, T C, T D) +{ + T x = A; + T y = B; + cmul (&x, &y, a, b, c, d); + A = x; + B = y; +} + +P vcmul (T X, T A, T C) { cmul(x,x+1,a,a+1,c,c+1); } +P vcmul2 (T A, T C) { cmul2(a,a+1,c,c+1); } + + +/* norm */ +static inline float vcnorm(T X) { return hypot(x[0], x[1]); } + + + +/* swap */ +P vcswap(T Y, T X) +{ + float t[2] = {x[0], x[1]}; + x[0] = y[0]; + x[1] = y[1]; + y[0] = t[0]; + y[1] = t[1]; +} + + +/* inverse */ +P vcinv(T Y, T X) +{ + float scale = 1.0f / vcnorm(x); + y[0] = scale * x[0]; + y[1] = scale * x[1]; +} + +P vcinv1(T X) +{ + float scale = 1.0f / vcnorm(x); + x[0] *= scale; + x[1] *= scale; +} + +/* exp */ + +/* X = exp(Y) */ +P vcexp2 (T Y, T X) +{ + T r = exp(x[0]); + y[0] = cos (x[1]); + y[1] = sin (x[1]); + y[0] *= r; + y[1] *= r; +} + +P vcexp1 (T X) +{ + T y[2]; + vcexp2(y,x); + x[0] = y[0]; + x[1] = y[1]; +} + +/* + FILTERS + + the transfer function is defined in terms of the "engineering" + bilateral z-transform of the discrete impulse response h[n]. + + H(z) = Sum{n = -inf -> inf} h[n] z^(-n) + + a unit delay operating on a singnal S(z) is therefore + represented as z^(-1) S(z) + +*/ + + + + + + + +/* biquads */ + +/* biquad, orthogonal (poles & zeros), real in, out, state, pole, output */ +P biq_orth_r (T X, T Y, T S, T A, T C) +{ + Y = X + c[0] * s[0] - c[1] * s[1]; /* mind sign of c[1] */ + vcmul2(s, a); + S += X; +} + + +/* biquad, orthogonal, complex one-pole, with scaling */ + +/* complex one pole: (output = s[0] + is[1]): C / (1-A z^(-1)) */ + +P one_pole_complex (T X, T Y, T S, T A, T C) +{ + vcmul(y, s, a); + vcadd2(y, x); + s[0] = y[0]; + s[1] = y[1]; + vcmul(y, s, c); +} + +/* complex conj two pole: (output = s[0] : (Re(C) - Re(C*Conj(A))) / (1 - A z^(-1)) (1 - Conj(A) z^(-1)) */ + +P two_pole_complex_conj (T X, T Y, T S, T A, T C) +{ + vcmul2(s, a); + s[0] += x[0]; + y[0] = s[0] * c[0] - s[1] * c[1]; +} + + + +/* support functions for IIR filter design */ + +/* evaluate pole and allzero TF in z^-1 given the complex zeros/poles: + p(z) (or p(z)^-1) = \product (1-z_i z^-1) */ +PP eval_zero_poly(float *val, float *arg, float *zeros, int nb_zeros) +{ + int i; + float a[2] = {arg[0], arg[1]}; + vcinv1(a); + val[0] = 1.0f; + val[1] = 0.0f; + a[0] *= -1; + a[1] *= -1; + for (i=0; i