externals/creb/modules++ DSPI.h, NONE, 1.1 DSPIcomplex.h, NONE, 1.1 DSPIfilters.h, NONE, 1.1 filters.h, NONE, 1.1 - Pd-cvs

22 Aug 2006


      Update of /cvsroot/pure-data/externals/creb/modules++
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv23744/creb/modules++
Added Files:
    DSPI.h DSPIcomplex.h DSPIfilters.h filters.h 
Log Message:
fixed things so that all but one of the objects compile into a libdir
--- NEW FILE: DSPIcomplex.h ---
/*
 *   DSPIcomplex.h - Quick and dirty inline class for complex numbers 
 *   (mainly to compute filter poles/zeros, not to be used inside loops)
 *   Copyright (c) 2000 by Tom Schouten
 *
 *   This program is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation; either version 2 of the License, or
 *   (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program; if not, write to the Free Software
 *   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
#ifndef DSPIcomplex_h
#define DSPIcomplex_h
#include <math.h>
#include <iostream>
class DSPIcomplex
{
    public:
        inline DSPIcomplex() {_r = _i = 0;}
        inline DSPIcomplex(const float &a, const float &b) {setCart(a, b);}
        inline DSPIcomplex(const float &phasor) {setAngle(phasor);}
inline void setAngle(const float &angle) {_r = cos(angle); _i = sin(angle);}
        inline void setPolar(const float &phasor, const float &norm) 
        {_r = norm * cos(phasor); _i = norm * sin(phasor);}
        inline void setCart(const float &a, const float &b) {_r = a; _i = b;}
inline const float& r() const {return _r;}
        inline const float& i() const {return _i;}
inline float norm2() const {return _r*_r+_i*_i;}
        inline float norm() const {return sqrt(norm2());}
        inline void normalize() {float n = 1.0f / norm(); _r *= n; _i *= n;}
inline DSPIcomplex conj() const {return DSPIcomplex(_r, -_i);}
inline float angle() const {return atan2(_i, _r);}
inline DSPIcomplex operator+ (const DSPIcomplex &a) const
        {
            return DSPIcomplex(_r + a.r(), _i + a.i());
        }
        inline DSPIcomplex operator+ (float f) const
        {
            return DSPIcomplex(_r + f, _i);
        }
        inline DSPIcomplex operator- (const DSPIcomplex &a) const
        {
            return DSPIcomplex(_r - a.r(), _i - a.i());
        }
        inline DSPIcomplex operator- (float f) const
        {
            return DSPIcomplex(_r - f, _i);
        }
inline DSPIcomplex operator* (const DSPIcomplex &a) const 
        {
            return DSPIcomplex(_r * a.r() - _i * a.i(), _i * a.r() + _r * a.i());
        }
        inline DSPIcomplex operator* (float f) const
        {
            return DSPIcomplex(_r * f, _i * f);
        }
        inline DSPIcomplex operator/ (const DSPIcomplex &a) const 
        {
            float n_t = 1.0f / a.norm2();
            return DSPIcomplex(n_t * (_r * a.r() + _i * a.i()), n_t * (_i * a.r() - _r * a.i()));
        }
        inline DSPIcomplex operator/ (float f) const 
        {
            float n_t = 1.0f / f;
            return DSPIcomplex(n_t * _r, n_t * _i);
        }
inline friend std::ostream& operator<< (std::ostream& o, DSPIcomplex& a)
        {
            return o << "(" << a.r() << "," << a.i() << ")";
        }
inline friend DSPIcomplex operator+ (float f, DSPIcomplex& a)
        {
            return(DSPIcomplex(a.r() + f, a.i()));
        }
inline friend DSPIcomplex operator- (float f, DSPIcomplex& a)
        {
            return(DSPIcomplex(f - a.r(), - a.i()));
        }
inline friend DSPIcomplex operator/ (float f, DSPIcomplex& a)
        {
            return(DSPIcomplex(f,0) / a);
        }
// ????
        inline friend DSPIcomplex operator* (float f, DSPIcomplex& a)
        {
            return(DSPIcomplex(f*a.r(), f*a.i()));
        }
inline DSPIcomplex& operator *= (float f)
        {
            _r *= f;
            _i *= f;
            return *this;
        }
inline DSPIcomplex& operator /= (float f)
        {
            _r /= f;
            _i /= f;
            return *this;
        }
inline DSPIcomplex& operator *= (DSPIcomplex& a)
        {
            float r_t = _r * a.r() - _i * a.i();
                   _i = _r * a.i() + _i * a.r();
                   _r = r_t;
return *this;
        }
inline DSPIcomplex& operator /= (DSPIcomplex& a)
        {
            float n_t = a.norm2();
            float r_t = n_t * (_r * a.r() + _i * a.i());
                   _i = n_t * (_i * a.r() - _r * a.i());
                   _r = r_t;
return *this;
        }
float _r;
        float _i;
};
// COMPLEX LOG
inline DSPIcomplex dspilog(DSPIcomplex a) /* complex log */
{
    float r_t = log(a.norm());
    float i_t = a.angle();
    return DSPIcomplex(r_t, i_t);
}
// COMPLEX EXP
inline DSPIcomplex dspiexp(DSPIcomplex a) /* complex exp */
{
    return (DSPIcomplex(a.i()) * exp(a.r()));
}
// BILINEAR TRANSFORM analog -> digital
inline DSPIcomplex bilin_stoz(DSPIcomplex a)
{
    DSPIcomplex a2 = a * 0.5f;
    return((1.0f + a2)/(1.0f - a2));
}
// BILINEAR TRANSFORM digital -> analog
inline DSPIcomplex bilin_ztos(DSPIcomplex a)
{
    return ((a - 1.0f) / (a + 1.0f))*2.0f;
}
// not really a complex function but a nice complement to the bilinear routines
inline float bilin_prewarp(float freq)
{
    return 2.0f * tan(M_PI * freq);
}
#endif //DSPIcomplex_h
--- NEW FILE: DSPI.h ---
#ifndef DSPI_h
#define DSPI_h
#define DSPImin(x,y)			(((x)<(y)) ? (x) : (y))
#define DSPImax(x,y)			(((x)>(y)) ? (x) : (y))
#define DSPIclip(min, x, max)	(DSPImin(DSPImax((min), (x)), max))
// test if floating point number is denormal
#define DSPI_IS_DENORMAL(f) (((*(unsigned int *)&(f))&0x7f800000) == 0)
// test if almost denormal, choose whichever is fastest
#define DSPI_IS_ALMOST_DENORMAL(f) (((*(unsigned int *)&(f))&0x7f800000) < 0x08000000)
//#define DSPI_IS_ALMOST_DENORMAL(f) (fabs(f) < 3.e-34)
#endif
--- NEW FILE: DSPIfilters.h ---
/*
 *   DSPIfilters.h - Inline classes for biquad filters 
 *   Copyright (c) 2000 by Tom Schouten
 *
 *   This program is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation; either version 2 of the License, or
 *   (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program; if not, write to the Free Software
 *   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
#ifndef DSPIfilters_h
#define DSPIfilters_h
#include "DSPIcomplex.h"
#include "DSPI.h"
//#include <stdio.h>
/* orthogonal biquad */
class DSPIfilterOrtho {
    public:
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;}
/*
         *		Biquad filters remarks
         *
         *		Q is defined with reference to the analog prototype:
         *		poles/zeros = w0 * (1/Q +- sqrt(1 - 1/Q^2))
         *
         *		the num/den polynomial then has the form:
         *		1 + 2s/Qw0 + (s/w0)^2
         *
         *		if Q < 1 => real valued poles/zeros
         *		if Q > 1 => complex values poles/zeros
         *		if Q = inf => imaginary poles/zeros
         *		if Q = sqrt(2) => 'maximally flat' poles/zeros
         *
         *		the analog prototypes are converted to the digital domain
         *		using the bilinear transform. hence the prewarping.
         */
// make sure freq and Q are positive and within bounds
        inline void checkBounds(float &freq, 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;
if (freq < fmin) freq = fmin; 
            if (freq > fmax) freq = fmax;
if (Q < Qmin) Q = Qmin;
}
inline void setAP(float freq, float Q) // allpass
        {
// prototype: H(s) = (1 - 2s/Qw0 + (s/w0)^2) / (1 + 2s/Qw0 + (s/w0)^2)
            // s_p = - s_z (analog: symmetric wrt. im axis)
            // z_p = 1/z_z (digiatl: summ wrt. unit circle)
            checkBounds(freq, Q);
// prewarp for bilin transfo
            freq = bilin_prewarp(freq);
            float zeta = 1.0f/Q;
DSPIcomplex p = bilin_stoz(DSPIcomplex(-zeta, (1.0f-zeta*zeta))*freq);
            DSPIcomplex z = 1.0f / p;
            setPoleZeroNormalized(p, z, DSPIcomplex(1,0));
}
        inline void setLP(float freq, 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;
DSPIcomplex p = bilin_stoz(DSPIcomplex(-zeta, (1.0f-zeta*zeta))*freq);
            setPoleZeroNormalized(p, DSPIcomplex(-1, 0), DSPIcomplex(1, 0));
}
        inline void setHP(float freq, 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;
DSPIcomplex p = bilin_stoz(DSPIcomplex(-zeta, (1.0f-zeta*zeta))*freq);
            setPoleZeroNormalized(p, DSPIcomplex(1, 0), DSPIcomplex(-1, 0));
}
inline void setBP(float freq, float Q) // band pass (1-allpass)
        {
            // prototype: 1/2 * (1 - H_allpass(z))
            setAP(freq, Q);
            float h = -0.5f;
            c0 *= h;
            c1 *= h;
            c2 *= h;
            c0 -= h;
}
inline void setBR(float freq, 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);
// prewarp for bilin transfo
            freq = bilin_prewarp(freq);
            float zeta = 1/Q;
DSPIcomplex p = bilin_stoz(DSPIcomplex(-zeta, (1.0f-zeta*zeta))*freq);
            setPoleZeroNormalized(p, z, DSPIcomplex(1,0)); 
        }
inline void setHS(float freq, float gain) // low shelf
        {
            // hi shelf = LP - g(LP-1)
            float Q = M_SQRT2;
            setLP(freq, Q); 
            c0 -= gain * (c0 - 1.0f);
            c1 -= gain * (c1);
            c2 -= gain * (c2);
        }
inline void setLS(float freq, float gain) // low shelf
        {
            // hi shelf = HP - g(HP-1)
            float Q = M_SQRT2;
            setHP(freq, Q); 
            c0 -= gain * (c0 - 1.0f);
            c1 -= gain * (c1);
            c2 -= gain * (c2);
        }
        inline void setEQ(float freq, float Q, 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);
            setAP(freq, Q);
            c0 *= a1;
            c1 *= a1;
            c2 *= a1;
            c0 += a0;
        }
inline void setPoleZero
        (
            const DSPIcomplex& a, // pole
            const DSPIcomplex& b  // zero
        )
        {
            ar = a.r();
            ai = a.i();
c0 = 1.0f;
            c1 = 2.0f * (a.r() - b.r());
            c2 = (a.norm2() - b.norm2() - c1 * a.r()) / a.i();
        }
inline void setPoleZeroNormalized
        (
            const DSPIcomplex& a, // pole
            const DSPIcomplex& b, // zero
            const DSPIcomplex& c  // gain = 1 at this freq
        )
        {
            setPoleZero(a, b);
            DSPIcomplex invComplexGain = ((c-a)*(c-a.conj()))/((c-b)*(c-b.conj()));
            float invGain = invComplexGain.norm();
            c0 *= invGain;
            c1 *= invGain;
            c2 *= invGain;
}
// one channel bang
        inline void Bang
        (
            float &input, 
            float &output
        )
        {
            float d1t = ar * d1A + ai * d2A + input;
            float d2t = ar * d2A - ai * d1A;
            output = c0 * input + c1 * d1A + c2 * d2A;
            d1A = d1t;
            d2A = d2t;
        }
// one channel bang smooth
        // 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
        )
        {
            float d1t = s_ar * d1A + s_ai * d2A + input;
            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;
            d1A = d1t;
            d2A = d2t;
            s_c0 += s * (c0 - s_c0);
            s_c1 += s * (c1 - s_c1);
            s_c2 += s * (c2 - s_c2);
        }
// two channel bang
        inline void Bang2
        (
            float &input1, 
            float &input2, 
            float &output1, 
            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;
            output1 = c0 * input1 + d1A * c1 + d2A * c2;
            output2 = c0 * input2 + d1B * c1 + d2B * c2;
            d1A = d1tA;
            d2A = d2tA;
            d1B = d1tB;
            d2B = d2tB;
        }
// two channel bang smooth
        inline void Bang2Smooth
        (
            float &input1, 
            float &input2, 
            float &output1, 
            float &output2,
            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;
            s_ar += s * (ar - s_ar);
            s_ai += s * (ai - s_ai);
            output1 = s_c0 * input1 + d1A * s_c1 + d2A * s_c2;
            output2 = s_c0 * input2 + d1B * s_c1 + d2B * s_c2;
            d1A = d1tA;
            d2A = d2tA;
            d1B = d1tB;
            d2B = d2tB;
            s_c0 += s * (c0 - s_c0);
            s_c1 += s * (c1 - s_c1);
            s_c2 += s * (c2 - s_c2);
        }
inline void killDenormals()
    {
// state data
        float zero = 0.0f;
d1A = DSPI_IS_DENORMAL(d1A) ? zero : d1A;
        d2A = DSPI_IS_DENORMAL(d2A) ? zero : d2A;
        d1B = DSPI_IS_DENORMAL(d1B) ? zero : d1B;
        d2B = DSPI_IS_DENORMAL(d2B) ? zero : d2B;
/* test on athlon showed nuking smooth data does not
         * present a noticable difference in performance however
         * nuking state data is really necessary
// 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;
s_ai = DSPI_IS_DENORMAL(dai) ? ai : s_ai;
        s_ar = DSPI_IS_DENORMAL(dar) ? ar : s_ar;
        s_c0 = DSPI_IS_DENORMAL(dc0) ? c0 : s_c0;
        s_c1 = DSPI_IS_DENORMAL(dc0) ? c1 : s_c1;
        s_c2 = DSPI_IS_DENORMAL(dc0) ? c2 : s_c2;
*/
}
private:
// state data
        float d1A;
        float d2A;
float d1B;
        float d2B;
// pole data
        float ai;
        float s_ai;
        float ar;
        float s_ar;
// zero data
        float c0;
        float s_c0;
        float c1;
        float s_c1;
        float c2;
        float s_c2;
};
class DSPIfilterSeries{
    public:
        inline DSPIfilterSeries() {DSPIfilterSeries(1);}
        inline ~DSPIfilterSeries() {delete [] biquad;};
inline DSPIfilterSeries(int numberOfSections)
        {
            // create a set of biquads
            sections = numberOfSections;
            biquad = new DSPIfilterOrtho[numberOfSections];
        }
inline void setButterHP(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;
if (freq < min) freq = min; 
            if (freq > max) freq = max;
// prewarp cutoff frequency
            float omega = bilin_prewarp(freq);
DSPIcomplex NY(-1,0);  //normalize at NY
            DSPIcomplex DC(1,0);  //all zeros will be at DC
            DSPIcomplex pole( (2*sections + 1) * M_PI / (4*sections)); // first pole of lowpass filter with omega == 1
            DSPIcomplex pole_inc(M_PI / (2*sections)); // phasor to get to next pole, see Porat p. 331
for (int i=0; i<sections; i++)
            {
                // setup the biquad with the computed pole and zero and unit gain at NY
                biquad[i].setPoleZeroNormalized(
                    bilin_stoz(omega/pole),	// LP -> HP -> digital transfo 
                    DC, 					// all zeros at DC
                    NY);					// normalized (gain == 1) at NY
                pole *= pole_inc;			// compe next (lowpass) pole
            }
}
inline void setButterLP(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.
             *  Every biquad section is normalized at DC.
             *  Doing it this way, only the pole locations need to be transformed.
             *  The constant gain factor can be computed by setting the DC gain of every section to 1.
             *  An analog butterworth is all-pole, meaning the bilinear transform has all zeros at -1
             */
float epsilon = .0001f; // stability guard
            float min = 0.0f + epsilon;
            float max = 0.5f - epsilon;
if (freq < min) freq = min; 
            if (freq > max) freq = max;
// prewarp cutoff frequency
            float omega = bilin_prewarp(freq);
DSPIcomplex DC(1,0);  //normalize at DC
            DSPIcomplex NY(-1,0); //all zeros will be at NY
            DSPIcomplex pole( (2*sections + 1) * M_PI / (4*sections)); 
            pole *= omega; // first pole, see Porat p. 331
            DSPIcomplex pole_inc(M_PI / (2*sections)); // phasor to get to next pole, see Porat p. 331
for (int i=0; i<sections; i++)
            {
                // setup the biquad with the computed pole and zero and unit gain at DC
                biquad[i].setPoleZeroNormalized(bilin_stoz(pole), NY, DC);
                pole *= pole_inc;
            }
        }
inline void resetState()
        {
            for (int i=0; i<sections; i++) biquad[i].resetState();
        }
inline void Bang(float &input, float &output)
        {
            float x = input;
            for (int i=0; i<sections; i++)
            {
                biquad[i].Bang(x, x);
            }
            output = x;
        }
        inline void Bang2(float &input1, float &input2, float &output1, float &output2)
        {
            float x = input1;
            float y = input2;
            for (int i=0; i<sections; i++)
            {
                biquad[i].Bang2(x, y, x, y);
            }
            output1 = x;
            output2 = y;
        }
inline void BangSmooth(float &input, float &output, float s)
        {
            float x = input;
            for (int i=0; i<sections; i++)
            {
                biquad[i].BangSmooth(x, x, s);
            }
            output = x;
        }
        inline void Bang2(float &input1, float &input2, float &output1, float &output2, float s)
        {
            float x = input1;
            float y = input2;
            for (int i=0; i<sections; i++)
            {
                biquad[i].Bang2Smooth(x, y, x, y, s);
            }
            output1 = x;
            output2 = y;
        }
private:
        int sections;
        DSPIfilterOrtho *biquad;
        float gain;
};
#endif //DSPIfilters_h
--- NEW FILE: filters.h ---
/* 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<nb_zeros; i++){
    float t[2];
    vcmul(t, a, zeros + 2*i);
    t[0] += 1.0f;
    vcmul2(val, t);
    }
}
PP eval_pole_poly(float *val, float *arg, float *poles, int nb_poles)
{
    eval_zero_poly(val, arg, poles, nb_poles);
    vcinv1(val);
}
/* since it's more efficient to store half of the poles for a real impulse
   response, these functions compute p(z) conj(p(conj(z)))  */
PP eval_conj_zero_poly(float *val, float *arg, float *zeros, int nb_zeros)
{
    float t[2];
    float a[2] = {arg[0], arg[1]};
    eval_zero_poly(t, a, zeros, nb_zeros);
    a[1] *= -1;
    eval_zero_poly(val, a, zeros, nb_zeros);
    val[1] *= -1;
    vcmul2(val, t);
}
PP eval_conj_pole_poly(float *val, float *arg, float *poles, int nb_poles)
{
    eval_conj_zero_poly(val, arg, poles, nb_poles);
    vcinv1(val);
}
PP eval_conj_pole_zero_ratfunc(float *val, float *arg, float *poles, float *zeros, int nb_poles, int nb_zeros)
{
    float t[2];
    eval_conj_zero_poly(t, arg, zeros, nb_zeros);
    eval_conj_pole_poly(val, arg, poles, nb_zeros);
    vcmul2(val, t);
}
/* bandlimited IIR impulse:
* design analog butterworth filter
   * obtain the partial fraction expansion of the transfer function
   * determine the state increment as a function of fractional delay of impulse location 
   (sample the impulse response)
*/