path: root/lib/fidlib/fidmkf.h

//
//      mkfilter-derived code
//      ---------------------
//
//        Copyright (c) 2002-2004 Jim Peters <http://uazu.net/>.  This
//        file is released under the GNU Lesser General Public License
//        (LGPL) version 2.1 as published by the Free Software
//        Foundation.  See the file COPYING_LIB for details, or visit
//        <http://www.fsf.org/licenses/licenses.html>.
//
//	This is all code derived from 'mkfilter' by Tony Fisher of the
//	University of York.  I've rewritten it all in C, and given it
//	a thorough overhaul, so there is actually none of his code
//	here any more, but it is all still very strongly based on the
//	algorithms and techniques that he used in 'mkfilter'.
//
//	For those who didn't hear, Tony Fisher died in February 2000
//	at the age of 43.  See his web-site for information and a
//	tribute:
//
//        http://www-users.cs.york.ac.uk/~fisher/
//        http://www-users.cs.york.ac.uk/~fisher/tribute.html
//
//      The original C++ sources and the rest of the mkfilter tool-set
//      are still available from his site:
//
//        http://www-users.cs.york.ac.uk/~fisher/mkfilter/
//
//	
//	I've made a number of improvements and changes whilst
//	rewriting the code in C.  For example, I halved the
//	calculations required in designing the filters by storing only
//	one complex pole/zero out of each conjugate pair.  This also
//	made it much easier to output the filter as a list of
//	sub-filters without lots of searching around to match up
//	conjugate pairs later on.  Outputting as a list of subfilters
//	permits greater accuracy in calculation of the response, and
//	also in the execution of the final filter.  Also, some FIR
//	coefficients can be marked as 'constant', allowing optimised
//	routines to be generated for whole classes of filters, with
//	just the variable coefficients filled in at run-time.
//
//	On the down-side, complex numbers are not portably available
//	in C before C99, so complex calculations here are done on
//	double[] arrays with inline functions, which ends up looking
//	more like assembly language than C.  Never mind.
//

//
//	LEGAL STUFF
//	-----------
//
//	Tony Fisher released his software on his University of York
//	pages for free use and free download.  The software itself has
//	no licence terms attached, nor copyright messages, just the
//	author's name, E-mail address and date.  Nor are there any
//	licence terms indicated on the website.  I understand that
//	under the Berne convention copyright does not have to be
//	claimed explicitly, so these are in fact copyright files by
//	legal default.  However, the intention was obviously that
//	these files should be used by others.
//
//	None of this really helps, though, if we're going to try to be
//	100% legally correct, so I wrote to Anthony Moulds who is the
//	contact name on Tony Fisher's pages now.  I explained what I
//	planned to do with the code, and he answered as follows:
//
//	(Note that I was planning to use it 'as-is' at that time,
//	rather than rewrite it as I have done now)
//
//      > To: "Jim Peters" <jim@uazu.net>
//      > From: "Anthony Moulds" <anthony@cs.york.ac.uk>
//      > Subject: RE: mkfilter source
//      > Date: Tue, 29 Oct 2002 15:30:19 -0000
//      > 
//      > Hi Jim,
//      > 
//      > Thanks for your email.
//      > 
//      > The University will be happy to let you use Dr Fisher's mkfilter
//      > code since your intention is not to profit financially from his work.
//      > 
//      > It would be nice if in some way you could acknowledge his contribution.
//      > 
//      > Best wishes and good luck with your work,
//      > 
//      > Anthony Moulds
//      > Senior Experimental Officer, 
//      > Computer Science Department,  University of York,
//      > York, England, UK. Tel: 44(0)1904 434758  Fax: 44(0)19042767
//      > ============================================================
//      > 
//      > 
//      > > -----Original Message-----
//      > > From: Jim Peters [mailto:jim@uazu.net]
//      > > Sent: Monday, October 28, 2002 12:36 PM
//      > > To: anthony@cs.york.ac.uk
//      > > Subject: mkfilter source
//      > > 
//      > > 
//      > > I'm very sorry to hear (rather late, I know) that Tony Fisher died --
//      > > I've always gone straight to the filter page, rather than through his
//      > > home page.  I hope his work remains available for the future.
//      > > 
//      > > Anyway, the reason I'm writing is to clarify the status of the
//      > > mkfilter source code.  Because copyright is not claimed on the web
//      > > page nor in the source distribution, I guess that Tony's intention was
//      > > that this code should be in the public domain.  However, I would like
//      > > to check this now to avoid complications later.
//      > > 
//      > > I am using his code, modified, to provide a library of filter-design
//      > > routines for a GPL'd filter design app, which is not yet released.
//      > > The library could also be used standalone, permitting apps to design
//      > > filters at run-time rather than using hard-coded compile-time filters.
//      > > My interest in filters is as a part of my work on the OpenEEG project
//
//	So this looks pretty clear to me.  I am not planning to profit
//	from the work, so everything is fine with the University.  I
//	guess others might profit from the work, indirectly, as with
//	any free software release, but so long as I don't, we're fine.
//
//	I hope this is watertight enough for Debian/etc.  Otherwise
//	I'll have to go back to Anthony Moulds for clarification.
//
//	Even though there is no code cut-and-pasted from 'mkfilter'
//	here, it is all very obviously based on that code, so it
//	probably counts as a derived work -- although as ever "I Am
//	Not A Lawyer".
//

#ifndef FIDMK_H
#define FIDMK_H

#ifndef T_MSVC
 #ifdef HUGE_VAL
  #define INF HUGE_VAL
 #else
  #define INF (1.0/0.0)
 #endif
#endif

//Hacks for crappy linker error in MSVC... - Albert
#ifdef T_MSVC
 #undef HUGE_VAL
 #define HUGE_VAL 1.797693134862315E+308
 #define INF HUGE_VAL
#endif

#define TWOPI (2*M_PI)


//
//	Complex square root: aa= aa^0.5
//

STATIC_INLINE double
my_sqrt(double aa) {
   return aa <= 0.0 ? 0.0 : sqrt(aa);
}

// 'csqrt' clashes with builtin in GCC 4, so call it 'c_sqrt'
STATIC_INLINE void 
c_sqrt(double *aa) {
   double mag= hypot(aa[0], aa[1]);
   double rr= my_sqrt((mag + aa[0]) * 0.5);
   double ii= my_sqrt((mag - aa[0]) * 0.5);
   if (aa[1] < 0.0) ii= -ii;
   aa[0]= rr;
   aa[1]= ii;
}

//
//	Complex imaginary exponent: aa= e^i.theta
//

STATIC_INLINE void 
cexpj(double *aa, double theta) {
   aa[0]= cos(theta);
   aa[1]= sin(theta);
}

//
//	Complex exponent: aa= e^aa
//

// 'cexp' clashes with builtin in GCC 4, so call it 'c_exp'
STATIC_INLINE void 
c_exp(double *aa) {
   double mag= exp(aa[0]);
   aa[0]= mag * cos(aa[1]);
   aa[1]= mag * sin(aa[1]);
}

//
//	Context for generating filters.
//
//	Note that the poles and zeros are stored in a strange way.
//	Rather than storing both a pole (or zero) and its complex
//	conjugate, I'm storing just one of the pair.  Also, for real
//	poles, I'm not storing the imaginary part (which is zero).
//	This results in a list of numbers exactly half the length you
//	might otherwise expect.  However, since some of these numbers
//	are in pairs, and some are single, we need a separate flag
//	array to indicate which is which.  poltyp[] serves this
//	purpose.  An entry is 1 if the corresponding offset is a real
//	pole, or 2 if it is the first of a pair of values making up a
//	complex pole.  The second value of the pair has an entry of 0
//	attached.  (Similarly for zeros in zertyp[])
//

#define MAXPZ 64 

struct mk_filter_context {
   int n_pol;		// Number of poles
   double pol[MAXPZ];	// Pole values (see above)
   char poltyp[MAXPZ];	// Pole value types: 1 real, 2 first of complex pair, 0 second
   int n_zer;		// Same for zeros ...
   double zer[MAXPZ];
   char zertyp[MAXPZ];
};

//
//	Pre-warp a frequency
//

STATIC_INLINE double 
prewarp(double val) {
   return tan(val * M_PI) / M_PI;
}


//
//	Bessel poles; final one is a real value for odd numbers of
//	poles
//

static double bessel_1[]= {
 -1.00000000000e+00
};

static double bessel_2[]= {
 -1.10160133059e+00, 6.36009824757e-01,
};

static double bessel_3[]= {
 -1.04740916101e+00, 9.99264436281e-01,
 -1.32267579991e+00,
};

static double bessel_4[]= {
 -9.95208764350e-01, 1.25710573945e+00,
 -1.37006783055e+00, 4.10249717494e-01,
};

static double bessel_5[]= {
 -9.57676548563e-01, 1.47112432073e+00,
 -1.38087732586e+00, 7.17909587627e-01,
 -1.50231627145e+00,
};

static double bessel_6[]= {
 -9.30656522947e-01, 1.66186326894e+00,
 -1.38185809760e+00, 9.71471890712e-01,
 -1.57149040362e+00, 3.20896374221e-01,
};

static double bessel_7[]= {
 -9.09867780623e-01, 1.83645135304e+00,
 -1.37890321680e+00, 1.19156677780e+00,
 -1.61203876622e+00, 5.89244506931e-01,
 -1.68436817927e+00,
};

static double bessel_8[]= {
 -8.92869718847e-01, 1.99832584364e+00,
 -1.37384121764e+00, 1.38835657588e+00,
 -1.63693941813e+00, 8.22795625139e-01,
 -1.75740840040e+00, 2.72867575103e-01,
};

static double bessel_9[]= {
 -8.78399276161e-01, 2.14980052431e+00,
 -1.36758830979e+00, 1.56773371224e+00,
 -1.65239648458e+00, 1.03138956698e+00,
 -1.80717053496e+00, 5.12383730575e-01,
 -1.85660050123e+00,
};

static double bessel_10[]= {
 -8.65756901707e-01,