r3972 - trunk/bse
- From: timj svn gnome org
- To: svn-commits-list gnome org
- Subject: r3972 - trunk/bse
- Date: Sun, 15 Oct 2006 21:05:58 -0400 (EDT)
Author: timj
Date: 2006-10-15 21:05:27 -0400 (Sun, 15 Oct 2006)
New Revision: 3972
Modified:
trunk/bse/ChangeLog
trunk/bse/bseellipticfilter.c
trunk/bse/bseellipticfilter.h
Log:
Sun Oct 15 23:17:39 2006 Tim Janik <timj gtk org>
* bseellipticfilter.h:
* bseellipticfilter.c: include math.h instead of declaring standard
math.h functions. print debugging/informative output through an
extra function. beefed up ARRSIZ once more for high order filters.
added custom function to print the transfer function coefficients
for debugging. fixed reading from stdin to use fgets().
reverted granularity changes in filter table printing. renamed all
complex functions to not clash with complex.h.
added some workarounds with fixme marks to stay value (error)
compatible with former versions of the ellf program, these can be
removed once we have filter tests in place. removed unecessary
comments, streamlined others. removed obsolete configuration
logic and other unecessary preprocessor cruft. removed unused
declarations, moved many variable and constant declarations into
an inner scope. renamed heavily abbreviated functions. renamed
and recoded main() into iir_filter_design() which does read all
design parameters from a structure. added a new main with simply
stdin-read logic to call iir_filter_design().
for the moment, moved structure declarations for filter design
parameters and filter calculation state into bseellipticfilter.h.
Modified: trunk/bse/ChangeLog
===================================================================
--- trunk/bse/ChangeLog 2006-10-15 22:48:08 UTC (rev 3971)
+++ trunk/bse/ChangeLog 2006-10-16 01:05:27 UTC (rev 3972)
@@ -1,11 +1,25 @@
Sun Oct 15 23:17:39 2006 Tim Janik <timj gtk org>
+ * bseellipticfilter.h:
* bseellipticfilter.c: include math.h instead of declaring standard
math.h functions. print debugging/informative output through an
extra function. beefed up ARRSIZ once more for high order filters.
added custom function to print the transfer function coefficients
for debugging. fixed reading from stdin to use fgets().
- reverted granularity changes in filter table printing.
+ reverted granularity changes in filter table printing. renamed all
+ complex functions to not clash with complex.h.
+ added some workarounds with fixme marks to stay value (error)
+ compatible with former versions of the ellf program, these can be
+ removed once we have filter tests in place. removed unecessary
+ comments, streamlined others. removed obsolete configuration
+ logic and other unecessary preprocessor cruft. removed unused
+ declarations, moved many variable and constant declarations into
+ an inner scope. renamed heavily abbreviated functions. renamed
+ and recoded main() into iir_filter_design() which does read all
+ design parameters from a structure. added a new main with simply
+ stdin-read logic to call iir_filter_design().
+ for the moment, moved structure declarations for filter design
+ parameters and filter calculation state into bseellipticfilter.h.
Sun Oct 15 22:52:53 2006 Tim Janik <timj gtk org>
Modified: trunk/bse/bseellipticfilter.c
===================================================================
--- trunk/bse/bseellipticfilter.c 2006-10-15 22:48:08 UTC (rev 3971)
+++ trunk/bse/bseellipticfilter.c 2006-10-16 01:05:27 UTC (rev 3972)
@@ -17,40 +17,16 @@
* otherwise) arising in any way out of the use of this software, even
* if advised of the possibility of such damage.
*/
+#include "bseellipticfilter.h"
+#define _ISOC99_SOURCE /* for INFINITY and NAN */
+#include <math.h>
#include <stdarg.h>
#include <stdio.h>
-#if 0
-extern double sqrt ( double );
-extern double fabs ( double );
-extern double log ( double );
-extern double tan ( double );
-extern double atan ( double );
-extern double floor ( double );
-extern double fabs ( double );
-extern double cabs ( cmplx * );
-extern double sqrt ( double );
-extern double atan2 ( double, double );
-extern double cos ( double );
-extern double sin ( double );
-extern double sqrt ( double );
-extern double frexp ( double, int * );
-extern double ldexp ( double, int );
-extern double exp ( double );
-extern double log ( double );
-extern double cos ( double );
-extern double sin ( double );
-extern double sqrt ( double );
-extern double fabs ( double );
-extern double asin ( double );
-extern double atan ( double );
-extern double atan2 ( double, double );
-extern double pow ( double, double );
-double sqrt(), fabs(), sin(), cos(), asin(), tanh();
-double sinh(), cosh(), atan(), exp();
-#else
-#include <math.h>
-#endif
+#include <stdlib.h>
+
+#define ARRSIZ (300) /* size of arrays used to store coefficients */
+
static void __attribute__ ((__format__ (__printf__, 1, 2)))
VERBOSE (const char *format,
...)
@@ -60,7 +36,6 @@
va_start (args, format);
vsnprintf (buffer, sizeof (buffer), format, args);
va_end (args);
- // printf (" #V-BSEmgc: %s", buffer);
printf ("# %s", buffer);
fflush (stdout);
}
@@ -68,47 +43,46 @@
#define EVERBOSE VERBOSE
//#define EVERBOSE(...) do{}while(0)
-/* === ellf.doc - start === */
-/* ellf.c
- * This program calculates design coefficients for
+#if 0 // FIXME: increase precision by using:
+
+//#include "bseieee754.h"
+#define PI (3.141592653589793238462643383279502884197) // pi
+#define BSE_PI_DIV_2 (1.570796326794896619231321691639751442099) // pi/2
+#define BSE_DOUBLE_EPSILON (1.1102230246251565404236316680908203125e-16) /* 2^-53, round-off error at 1.0 */
+#define BSE_DOUBLE_MAX_NORMAL (1.7976931348623157e+308) /* 7fefffff ffffffff, 2^1024 * (1 - BSE_DOUBLE_EPSILON) */
+
+
+
+
+#define DECIBELL_FACTOR (4.3429448190325182765112891891661) /* 10.0 / ln (10.0) */
+#define MACHEP (BSE_DOUBLE_EPSILON) /* the machine roundoff error */
+// #define PI /* PI is defined in bseieee754.h */
+#define PIO2 (BSE_PI_DIV_2) /* pi/2 */
+#define MAXNUM (BSE_DOUBLE_MAX_NORMAL) /* 2**1024*(1-MACHEP) */
+static void init_constants (void) {}
+
+#else
+
+static double DECIBELL_FACTOR = -1;
+static void
+init_constants (void)
+{
+ DECIBELL_FACTOR = 10.0/log(10.0);
+}
+static const double MAXNUM = 1.79769313486231570815E308; /* 2**1024*(1-MACHEP) */
+static const double PI = 3.14159265358979323846; /* pi */
+static const double PIO2 = 1.57079632679489661923; /* pi/2 */
+static const double MACHEP = 1.11022302462515654042E-16; /* 2**-53 */
+#endif
+
+static double fixme2local_a, fixme2local_k, fixme2local_m;
+static double fixme2local_1; // k
+static double fixme2local_2; // u
+
+/* This code calculates design coefficients for
* digital filters of the Butterworth, Chebyshev, or
* elliptic varieties.
- *
- *
- *
- * Usage:
- *
- * Inputs are entered by keyboard, or are redirected to come from
- * a command file, as follows:
- *
- * Kind of filter (1: Butterworth, 2: Chebyshev, 3: Elliptic,
- * 0: exit to monitor)
- *
- * Shape of filter (1: low pass, 2: band pass, 3: high pass,
- * 4: band reject, 0: exit to monitor)
- *
- * Order of filter (an integer)
- *
- * Passband ripple (peak to peak decibels)
- *
- * Sampling frequency (Hz)
- *
- * Passband edge frequency (Hz)
- *
- * Second passband edge frequency (for band pass or reject filters)
- *
- * Stop band edge frequency (Hz)
- * or stop band attenuation (entered as -decibels)
- *
- * The "exit to monitor" type 0 may be used to terminate the
- * program when input is redirected to come from a command file.
- *
- * If your specification is illegal, e.g. the stop band edge
- * is in the middle of the passband, the program will make you
- * start over. However, it remembers and displays the last
- * value of each parameter entered. To use the same value, just
- * hit carriage return instead of typing it in again.
- *
+ *
* The program displays relevant pass band and stop band edge
* frequencies and stop band attenuation. The z-plane coefficients
* are printed in these forms:
@@ -120,8 +94,6 @@
* table of the frequency response of the filter. You can
* get a picture by reading the table into gnuplot.
*
- *
- *
* Filter design:
*
* The output coefficients of primary interest are shown as follows:
@@ -151,46 +123,6 @@
* Thus the two coefficients for the pole would normally be
* negated in a typical implementation of the filter.
*
- *
- *
- * Compilation:
- *
- * This program has been compiled successfully on many different
- * computers. See the accompanying output listing file ellf.ans,
- * for a set of correct answers. Use the batch file test.bat to
- * check your executable program. If the low pass and high pass
- * options work but the others don't, then examine your atan2()
- * function carefully for reversed arguments or perhaps an offest of
- * pi. On most systems, define ANSIC to be 1. This sets the
- * expected atan2() arguments but does not otherwise imply anything
- * about the ANSI-ness of the program.
- *
- *
- *
- * Files:
- *
- * mconf.h system configuration include file
- * Be sure to define type of computer here!
- * cmplx.c complex arithmetic subroutine package
- * ellf.ans right answer file for some elliptic filters
- * ellf.que elliptic filter questions
- * ellf.c main program
- * ellf.doc this file
- * ellf.mak Microsoft MSDOS makefile
- * ellfu.mak Unix makefile
- * ellik.c incomplete elliptic integral of the first kind
- * ellpe.c complete elliptic integral of the second kind
- * ellpj.c Jacobian Elliptic Functions
- * ellpk.c complete elliptic integral of the first kind
- * makefile Unix makefile
- * mtherr.c common math function error handler
- * polevl.c evaluates polynomials
- * test.bat batch file to run a test
- * descrip.mms VAX makefile
- * ellf.opt VAX makefile
- * testvax.bat VAX test
- *
- *
* References:
*
* A. H. Gray, Jr., and J. D. Markel, "A Computer Program for
@@ -204,482 +136,235 @@
* M. Abramowitz and I. A. Stegun, eds., Handbook of Mathematical
* Functions, National Bureau of Standards AMS 55, 1964,
* Chapters 16 and 17
- *
- *
- * - Steve Moshier, December 1986
- * Last rev: November, 1992
*/
-/* === ellf.doc - end === */
-/* === mconf.h - start === */
-/* mconf.h
- *
- * Common include file for math routines
- *
- *
- *
- * SYNOPSIS:
- *
- * #include "mconf.h"
- *
- *
- *
- * DESCRIPTION:
- *
- * This file contains definitions for error codes that are
- * passed to the common error handling routine mtherr()
- * (which see).
- *
- * The file also includes a conditional assembly definition
- * for the type of computer arithmetic (IEEE, DEC, Motorola
- * IEEE, or UNKnown).
- *
- * For Digital Equipment PDP-11 and VAX computers, certain
- * IBM systems, and others that use numbers with a 56-bit
- * significand, the symbol DEC should be defined. In this
- * mode, most floating point constants are given as arrays
- * of octal integers to eliminate decimal to binary conversion
- * errors that might be introduced by the compiler.
- *
- * For little-endian computers, such as IBM PC, that follow the
- * IEEE Standard for Binary Floating Point Arithmetic (ANSI/IEEE
- * Std 754-1985), the symbol IBMPC should be defined. These
- * numbers have 53-bit significands. In this mode, constants
- * are provided as arrays of hexadecimal 16 bit integers.
- *
- * Big-endian IEEE format is denoted MIEEE. On some RISC
- * systems such as Sun SPARC, double precision constants
- * must be stored on 8-byte address boundaries. Since integer
- * arrays may be aligned differently, the MIEEE configuration
- * may fail on such machines.
- *
- * To accommodate other types of computer arithmetic, all
- * constants are also provided in a normal decimal radix
- * which one can hope are correctly converted to a suitable
- * format by the available C language compiler. To invoke
- * this mode, define the symbol UNK.
- *
- * An important difference among these modes is a predefined
- * set of machine arithmetic constants for each. The numbers
- * MACHEP (the machine roundoff error), MAXNUM (largest number
- * represented), and several other parameters are preset by
- * the configuration symbol. Check the file const.c to
- * ensure that these values are correct for your computer.
- *
- * Configurations NANS, INFINITIES, MINUSZERO, and DENORMAL
- * may fail on many systems. Verify that they are supposed
- * to work on your computer.
- */
-/* Constant definitions for math error conditions
- */
+/* --- Complex numeral --- */
+typedef struct {
+ double r; // real part
+ double i; // imaginary part
+} Complex;
-#define DOMAIN 1 /* argument domain error */
-#define SING 2 /* argument singularity */
-#define OVERFLOW 3 /* overflow range error */
-#define UNDERFLOW 4 /* underflow range error */
-#define TLOSS 5 /* total loss of precision */
-#define PLOSS 6 /* partial loss of precision */
+/* --- prototypes --- */
+static const Complex COMPLEX_ONE = {1.0, 0.0};
+static double Cabs (const Complex *z);
+static void Cadd (const Complex *a, const Complex *b, Complex *c);
+static void Cdiv (const Complex *a, const Complex *b, Complex *c);
+static void Cmov (const Complex *a, Complex *b);
+static void Cmul (const Complex *a, const Complex *b, Complex *c);
+static void Csqrt (const Complex *z, Complex *w);
+static void Csub (const Complex *a, const Complex *b, Complex *c);
+static double polevl (double x, const double coef[], int N);
+static double ellik (double phi, double m); // incomplete elliptic integral of the first kind
+static double ellpk (double x); // complete elliptic integral of the first kind
+static int ellpj (double u, double m, double *sn, double *cn, double *dn, double *ph); // Jacobian Elliptic Functions
+static int math_set_error (char *name, int code);
-#define EDOM 33
-#define ERANGE 34
-/* Complex numeral. */
-typedef struct
-{
- double r;
- double i;
-} cmplx;
+/* --- math errors --- */
+static int math_global_error = 0;
+#define MATH_ERROR_DOMAIN 1 /* argument domain error */
+#define MATH_ERROR_SING 2 /* argument singularity */
+#define MATH_ERROR_OVERFLOW 3 /* overflow range error */
+#define MATH_ERROR_UNDERFLOW 4 /* underflow range error */
+#define MATH_ERROR_TOTAL_LOSS 5 /* total loss of precision */
+#define MATH_ERROR_PARTIAL_LOSS 6 /* partial loss of precision */
-/* Type of computer arithmetic is
- * UNKnown arithmetic, invokes coefficients given in
- * normal decimal format. Beware of range boundary
- * problems (MACHEP, MAXLOG, etc. in const.c) and
- * roundoff problems in pow.c:
- * (Sun SPARCstation, i386)
- */
-
-/* Define to support tiny denormal numbers, else undefine. */
-#define DENORMAL 1
-
-/* Define to ask for infinity support, else undefine. */
-/* #define INFINITIES 1 */
-
-/* Define to ask for support of numbers that are Not-a-Number,
- else undefine. This may automatically define INFINITIES in some files. */
-/* #define NANS 1 */
-
-/* Define to distinguish between -0.0 and +0.0. */
-#define MINUSZERO 1
-
-/* Define 1 for ANSI C atan2() function
- See atan.c and clog.c. */
-#define ANSIC 1
-
-int mtherr ( char *, int );
-
-/* Variable for error reporting. See mtherr.c. */
-extern int merror;
-/* === mconf.h - end === */
-/* === const.c - start === */
-/* const.c
+/* Common error handling routine
*
- * Globally declared constants
- *
- *
- *
* SYNOPSIS:
+ * char *fctnam;
+ * int code;
+ * int math_set_error();
*
- * extern double nameofconstant;
+ * math_set_error( fctnam, code );
*
- *
- *
- *
* DESCRIPTION:
+ * This routine may be called to report one of the following
+ * error conditions (in the include file mconf.h).
+ * Mnemonic Value Significance
+ * MATH_ERROR_DOMAIN 1 argument domain error
+ * MATH_ERROR_SING 2 function singularity
+ * MATH_ERROR_OVERFLOW 3 overflow range error
+ * MATH_ERROR_UNDERFLOW 4 underflow range error
+ * MATH_ERROR_TOTAL_LOSS 5 total loss of precision
+ * MATH_ERROR_PARTIAL_LOSS 6 partial loss of precision
*
- * This file contains a number of mathematical constants and
- * also some needed size parameters of the computer arithmetic.
- * The values are supplied as arrays of hexadecimal integers
- * for IEEE arithmetic; arrays of octal constants for DEC
- * arithmetic; and in a normal decimal scientific notation for
- * other machines. The particular notation used is determined
- * by a symbol (DEC, IBMPC, or UNK) defined in the include file
- * mconf.h.
+ * The default version of the file prints the function name,
+ * passed to it by the pointer fctnam, followed by the
+ * error condition. The display is directed to the standard
+ * output device. The routine then returns to the calling
+ * program. Users may wish to modify the program to abort by
+ * calling exit() under severe error conditions such as domain
+ * errors.
*
- * The default size parameters are as follows.
- *
- * For DEC and UNK modes:
- * MACHEP = 1.38777878078144567553E-17 2**-56
- * MAXLOG = 8.8029691931113054295988E1 log(2**127)
- * MINLOG = -8.872283911167299960540E1 log(2**-128)
- * MAXNUM = 1.701411834604692317316873e38 2**127
- *
- * For IEEE arithmetic (IBMPC):
- * MACHEP = 1.11022302462515654042E-16 2**-53
- * MAXLOG = 7.09782712893383996843E2 log(2**1024)
- * MINLOG = -7.08396418532264106224E2 log(2**-1022)
- * MAXNUM = 1.7976931348623158E308 2**1024
- *
- * The global symbols for mathematical constants are
- * PI = 3.14159265358979323846 pi
- * PIO2 = 1.57079632679489661923 pi/2
- * PIO4 = 7.85398163397448309616E-1 pi/4
- * SQRT2 = 1.41421356237309504880 sqrt(2)
- * SQRTH = 7.07106781186547524401E-1 sqrt(2)/2
- * LOG2E = 1.4426950408889634073599 1/log(2)
- * SQ2OPI = 7.9788456080286535587989E-1 sqrt( 2/pi )
- * LOGE2 = 6.93147180559945309417E-1 log(2)
- * LOGSQ2 = 3.46573590279972654709E-1 log(2)/2
- * THPIO4 = 2.35619449019234492885 3*pi/4
- * TWOOPI = 6.36619772367581343075535E-1 2/pi
- *
- * These lists are subject to change.
+ * Since all error conditions pass control to this function,
+ * the display may be easily changed, eliminated, or directed
+ * to an error logging device.
*/
-
-/* const.c */
+static int
+math_set_error (char *name, int code)
+{
+ /* Notice: the order of appearance of the following
+ * messages is bound to the error codes defined
+ * in mconf.h.
+ */
+ static const char *ermsg[7] = {
+ "unknown", /* error code 0 */
+ "domain", /* error code 1 */
+ "singularity", /* et seq. */
+ "overflow",
+ "underflow",
+ "total loss of precision",
+ "partial loss of precision"
+ };
+
+ /* Display string passed by calling program,
+ * which is supposed to be the name of the
+ * function in which the error occurred:
+ */
+ printf ("\n%s ", name); // FIXME
+
+ /* Set global error message word */
+ math_global_error = code;
+
+ /* Display error message defined
+ * by the code argument.
+ */
+ if( (code <= 0) || (code >= 7) )
+ code = 0;
+ printf( "%s error\n", ermsg[code] );
+
+ /* Return to calling
+ * program
+ */
+ return( 0 );
+}
-#if 1
-double MACHEP = 1.11022302462515654042E-16; /* 2**-53 */
-#else
-double MACHEP = 1.38777878078144567553E-17; /* 2**-56 */
-#endif
-double UFLOWTHRESH = 2.22507385850720138309E-308; /* 2**-1022 */
-#ifdef DENORMAL
-double MAXLOG = 7.09782712893383996732E2; /* log(MAXNUM) */
-/* double MINLOG = -7.44440071921381262314E2; */ /* log(2**-1074) */
-double MINLOG = -7.451332191019412076235E2; /* log(2**-1075) */
-#else
-double MAXLOG = 7.08396418532264106224E2; /* log 2**1022 */
-double MINLOG = -7.08396418532264106224E2; /* log 2**-1022 */
-#endif
-double MAXNUM = 1.79769313486231570815E308; /* 2**1024*(1-MACHEP) */
-double PI = 3.14159265358979323846; /* pi */
-double PIO2 = 1.57079632679489661923; /* pi/2 */
-double PIO4 = 7.85398163397448309616E-1; /* pi/4 */
-double SQRT2 = 1.41421356237309504880; /* sqrt(2) */
-double SQRTH = 7.07106781186547524401E-1; /* sqrt(2)/2 */
-double LOG2E = 1.4426950408889634073599; /* 1/log(2) */
-double SQ2OPI = 7.9788456080286535587989E-1; /* sqrt( 2/pi ) */
-double LOGE2 = 6.93147180559945309417E-1; /* log(2) */
-double LOGSQ2 = 3.46573590279972654709E-1; /* log(2)/2 */
-double THPIO4 = 2.35619449019234492885; /* 3*pi/4 */
-double TWOOPI = 6.36619772367581343075535E-1; /* 2/pi */
-#ifdef INFINITIES
-double INFINITY = 1.0/0.0; /* 99e999; */
-#else
-double INFINITY = 1.79769313486231570815E308; /* 2**1024*(1-MACHEP) */
-#endif
-#ifdef NANS
-double NAN = 1.0/0.0 - 1.0/0.0;
-#else
-double NAN = 0.0;
-#endif
-#ifdef MINUSZERO
-double NEGZERO = -0.0;
-#else
-double NEGZERO = 0.0;
-#endif
-
-
-/* === const.c - end === */
-/* === protos.h - start === */
-/*
- * This file was automatically generated by version 1.7 of cextract.
- * Manual editing not recommended.
- *
- * Created: Sun Jan 9 15:07:08 2000
- */
-
-extern double cabs ( cmplx *z );
-extern void cadd ( cmplx *a, cmplx *b, cmplx *c );
-extern double cay ( double q );
-extern void cdiv ( cmplx *a, cmplx *b, cmplx *c );
-extern void cmov ( void *a, void *b );
-extern void cmul ( cmplx *a, cmplx *b, cmplx *c );
-extern void cneg ( cmplx *a );
-extern void csqrt ( cmplx *z, cmplx *w );
-extern void csub ( cmplx *a, cmplx *b, cmplx *c );
-extern double ellie ( double phi, double m );
-extern double ellik ( double phi, double m );
-extern double ellpe ( double x );
-extern int ellpj ( double u, double m, double *sn, double *cn, double *dn, double *ph );
-extern double ellpk ( double x );
-extern int getnum ( char *line, double *val );
-extern int lampln ( void );
-extern int main ( void );
-extern int mtherr ( char *name, int code );
-extern double p1evl ( double x, double coef[], int N );
-extern double polevl ( double x, double coef[], int N );
-extern int quadf ( double x, double y, int pzflg );
-extern double response ( double f, double amp );
-extern int spln ( void );
-extern int zplna ( void );
-extern int zplnb ( void );
-extern int zplnc ( void );
-
-/* === protos.h - end === */
-/* === cmplx.c - start === */
-/* cmplx.c
- *
- * Complex number arithmetic
- *
- *
- *
- * SYNOPSIS:
- *
- * typedef struct {
- * double r; real part
- * double i; imaginary part
- * }cmplx;
- *
- * cmplx *a, *b, *c;
- *
- * cadd( a, b, c ); c = b + a
- * csub( a, b, c ); c = b - a
- * cmul( a, b, c ); c = b * a
- * cdiv( a, b, c ); c = b / a
- * cneg( c ); c = -c
- * cmov( b, c ); c = b
- *
- *
- *
- * DESCRIPTION:
- *
- * Addition:
- * c.r = b.r + a.r
- * c.i = b.i + a.i
- *
- * Subtraction:
- * c.r = b.r - a.r
- * c.i = b.i - a.i
- *
- * Multiplication:
- * c.r = b.r * a.r - b.i * a.i
- * c.i = b.r * a.i + b.i * a.r
- *
- * Division:
- * d = a.r * a.r + a.i * a.i
- * c.r = (b.r * a.r + b.i * a.i)/d
- * c.i = (b.i * a.r - b.r * a.i)/d
- * ACCURACY:
- *
- * In DEC arithmetic, the test (1/z) * z = 1 had peak relative
- * error 3.1e-17, rms 1.2e-17. The test (y/z) * (z/y) = 1 had
- * peak relative error 8.3e-17, rms 2.1e-17.
- *
- * Tests in the rectangle {-10,+10}:
- * Relative error:
- * arithmetic function # trials peak rms
- * DEC cadd 10000 1.4e-17 3.4e-18
- * IEEE cadd 100000 1.1e-16 2.7e-17
- * DEC csub 10000 1.4e-17 4.5e-18
- * IEEE csub 100000 1.1e-16 3.4e-17
- * DEC cmul 3000 2.3e-17 8.7e-18
- * IEEE cmul 100000 2.1e-16 6.9e-17
- * DEC cdiv 18000 4.9e-17 1.3e-17
- * IEEE cdiv 100000 3.7e-16 1.1e-16
- */
- /* cmplx.c
- * complex number arithmetic
- */
-
-
-void cdiv ( cmplx *, cmplx *, cmplx * );
-void cadd ( cmplx *, cmplx *, cmplx * );
-
-extern double MAXNUM, MACHEP, PI, PIO2, INFINITY, NAN;
-
-cmplx czero = {0.0, 0.0};
-extern cmplx czero;
-cmplx cone = {1.0, 0.0};
-extern cmplx cone;
-
-/* c = b + a */
-void cadd( a, b, c )
- register cmplx *a, *b;
- cmplx *c;
+/* --- complex number arithmetic --- */
+/* c = b + a */
+static void
+Cadd (const Complex *a, const Complex *b, Complex *c)
{
c->r = b->r + a->r;
c->i = b->i + a->i;
}
-
-/* c = b - a */
-void csub( a, b, c )
- register cmplx *a, *b;
- cmplx *c;
+/* c = b - a */
+static void
+Csub (const Complex *a, const Complex *b, Complex *c)
{
c->r = b->r - a->r;
c->i = b->i - a->i;
}
-/* c = b * a */
-void cmul( a, b, c )
- register cmplx *a, *b;
- cmplx *c;
+/* c = b * a */
+static void
+Cmul (const Complex *a, const Complex *b, Complex *c)
{
+ /* Multiplication:
+ * c.r = b.r * a.r - b.i * a.i
+ * c.i = b.r * a.i + b.i * a.r
+ */
double y;
y = b->r * a->r - b->i * a->i;
c->i = b->r * a->i + b->i * a->r;
c->r = y;
+ /* see Cdiv() for accuracy comments */
}
-
-
-/* c = b / a */
-void cdiv( a, b, c )
- register cmplx *a, *b;
- cmplx *c;
+/* c = b / a */
+static void
+Cdiv (const Complex *a, const Complex *b, Complex *c)
{
- double y, p, q, w;
+ /* Division:
+ * d = a.r * a.r + a.i * a.i
+ * c.r = (b.r * a.r + b.i * a.i)/d
+ * c.i = (b.i * a.r - b.r * a.i)/d
+ */
+ double y = a->r * a->r + a->i * a->i;
+ double p = b->r * a->r + b->i * a->i;
+ double q = b->i * a->r - b->r * a->i;
-
- y = a->r * a->r + a->i * a->i;
- p = b->r * a->r + b->i * a->i;
- q = b->i * a->r - b->r * a->i;
-
if( y < 1.0 )
{
- w = MAXNUM * y;
+ double w = MAXNUM * y;
if( (fabs(p) > w) || (fabs(q) > w) || (y == 0.0) )
{
c->r = MAXNUM;
c->i = MAXNUM;
- mtherr( "cdiv", OVERFLOW );
+ math_set_error( "Cdiv", MATH_ERROR_OVERFLOW );
return;
}
}
c->r = p/y;
c->i = q/y;
+ /* ACCURACY:
+ * In DEC arithmetic, the test (1/z) * z = 1 had peak relative
+ * error 3.1e-17, rms 1.2e-17. The test (y/z) * (z/y) = 1 had
+ * peak relative error 8.3e-17, rms 2.1e-17.
+ * Tests in the rectangle {-10,+10}:
+ * Relative error:
+ * arithmetic function # trials peak rms
+ * DEC Cadd 10000 1.4e-17 3.4e-18
+ * IEEE Cadd 100000 1.1e-16 2.7e-17
+ * DEC Csub 10000 1.4e-17 4.5e-18
+ * IEEE Csub 100000 1.1e-16 3.4e-17
+ * DEC Cmul 3000 2.3e-17 8.7e-18
+ * IEEE Cmul 100000 2.1e-16 6.9e-17
+ * DEC Cdiv 18000 4.9e-17 1.3e-17
+ * IEEE Cdiv 100000 3.7e-16 1.1e-16
+ */
}
-
-/* b = a
- Caution, a `short' is assumed to be 16 bits wide. */
-
-void cmov( a, b )
- void *a, *b;
+/* b = a */
+static void
+Cmov (const Complex *a, Complex *b)
{
- register short *pa, *pb;
- int i;
-
- pa = (short *) a;
- pb = (short *) b;
- i = 8;
- do
- *pb++ = *pa++;
- while( --i );
+ *b = *a;
}
-
-void cneg( a )
- register cmplx *a;
-{
-
- a->r = -a->r;
- a->i = -a->i;
-}
-
-/* cabs()
+/* Cabs() - Complex absolute value
*
- * Complex absolute value
- *
- *
- *
* SYNOPSIS:
- *
- * double cabs();
- * cmplx z;
+ * double Cabs();
+ * Complex z;
* double a;
*
- * a = cabs( &z );
+ * a = Cabs( &z );
*
- *
- *
* DESCRIPTION:
- *
- *
* If z = x + iy
- *
* then
- *
* a = sqrt( x**2 + y**2 ).
- *
* Overflow and underflow are avoided by testing the magnitudes
* of x and y before squaring. If either is outside half of
* the floating point full scale range, both are rescaled.
*
- *
* ACCURACY:
- *
* Relative error:
* arithmetic domain # trials peak rms
* DEC -30,+30 30000 3.2e-17 9.2e-18
* IEEE -10,+10 100000 2.7e-16 6.9e-17
*/
-
-#define PREC 27
-#define MAXEXP 1024
-#define MINEXP -1077
-
-
-double cabs( z )
- register cmplx *z;
+static double
+Cabs (const Complex *z)
{
+ /* exponent thresholds for IEEE doubles */
+ const double PREC = 27;
+ const double MAXEXP = 1024;
+ const double MINEXP = -1077;
+
double x, y, b, re, im;
int ex, ey, e;
-#ifdef INFINITIES
- /* Note, cabs(INFINITY,NAN) = INFINITY. */
+ /* Note, Cabs(INFINITY,NAN) = INFINITY. */
if( z->r == INFINITY || z->i == INFINITY
|| z->r == -INFINITY || z->i == -INFINITY )
return( INFINITY );
-#endif
-#ifdef NANS
if( isnan(z->r) )
return(z->r);
if( isnan(z->i) )
return(z->i);
-#endif
re = fabs( z->r );
im = fabs( z->i );
@@ -717,7 +402,7 @@
/* Check it for overflow and underflow. */
if( ey > MAXEXP )
{
- mtherr( "cabs", OVERFLOW );
+ math_set_error( "Cabs", MATH_ERROR_OVERFLOW );
return( INFINITY );
}
if( ey < MINEXP )
@@ -728,57 +413,38 @@
return( b );
}
-/* csqrt()
+/* Csqrt() - Complex square root
*
- * Complex square root
- *
- *
- *
* SYNOPSIS:
+ * void Csqrt();
+ * Complex z, w;
+ * Csqrt( &z, &w );
*
- * void csqrt();
- * cmplx z, w;
- *
- * csqrt( &z, &w );
- *
- *
- *
* DESCRIPTION:
- *
- *
* If z = x + iy, r = |z|, then
- *
* 1/2
* Im w = [ (r - x)/2 ] ,
- *
* Re w = y / 2 Im w.
*
- *
* Note that -w is also a square root of z. The root chosen
* is always in the upper half plane.
- *
* Because of the potential for cancellation error in r - x,
* the result is sharpened by doing a Heron iteration
* (see sqrt.c) in complex arithmetic.
*
- *
- *
* ACCURACY:
- *
* Relative error:
* arithmetic domain # trials peak rms
* DEC -10,+10 25000 3.2e-17 9.6e-18
* IEEE -10,+10 100000 3.2e-16 7.7e-17
- *
* 2
- * Also tested by csqrt( z ) = z, and tested by arguments
+ * Also tested by Csqrt( z ) = z, and tested by arguments
* close to the real axis.
*/
-
-void csqrt( z, w )
- cmplx *z, *w;
+static void
+Csqrt (const Complex *z, Complex *w)
{
- cmplx q, s;
+ Complex q, s;
double x, y, r, t;
x = z->r;
@@ -823,7 +489,7 @@
}
else
{
- r = cabs(z);
+ r = Cabs(z);
t = 0.5*(r - x);
}
@@ -831,45 +497,21 @@
q.i = r;
q.r = y/(2.0*r);
/* Heron iteration in complex arithmetic */
- cdiv( &q, z, &s );
- cadd( &q, &s, w );
+ Cdiv( &q, z, &s );
+ Cadd( &q, &s, w );
w->r *= 0.5;
w->i *= 0.5;
}
-
-double hypot( x, y )
- double x, y;
-{
- cmplx z;
-
- z.r = x;
- z.i = y;
- return( cabs(&z) );
-}
-
-/* === cmplx.c - end === */
-/* === ellik.c - start === */
-/* ellik.c
+/* --- elliptic functions --- */
+/* ellik.c - Incomplete elliptic integral of the first kind
*
- * Incomplete elliptic integral of the first kind
- *
- *
- *
* SYNOPSIS:
- *
* double phi, m, y, ellik();
- *
* y = ellik( phi, m );
*
- *
- *
* DESCRIPTION:
- *
* Approximates the integral
- *
- *
- *
* phi
* -
* | |
@@ -879,32 +521,17 @@
* | | sqrt( 1 - m sin t )
* -
* 0
- *
* of amplitude phi and modulus m, using the arithmetic -
* geometric mean algorithm.
*
- *
- *
- *
* ACCURACY:
- *
* Tested at random points with m in [0, 1] and phi as indicated.
- *
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -10,10 200000 7.4e-16 1.0e-16
- *
- *
*/
-
-/* Incomplete elliptic integral of first kind */
-
-extern double ellpk ( double );
-double ellik ( double, double );
-extern double PI, PIO2, MACHEP, MAXNUM;
-
-double ellik( phi, m )
- double phi, m;
+static double
+ellik (double phi, double m)
{
double a, b, c, e, temp, t, K;
int d, mod, sign, npio2;
@@ -916,7 +543,7 @@
{
if( fabs(phi) >= PIO2 )
{
- mtherr( "ellik", SING );
+ math_set_error( "ellik", MATH_ERROR_SING );
return( MAXNUM );
}
return( log( tan( (PIO2 + phi)/2.0 ) ) );
@@ -981,126 +608,14 @@
return( temp );
}
-/* === ellik.c - end === */
-/* === ellpe.c - start === */
-/* ellpe.c
+/* ellpj - Jacobian Elliptic Functions
*
- * Complete elliptic integral of the second kind
- *
- *
- *
* SYNOPSIS:
- *
- * double m1, y, ellpe();
- *
- * y = ellpe( m1 );
- *
- *
- *
- * DESCRIPTION:
- *
- * Approximates the integral
- *
- *
- * pi/2
- * -
- * | | 2
- * E(m) = | sqrt( 1 - m sin t ) dt
- * | |
- * -
- * 0
- *
- * Where m = 1 - m1, using the approximation
- *
- * P(x) - x log x Q(x).
- *
- * Though there are no singularities, the argument m1 is used
- * rather than m for compatibility with ellpk().
- *
- * E(1) = 1; E(0) = pi/2.
- *
- *
- * ACCURACY:
- *
- * Relative error:
- * arithmetic domain # trials peak rms
- * DEC 0, 1 13000 3.1e-17 9.4e-18
- * IEEE 0, 1 10000 2.1e-16 7.3e-17
- *
- *
- * ERROR MESSAGES:
- *
- * message condition value returned
- * ellpe domain x<0, x>1 0.0
- *
- */
-/* ellpe.c */
-
-/* Elliptic integral of second kind */
-
-static double P_ellpe[] = {
- 1.53552577301013293365E-4,
- 2.50888492163602060990E-3,
- 8.68786816565889628429E-3,
- 1.07350949056076193403E-2,
- 7.77395492516787092951E-3,
- 7.58395289413514708519E-3,
- 1.15688436810574127319E-2,
- 2.18317996015557253103E-2,
- 5.68051945617860553470E-2,
- 4.43147180560990850618E-1,
- 1.00000000000000000299E0
-};
-static double Q_ellpe[] = {
- 3.27954898576485872656E-5,
- 1.00962792679356715133E-3,
- 6.50609489976927491433E-3,
- 1.68862163993311317300E-2,
- 2.61769742454493659583E-2,
- 3.34833904888224918614E-2,
- 4.27180926518931511717E-2,
- 5.85936634471101055642E-2,
- 9.37499997197644278445E-2,
- 2.49999999999888314361E-1
-};
-
-
-extern double polevl ( double, double[], int );
-extern double log ( double );
-
-double ellpe(x)
- double x;
-{
-
- if( (x <= 0.0) || (x > 1.0) )
- {
- if( x == 0.0 )
- return( 1.0 );
- mtherr( "ellpe", DOMAIN );
- return( 0.0 );
- }
- return( polevl(x,P_ellpe,10) - log(x) * (x * polevl(x,Q_ellpe,9)) );
-}
-/* === ellpe.c - end === */
-/* === ellpj.c - start === */
-/* ellpj.c
- *
- * Jacobian Elliptic Functions
- *
- *
- *
- * SYNOPSIS:
- *
* double u, m, sn, cn, dn, phi;
* int ellpj();
- *
* ellpj( u, m, _&sn, _&cn, _&dn, _&phi );
*
- *
- *
* DESCRIPTION:
- *
- *
* Evaluates the Jacobian elliptic functions sn(u|m), cn(u|m),
* and dn(u|m) of parameter m between 0 and 1, and real
* argument u.
@@ -1119,10 +634,8 @@
* only for phi < pi/2.
*
* ACCURACY:
- *
* Tested at random points with u between 0 and 10, m between
* 0 and 1.
- *
* Absolute error (* = relative error):
* arithmetic function # trials peak rms
* DEC sn 1800 4.5e-16 8.7e-17
@@ -1130,33 +643,23 @@
* IEEE sn 50000 4.1e-15 4.6e-16
* IEEE cn 40000 3.6e-15 4.4e-16
* IEEE dn 10000 1.3e-12 1.8e-14
- *
- * Peak error observed in consistency check using addition
+ * Peak error observed in consistency check using addition
* theorem for sn(u+v) was 4e-16 (absolute). Also tested by
* the above relation to the incomplete elliptic integral.
* Accuracy deteriorates when u is large.
- *
*/
-
-/* ellpj.c */
-
-
-extern double PIO2, MACHEP;
-
-int ellpj( u, m, sn, cn, dn, ph )
- double u, m;
- double *sn, *cn, *dn, *ph;
+static int
+ellpj (double u, double m,
+ double *sn, double *cn, double *dn,
+ double *ph)
{
double ai, b, phi, t, twon;
double a[9], c[9];
int i;
-
-
/* Check for special cases */
-
if( m < 0.0 || m > 1.0 )
{
- mtherr( "ellpj", DOMAIN );
+ math_set_error( "ellpj", MATH_ERROR_DOMAIN );
*sn = 0.0;
*cn = 0.0;
*ph = 0.0;
@@ -1189,8 +692,6 @@
*dn = phi + ai * (twon + u);
return(0);
}
-
-
/* A. G. M. scale */
a[0] = 1.0;
b = sqrt(1.0 - m);
@@ -1202,7 +703,7 @@
{
if( i > 7 )
{
- mtherr( "ellpj", OVERFLOW );
+ math_set_error( "ellpj", MATH_ERROR_OVERFLOW );
goto done;
}
ai = a[i];
@@ -1215,7 +716,6 @@
}
done:
-
/* backward recurrence */
phi = twon * a[i] * u;
do
@@ -1233,28 +733,15 @@
*ph = phi;
return(0);
}
-/* === ellpj.c - end === */
-/* === ellpk.c - start === */
-/* ellpk.c
+
+/* ellpk - Complete elliptic integral of the first kind
*
- * Complete elliptic integral of the first kind
- *
- *
- *
* SYNOPSIS:
- *
* double m1, y, ellpk();
- *
* y = ellpk( m1 );
*
- *
- *
* DESCRIPTION:
- *
* Approximates the integral
- *
- *
- *
* pi/2
* -
* | |
@@ -1264,78 +751,57 @@
* | | sqrt( 1 - m sin t )
* -
* 0
- *
* where m = 1 - m1, using the approximation
- *
* P(x) - log x Q(x).
*
* The argument m1 is used rather than m so that the logarithmic
* singularity at m = 1 will be shifted to the origin; this
* preserves maximum accuracy.
- *
* K(0) = pi/2.
*
* ACCURACY:
- *
* Relative error:
* arithmetic domain # trials peak rms
* DEC 0,1 16000 3.5e-17 1.1e-17
* IEEE 0,1 30000 2.5e-16 6.8e-17
- *
* ERROR MESSAGES:
- *
* message condition value returned
* ellpk domain x<0, x>1 0.0
- *
*/
-
-/* ellpk.c */
-
-
-
-static double P_ellpk[] =
+static double
+ellpk (double x)
{
- 1.37982864606273237150E-4,
- 2.28025724005875567385E-3,
- 7.97404013220415179367E-3,
- 9.85821379021226008714E-3,
- 6.87489687449949877925E-3,
- 6.18901033637687613229E-3,
- 8.79078273952743772254E-3,
- 1.49380448916805252718E-2,
- 3.08851465246711995998E-2,
- 9.65735902811690126535E-2,
- 1.38629436111989062502E0
-};
+ static const double P_ellpk[] = {
+ 1.37982864606273237150E-4,
+ 2.28025724005875567385E-3,
+ 7.97404013220415179367E-3,
+ 9.85821379021226008714E-3,
+ 6.87489687449949877925E-3,
+ 6.18901033637687613229E-3,
+ 8.79078273952743772254E-3,
+ 1.49380448916805252718E-2,
+ 3.08851465246711995998E-2,
+ 9.65735902811690126535E-2,
+ 1.38629436111989062502E0
+ };
+ static const double Q_ellpk[] = {
+ 2.94078955048598507511E-5,
+ 9.14184723865917226571E-4,
+ 5.94058303753167793257E-3,
+ 1.54850516649762399335E-2,
+ 2.39089602715924892727E-2,
+ 3.01204715227604046988E-2,
+ 3.73774314173823228969E-2,
+ 4.88280347570998239232E-2,
+ 7.03124996963957469739E-2,
+ 1.24999999999870820058E-1,
+ 4.99999999999999999821E-1
+ };
+ static const double C1_ellpk = 1.3862943611198906188E0; /* log(4) */
-static double Q_ellpk[] =
-{
- 2.94078955048598507511E-5,
- 9.14184723865917226571E-4,
- 5.94058303753167793257E-3,
- 1.54850516649762399335E-2,
- 2.39089602715924892727E-2,
- 3.01204715227604046988E-2,
- 3.73774314173823228969E-2,
- 4.88280347570998239232E-2,
- 7.03124996963957469739E-2,
- 1.24999999999870820058E-1,
- 4.99999999999999999821E-1
-};
-static double C1 = 1.3862943611198906188E0; /* log(4) */
-
-extern double polevl ( double, double[], int );
-extern double p1evl ( double, double[], int );
-extern double log ( double );
-extern double MACHEP, MAXNUM;
-
-double ellpk(x)
- double x;
-{
-
if( (x < 0.0) || (x > 1.0) )
{
- mtherr( "ellpk", DOMAIN );
+ math_set_error( "ellpk", MATH_ERROR_DOMAIN );
return( 0.0 );
}
@@ -1347,168 +813,46 @@
{
if( x == 0.0 )
{
- mtherr( "ellpk", SING );
+ math_set_error( "ellpk", MATH_ERROR_SING );
return( MAXNUM );
}
else
{
- return( C1 - 0.5 * log(x) );
+ return( C1_ellpk - 0.5 * log(x) );
}
}
}
-/* === ellpk.c - end === */
-/* === mtherr.c - start === */
-/* mtherr.c
- *
- * Library common error handling routine
- *
- *
- *
- * SYNOPSIS:
- *
- * char *fctnam;
- * int code;
- * int mtherr();
- *
- * mtherr( fctnam, code );
- *
- *
- *
- * DESCRIPTION:
- *
- * This routine may be called to report one of the following
- * error conditions (in the include file mconf.h).
- *
- * Mnemonic Value Significance
- *
- * DOMAIN 1 argument domain error
- * SING 2 function singularity
- * OVERFLOW 3 overflow range error
- * UNDERFLOW 4 underflow range error
- * TLOSS 5 total loss of precision
- * PLOSS 6 partial loss of precision
- * EDOM 33 Unix domain error code
- * ERANGE 34 Unix range error code
- *
- * The default version of the file prints the function name,
- * passed to it by the pointer fctnam, followed by the
- * error condition. The display is directed to the standard
- * output device. The routine then returns to the calling
- * program. Users may wish to modify the program to abort by
- * calling exit() under severe error conditions such as domain
- * errors.
- *
- * Since all error conditions pass control to this function,
- * the display may be easily changed, eliminated, or directed
- * to an error logging device.
- *
- * SEE ALSO:
- *
- * mconf.h
- *
- */
-
-#include <stdio.h>
-
-int merror = 0;
-
-/* Notice: the order of appearance of the following
- * messages is bound to the error codes defined
- * in mconf.h.
- */
-static char *ermsg[7] = {
- "unknown", /* error code 0 */
- "domain", /* error code 1 */
- "singularity", /* et seq. */
- "overflow",
- "underflow",
- "total loss of precision",
- "partial loss of precision"
-};
-
-
-int mtherr( name, code )
- char *name;
- int code;
-{
-
- /* Display string passed by calling program,
- * which is supposed to be the name of the
- * function in which the error occurred:
- */
- printf( "\n%s ", name );
-
- /* Set global error message word */
- merror = code;
-
- /* Display error message defined
- * by the code argument.
- */
- if( (code <= 0) || (code >= 7) )
- code = 0;
- printf( "%s error\n", ermsg[code] );
-
- /* Return to calling
- * program
- */
- return( 0 );
-}
-/* === mtherr.c - end === */
-/* === polevl.c - start === */
-/* polevl.c
- * p1evl.c
+/* --- misc utilities --- */
+/* polevl - Evaluate polynomial
*
- * Evaluate polynomial
- *
- *
- *
* SYNOPSIS:
- *
* int N;
* double x, y, coef[N+1], polevl[];
- *
* y = polevl( x, coef, N );
*
- *
- *
* DESCRIPTION:
- *
* Evaluates polynomial of degree N:
- *
* 2 N
* y = C + C x + C x +...+ C x
* 0 1 2 N
*
* Coefficients are stored in reverse order:
- *
* coef[0] = C , ..., coef[N] = C .
* N 0
- *
- * The function p1evl() assumes that coef[N] = 1.0 and is
- * omitted from the array. Its calling arguments are
- * otherwise the same as polevl().
- *
- *
* SPEED:
- *
* In the interest of speed, there are no checks for out
* of bounds arithmetic. This routine is used by most of
* the functions in the library. Depending on available
* equipment features, the user may wish to rewrite the
* program in microcode or assembly language.
- *
*/
-
-
-double polevl( x, coef, N )
- double x;
- double coef[];
- int N;
+static double
+polevl (double x, const double coef[], int N)
{
double ans;
int i;
- double *p;
+ const double *p;
p = coef;
ans = *p++;
@@ -1521,82 +865,14 @@
return( ans );
}
-/* p1evl() */
-/* N
- * Evaluate polynomial when coefficient of x is 1.0.
- * Otherwise same as polevl.
- */
-
-double p1evl( x, coef, N )
- double x;
- double coef[];
- int N;
-{
- double ans;
- double *p;
- int i;
-
- p = coef;
- ans = x + *p++;
- i = N-1;
-
- do
- ans = ans * x + *p++;
- while( --i );
-
- return( ans );
-}
-/* === polevl.c - end === */
-/* === ellf.c - start === */
-/* ellf.c
- *
- * Read ellf.doc before attempting to compile this program.
- */
-
-
-#include <stdio.h>
-#include <stdlib.h>
-
-/* size of arrays: */
-#define ARRSIZ 300
-
-
-/* System configurations */
-
-
-extern double PI, PIO2, MACHEP, MAXNUM;
-
+#if 1
static double aa[ARRSIZ];
static double pp[ARRSIZ];
static double y[ARRSIZ];
static double zs[ARRSIZ];
-cmplx z[ARRSIZ];
-static double wr = 0.0;
-static double cbp = 0.0;
-static double wc = 0.0;
-static double rn = 8.0;
-static double c = 0.0;
-static double cgam = 0.0;
-static double scale = 0.0;
-double fs = 1.0e4; /* sampling frequency -- stw */
-static double dbr = 0.5;
-static double dbd = -40.0;
-static double f1 = 1.5e3;
-static double f2 = 2.0e3;
-static double f3 = 2.4e3;
-double dbfac = 0.0;
-static double a = 0.0;
-static double b = 0.0;
-static double q = 0.0;
-static double r = 0.0;
-static double u = 0.0;
-static double k = 0.0;
-static double m = 0.0;
+static Complex z[ARRSIZ];
static double Kk = 0.0;
-static double Kk1 = 0.0;
static double Kpk = 0.0;
-static double Kpk1 = 0.0;
-static double eps = 0.0;
static double rho = 0.0;
static double phi = 0.0;
static double sn = 0.0;
@@ -1606,308 +882,262 @@
static double cn1 = 0.0;
static double dn1 = 0.0;
static double phi1 = 0.0;
-static double m1 = 0.0;
-static double m1p = 0.0;
-static double cang = 0.0;
-static double sang = 0.0;
-static double bw = 0.0;
-static double ang = 0.0;
-double fnyq = 0.0; /* nyquist frequency -- stw */
-static double ai = 0.0;
static double pn = 0.0;
static double an = 0.0;
static double gam = 0.0;
-static double cng = 0.0;
-double gain = 0.0;
-static int lr = 0;
-static int nt = 0;
-static int i = 0;
-static int j = 0;
+static double gain = 0.0;
static int jt = 0;
-static int nc = 0;
static int ii = 0;
static int ir = 0;
-int zord = 0;
+static int zord = 0;
static int icnt = 0;
-static int mh = 0;
-static int jj = 0;
-static int jh = 0;
-static int jl = 0;
-static int n = 8;
-static int np = 0;
-static int nz = 0;
-static int type = 1;
-static int kind = 1;
+#endif
-static char wkind[] =
-{"Filter kind:\n1 Butterworth\n2 Chebyshev\n3 Elliptic\n"};
+/* --- prototypes --- */
+static int find_elliptic_locations_in_lambda_plane (const BseIIRFilterRequirements *ifr,
+ DesignState *ds);
+static int find_s_plane_poles_and_zeros (const BseIIRFilterRequirements *ifr,
+ DesignState *ds);
+static int convert_s_plane_to_z_plane (const BseIIRFilterRequirements *ifr,
+ DesignState *ds);
+static double jacobi_theta_by_nome (double q);
+static int z_plane_zeros_poles_to_numerator_denomerator (const BseIIRFilterRequirements *ifr,
+ DesignState *ds);
+static int gainscale_and_print_deno_nume_zeros2_poles2 (const BseIIRFilterRequirements *ifr,
+ DesignState *ds);
+static int print_quadratic_factors (const BseIIRFilterRequirements *ifr,
+ DesignState *ds,
+ double x, double y, int pzflg);
+static void print_filter_table (const BseIIRFilterRequirements *ifr,
+ DesignState *ds);
+static double response (const BseIIRFilterRequirements *ifr,
+ DesignState *ds,
+ double f, double amp);
-static char salut[] =
-{"Filter shape:\n1 low pass\n2 band pass\n3 high pass\n4 band stop\n"};
-
-extern double cabs ( cmplx *z );
-extern void cadd ( cmplx *a, cmplx *b, cmplx *c );
-extern void cdiv ( cmplx *a, cmplx *b, cmplx *c );
-extern void cmov ( void *a, void *b );
-extern void cmul ( cmplx *a, cmplx *b, cmplx *c );
-extern void cneg ( cmplx *a );
-extern void csqrt ( cmplx *z, cmplx *w );
-extern void csub ( cmplx *a, cmplx *b, cmplx *c );
-extern double ellie ( double phi, double m );
-extern double ellik ( double phi, double m );
-extern double ellpe ( double x );
-extern int ellpj ( double, double, double *, double *, double *, double * );
-extern double ellpk ( double x );
-int getnum ( char *line, double *val );
-double cay ( double q );
-int lampln ( void );
-int spln ( void );
-void print_filter_table (void);
-int zplna ( void );
-int zplnb ( void );
-int zplnc ( void );
-int quadf ( double, double, int );
-double response ( double, double );
-
static void
-print_z_fraction_before_zplnc (void) /* must be called *before* zplnc() */
+print_z_fraction_before_zplnc (const BseIIRFilterRequirements *ifr,
+ DesignState *ds) /* must be called *before* zplnc() */
{
double zgain;
- if (kind != 3 && pn == 0)
+ if (ifr->kind != 3 && pn == 0)
zgain = 1.0;
else
- zgain = an / (pn * scale);
+ zgain = an / (pn * ds->gain_scale);
VERBOSE ("# constant mygain factor %23.13E\n", zgain); // BSE info
VERBOSE ("# z plane Denominator Numerator\n" ); // BSE info
int j;
for (j = 0; j <= zord; j++)
- VERBOSE ("%2d %17.9E %17.9E #BSEmgc\n", j, aa[j], pp[j] * zgain); // BSE info
+ VERBOSE ("%2d %17.9E %17.9E\n", j, aa[j], pp[j] * zgain); // BSE info
}
-int main()
+/* --- main IIR filter design function --- */
+static const char*
+iir_filter_design (const BseIIRFilterRequirements *ifr,
+ DesignState *ds)
{
- char str[80];
+ double passband_edge1 = ifr->passband_edge;
+ double passband_edge0 = ifr->passband_edge2;
- dbfac = 10.0/log(10.0);
-
- top:
-
- printf( "%s ? ", wkind ); /* ask for filter kind */
- if (!fgets (str, sizeof (str), stdin))
- exit (0);
- sscanf( str, "%d", &kind );
- printf( "%d\n", kind );
- if( (kind <= 0) || (kind > 3) )
- exit(0);
-
- printf( "%s ? ", salut ); /* ask for filter type */
- if (!fgets (str, sizeof (str), stdin))
- exit (0);
- sscanf( str, "%d", &type );
- printf( "%d\n", type );
- if( (type <= 0) || (type > 4) )
- exit(0);
-
- getnum( "Order of filter", &rn ); /* see below for getnum() */
- n = rn;
- if( n <= 0 )
+ if (ifr->kind <= 0 || ifr->kind > 3)
+ return "unknown kind";
+ if (ifr->type <= 0 || ifr->type > 4)
+ return "unknown type";
+ if (ifr->order <= 0)
+ return "order too small";
+
+ if (ifr->kind > 1) /* not Butterworth */
{
- specerr:
- printf( "? Specification error\n" );
- goto top;
- }
- rn = n; /* ensure it is an integer */
- if( kind > 1 ) /* not Butterworth */
- {
- getnum( "Passband ripple, db", &dbr );
- if( dbr <= 0.0 )
- goto specerr;
- if( kind == 2 )
+ if (ifr->passband_ripple_db <= 0.0)
+ return "passband_ripple_db too small";
+ if (ifr->kind == 2)
{
- /* For Chebyshev filter, ripples go from 1.0 to 1/sqrt(1+eps^2) */
- phi = exp( 0.5*dbr/dbfac );
-
- if( (n & 1) == 0 )
- scale = phi;
+ /* For Chebyshev filter, ripples go from 1.0 to 1/sqrt(1+ds->ripple_epsilon^2) */
+ phi = exp( 0.5*ifr->passband_ripple_db / DECIBELL_FACTOR );
+
+ if( (ifr->order & 1) == 0 )
+ ds->gain_scale = phi;
else
- scale = 1.0;
+ ds->gain_scale = 1.0;
}
else
{ /* elliptic */
- eps = exp( dbr/dbfac );
- scale = 1.0;
- if( (n & 1) == 0 )
- scale = sqrt( eps );
- eps = sqrt( eps - 1.0 );
+ ds->ripple_epsilon = exp( ifr->passband_ripple_db / DECIBELL_FACTOR );
+ ds->gain_scale = 1.0;
+ if( (ifr->order & 1) == 0 )
+ ds->gain_scale = sqrt( ds->ripple_epsilon );
+ ds->ripple_epsilon = sqrt( ds->ripple_epsilon - 1.0 );
}
}
- getnum( "Sampling frequency", &fs );
- if( fs <= 0.0 )
- goto specerr;
+ if (ifr->sampling_frequency <= 0.0)
+ return "sampling_frequency too small";
- fnyq = 0.5 * fs;
+ ds->nyquist_frequency = 0.5 * ifr->sampling_frequency;
- getnum( "Passband edge", &f2 );
- if( (f2 <= 0.0) || (f2 >= fnyq) )
- goto specerr;
+ if (passband_edge1 <= 0.0)
+ return "passband_edge1 too small";
+ if (passband_edge1 >= ds->nyquist_frequency)
+ return "passband_edge1 too high";
- if( (type & 1) == 0 )
+ if( (ifr->type & 1) == 0 )
{
- getnum( "Other passband edge", &f1 );
- if( (f1 <= 0.0) || (f1 >= fnyq) )
- goto specerr;
+ if (passband_edge0 <= 0.0)
+ return "passband_edge too small";
+ if (passband_edge0 >= ds->nyquist_frequency)
+ return "passband_edge too high";
}
else
{
- f1 = 0.0;
+ passband_edge0 = 0.0;
}
-
- if( f2 < f1 )
+ if (passband_edge1 < passband_edge0)
{
- a = f2;
- f2 = f1;
- f1 = a;
+ double tmp = passband_edge1;
+ passband_edge1 = passband_edge0;
+ passband_edge0 = tmp;
}
- if( type == 3 ) /* high pass */
+ double high_edge, band_width;
+ if (ifr->type == 3) /* high pass */
{
- bw = f2;
- a = fnyq;
+ band_width = passband_edge1;
+ high_edge = ds->nyquist_frequency;
}
else
{
- bw = f2 - f1;
- a = f2;
+ band_width = passband_edge1 - passband_edge0;
+ high_edge = passband_edge1;
}
/* Frequency correspondence for bilinear transformation
*
* Wanalog = tan( 2 pi Fdigital T / 2 )
*
- * where T = 1/fs
+ * where T = 1/ifr->sampling_frequency
*/
- ang = bw * PI / fs;
- cang = cos( ang );
- c = sin(ang) / cang; /* Wanalog */
- if( kind != 3 )
+ double ang = band_width * PI / ifr->sampling_frequency; /* angle frequency */
+ double sang;
+ double cang = cos (ang);
+ ds->tan_angle_frequency = sin (ang) / cang; /* Wanalog */
+ if( ifr->kind != 3 )
{
- wc = c;
- /*printf( "cos( 1/2 (Whigh-Wlow) T ) = %.5e, wc = %.5e\n", cang, wc );*/
+ ds->wc = ds->tan_angle_frequency;
+ /*printf( "cos( 1/2 (Whigh-Wlow) T ) = %.5e, wc = %.5e\n", cang, ds->wc );*/
}
- if( kind == 3 )
+ if( ifr->kind == 3 )
{ /* elliptic */
- cgam = cos( (a+f1) * PI / fs ) / cang;
- getnum( "Stop band edge or -(db down)", &dbd );
- if( dbd > 0.0 )
- f3 = dbd;
- else
+ double tmp_cgam = cos( (high_edge+passband_edge0) * PI / ifr->sampling_frequency ) / cang;
+ ds->cgam = tmp_cgam;
+ if (ifr->stopband_edge > 0.0)
+ ds->stopband_edge = ifr->stopband_edge;
+ else if (ifr->stopband_db >= 0.0)
+ return "need stopband_db or stopband_edge";
+ else /* stopband_db < 0.0 */
{ /* calculate band edge from db down */
- a = exp( -dbd/dbfac );
- m1 = eps/sqrt( a - 1.0 );
+ double a = exp (-ifr->stopband_db / DECIBELL_FACTOR);
+ double m1 = ds->ripple_epsilon/sqrt( a - 1.0 );
m1 *= m1;
- m1p = 1.0 - m1;
- Kk1 = ellpk( m1p );
- Kpk1 = ellpk( m1 );
- q = exp( -PI * Kpk1 / (rn * Kk1) );
- k = cay(q);
- if( type >= 3 )
- wr = k;
+ double m1p = 1.0 - m1;
+ double Kk1 = ellpk( m1p );
+ double Kpk1 = ellpk( m1 );
+ double q = exp( -PI * Kpk1 / (ifr->order * Kk1) );
+ fixme2local_1 = jacobi_theta_by_nome (q);
+ if( ifr->type >= 3 )
+ ds->wr = fixme2local_1;
else
- wr = 1.0/k;
- if( type & 1 )
+ ds->wr = 1.0 / fixme2local_1;
+ if( ifr->type & 1 )
{
- f3 = atan( c * wr ) * fs / PI;
+ ds->stopband_edge = atan( ds->tan_angle_frequency * ds->wr ) * ifr->sampling_frequency / PI;
}
else
{
- a = c * wr;
- a *= a;
- b = a * (1.0 - cgam * cgam) + a * a;
- b = (cgam + sqrt(b))/(1.0 + a);
- f3 = (PI/2.0 - asin(b)) * fs / (2.0*PI);
+ // FIXME: using tmp_cgam here increases precision
+ fixme2local_a = ds->tan_angle_frequency * ds->wr;
+ fixme2local_a *= fixme2local_a;
+ double b = fixme2local_a * (1.0 - ds->cgam * ds->cgam) + fixme2local_a * fixme2local_a;
+ b = (ds->cgam + sqrt(b))/(1.0 + fixme2local_a);
+ ds->stopband_edge = (PI/2.0 - asin(b)) * ifr->sampling_frequency / (2.0*PI);
}
}
- switch( type )
+ switch( ifr->type )
{
case 1:
- if( f3 <= f2 )
- goto specerr;
+ if (ds->stopband_edge <= passband_edge1)
+ return "need stopband_edge > passband_edge";
break;
-
case 2:
- if( (f3 > f2) || (f3 < f1) )
- break;
- goto specerr;
-
+ if (ds->stopband_edge >= passband_edge0 && ds->stopband_edge <= passband_edge1)
+ return "need stopband_edge < passband_edge or stopband_edge > passband_edge2";
+ break;
case 3:
- if( f3 >= f2 )
- goto specerr;
+ if (ds->stopband_edge >= passband_edge1)
+ return "need stopband_edge < passband_edge";
break;
-
case 4:
- if( (f3 <= f1) || (f3 >= f2) )
- goto specerr;
+ if (ds->stopband_edge <= passband_edge0)
+ return "need stopband_edge > passband_edge2";
+ if (ds->stopband_edge >= passband_edge1)
+ return "need stopband_edge < passband_edge";
break;
}
- ang = f3 * PI / fs;
+ ang = ds->stopband_edge * PI / ifr->sampling_frequency;
cang = cos(ang);
sang = sin(ang);
-
- if( type & 1 )
+
+ if( ifr->type & 1 )
{
- wr = sang/(cang*c);
+ ds->wr = sang/(cang*ds->tan_angle_frequency);
}
else
{
- q = cang * cang - sang * sang;
+ double q = cang * cang - sang * sang;
sang = 2.0 * cang * sang;
cang = q;
- wr = (cgam - cang)/(sang * c);
+ ds->wr = (ds->cgam - cang)/(sang * ds->tan_angle_frequency);
}
-
- if( type >= 3 )
- wr = 1.0/wr;
- if( wr < 0.0 )
- wr = -wr;
+
+ if( ifr->type >= 3 )
+ ds->wr = 1.0/ds->wr;
+ if( ds->wr < 0.0 )
+ ds->wr = -ds->wr;
y[0] = 1.0;
- y[1] = wr;
- cbp = wr;
+ y[1] = ds->wr;
+ /* ds->chebyshev_band_cbp = ds->wr; */
- if( type >= 3 )
+ if( ifr->type >= 3 )
y[1] = 1.0/y[1];
- if( type & 1 )
+ if( ifr->type & 1 )
{
- for( i=1; i<=2; i++ )
+ int i;
+ for (i = 1; i <= 2; i++)
{
- aa[i] = atan( c * y[i-1] ) * fs / PI ;
+ aa[i] = atan( ds->tan_angle_frequency * y[i-1] ) * ifr->sampling_frequency / PI ;
}
printf( "pass band %.9E\n", aa[1] );
printf( "stop band %.9E\n", aa[2] );
}
else
{
- for( i=1; i<=2; i++ )
+ int i;
+ for (i = 1; i <= 2; i++)
{
- a = c * y[i-1];
- b = atan(a);
- q = sqrt( 1.0 + a * a - cgam * cgam );
-#ifdef ANSIC
- q = atan2( q, cgam );
-#else
- q = atan2( cgam, q );
-#endif
- aa[i] = (q + b) * fnyq / PI;
- pp[i] = (q - b) * fnyq / PI;
+ double a = ds->tan_angle_frequency * y[i-1];
+ double b = atan(a);
+ double q = sqrt( 1.0 + a * a - ds->cgam * ds->cgam );
+ q = atan2( q, ds->cgam );
+ aa[i] = (q + b) * ds->nyquist_frequency / PI;
+ pp[i] = (q - b) * ds->nyquist_frequency / PI;
}
printf( "pass band %.9E %.9E\n", pp[1], aa[1] );
printf( "stop band %.9E %.9E\n", pp[2], aa[2] );
}
- lampln(); /* find locations in lambda plane */
- if( (2*n+2) > ARRSIZ )
+ ds->wc = 1.0;
+ find_elliptic_locations_in_lambda_plane (ifr, ds); /* find locations in lambda plane */
+ if( (2*ifr->order+2) > ARRSIZ )
goto toosml;
- }
+ } /* elliptic */
/* Transformation from low-pass to band-pass critical frequencies
*
@@ -1923,135 +1153,136 @@
* sin( Wdigital T )
*/
- if( kind == 2 )
+ if( ifr->kind == 2 )
{ /* Chebyshev */
- a = PI * (a+f1) / fs ;
- cgam = cos(a) / cang;
- a = 2.0 * PI * f2 / fs;
- cbp = (cgam - cos(a))/sin(a);
+ double a = PI * (high_edge + passband_edge0) / ifr->sampling_frequency ;
+ ds->cgam = cos(a) / cang;
+ a = 2.0 * PI * passband_edge1 / ifr->sampling_frequency;
+ ds->chebyshev_band_cbp = (ds->cgam - cos (a)) / sin (a);
}
- if( kind == 1 )
+ if( ifr->kind == 1 )
{ /* Butterworth */
- a = PI * (a+f1) / fs ;
- cgam = cos(a) / cang;
- a = 2.0 * PI * f2 / fs;
- cbp = (cgam - cos(a))/sin(a);
- scale = 1.0;
+ double a = PI * (high_edge + passband_edge0) / ifr->sampling_frequency ;
+ ds->cgam = cos(a) / cang;
+ a = 2.0 * PI * passband_edge1 / ifr->sampling_frequency;
+ /* ds->chebyshev_band_cbp = (ds->cgam - cos (a)) / sin (a); */
+ ds->gain_scale = 1.0;
}
- EVERBOSE ("State: gain_scale=%.20g ripple_epsilon=%.20g nyquist_frequency=%.20g " // BSE info
+ EVERBOSE ("State: gain_scale=%.20g ripple_epsilon=%.20g nyquist_frequency=%.20g " // BSE info
"tan_angle_frequency=%.20g stopband_edge=%.20g wc=%.20g wr=%.20g cgam=%.20g\n",
- scale, eps, fnyq,
- c, f3, wc, wr, cgam);
+ ds->gain_scale, ds->ripple_epsilon, ds->nyquist_frequency,
+ ds->tan_angle_frequency, ds->stopband_edge, ds->wc, ds->wr, ds->cgam);
+
+ find_s_plane_poles_and_zeros (ifr, ds); /* find s plane poles and zeros */
- spln(); /* find s plane poles and zeros */
-
- if( ((type & 1) == 0) && ((4*n+2) > ARRSIZ) )
+ if( ((ifr->type & 1) == 0) && ((4*ifr->order+2) > ARRSIZ) )
goto toosml;
- zplna(); /* convert s plane to z plane */
- zplnb();
- EVERBOSE ("an=%.20g pn=%.20g scale=%.20g\n", an, pn, scale); // BSE info
- print_z_fraction_before_zplnc ();
- zplnc();
- print_filter_table(); /* tabulate transfer function */
- goto top;
-
+ convert_s_plane_to_z_plane (ifr, ds); /* convert s plane to z plane */
+ // volatile_sink ("x");
+ z_plane_zeros_poles_to_numerator_denomerator (ifr, ds);
+ EVERBOSE ("an=%.20g pn=%.20g scale=%.20g\n", an, pn, ds->gain_scale); // BSE info
+ print_z_fraction_before_zplnc (ifr, ds);
+ gainscale_and_print_deno_nume_zeros2_poles2 (ifr, ds);
+ print_filter_table (ifr, ds); /* tabulate transfer function */
+ return NULL;
+
toosml:
- printf( "Cannot continue, storage arrays too small\n" );
- goto top;
+ return "storage arrays too small";
}
-int lampln()
+static int
+find_elliptic_locations_in_lambda_plane (const BseIIRFilterRequirements *ifr,
+ DesignState *ds)
{
-
- wc = 1.0;
- k = wc/wr;
- m = k * k;
- Kk = ellpk( 1.0 - m );
- Kpk = ellpk( m );
- EVERBOSE ("check: k=%.20g m=%.20g Kk=%.20g Kpk=%.20g\n", k, m, Kk, Kpk); // BSE info
- q = exp( -PI * rn * Kpk / Kk ); /* the nome of k1 */
- m1 = cay(q); /* see below */
- /* Note m1 = eps / sqrt( A*A - 1.0 ) */
- a = eps/m1;
+ fixme2local_k = ds->wc / ds->wr;
+ fixme2local_m = fixme2local_k * fixme2local_k;
+ Kk = ellpk (1.0 - fixme2local_m);
+ Kpk = ellpk (fixme2local_m);
+ EVERBOSE ("check: k=%.20g m=%.20g Kk=%.20g Kpk=%.20g\n", fixme2local_k, fixme2local_m, Kk, Kpk); // BSE info
+ double q = exp( -PI * ifr->order * Kpk / Kk ); /* the nome of k1 */
+ double m1 = jacobi_theta_by_nome (q); /* see below */
+ /* Note m1 = ds->ripple_epsilon / sqrt( A*A - 1.0 ) */
+ double a = ds->ripple_epsilon/m1;
a = a * a + 1;
a = 10.0 * log(a) / log(10.0);
printf( "dbdown %.9E\n", a );
- a = 180.0 * asin( k ) / PI;
- b = 1.0/(1.0 + eps*eps);
+ a = 180.0 * asin (fixme2local_k) / PI;
+ double b = 1.0/(1.0 + ds->ripple_epsilon*ds->ripple_epsilon);
b = sqrt( 1.0 - b );
printf( "theta %.9E, rho %.9E\n", a, b );
m1 *= m1;
- m1p = 1.0 - m1;
- Kk1 = ellpk( m1p );
- Kpk1 = ellpk( m1 );
- r = Kpk1 * Kk / (Kk1 * Kpk);
+ double m1p = 1.0 - m1;
+ double Kk1 = ellpk( m1p );
+ double Kpk1 = ellpk( m1 );
+ double r = Kpk1 * Kk / (Kk1 * Kpk);
printf( "consistency check: n= %.14E\n", r );
EVERBOSE ("consistency check: r=%.20g Kpk1=%.20g Kk1=%.20g m1=%.20g m1p=%.20g\n", r, Kpk1, Kk1, m1, m1p); // BSE info
/* -1
- * sn j/eps\m = j ellik( atan(1/eps), m )
+ * sn j/ds->ripple_epsilon\m = j ellik( atan(1/ds->ripple_epsilon), m )
*/
- b = 1.0/eps;
+ b = 1.0/ds->ripple_epsilon;
phi = atan( b );
- u = ellik( phi, m1p );
- EVERBOSE ("phi=%.20g m=%.20g u=%.20g\n", phi, m1p, u);
+ fixme2local_2 = ellik (phi, m1p);
+ EVERBOSE ("phi=%.20g m=%.20g u=%.20g\n", phi, m1p, fixme2local_2);
/* consistency check on inverse sn */
- ellpj( u, m1p, &sn, &cn, &dn, &phi );
+ ellpj (fixme2local_2, m1p, &sn, &cn, &dn, &phi);
a = sn/cn;
EVERBOSE ("consistency check: sn/cn = %.20g = %.20g = 1/ripple\n", a, b);
- u = u * Kk / (rn * Kk1); /* or, u = u * Kpk / Kpk1 */
+ ds->elliptic_k = fixme2local_k;
+ ds->elliptic_u = fixme2local_2 * Kk / (ifr->order * Kk1); /* or, u = u * Kpk / Kpk1 */
+ ds->elliptic_m = fixme2local_m;
return 0;
}
-
-
-
/* calculate s plane poles and zeros, normalized to wc = 1 */
-int spln()
+static int
+find_s_plane_poles_and_zeros (const BseIIRFilterRequirements *ifr,
+ DesignState *ds)
{
- for( i=0; i<ARRSIZ; i++ )
+ int i, j;
+ for (i = 0; i < ARRSIZ; i++)
zs[i] = 0.0;
- np = (n+1)/2;
- nz = 0;
- if( kind == 1 )
+ ds->n_poles = (ifr->order + 1) / 2;
+ ds->n_zeros = 0;
+ if (ifr->kind == 1)
{
- /* Butterworth poles equally spaced around the unit circle
- */
- if( n & 1 )
+ double m;
+ /* Butterworth poles equally spaced around the unit circle */
+ if (ifr->order & 1)
m = 0.0;
else
- m = PI / (2.0*n);
- for( i=0; i<np; i++ )
+ m = PI / (2.0 * ifr->order);
+ for (i = 0; i < ds->n_poles; i++)
{ /* poles */
- lr = i + i;
+ int lr = i + i;
zs[lr] = -cos(m);
zs[lr+1] = sin(m);
- m += PI / n;
+ m += PI / ifr->order;
}
/* high pass or band reject
*/
- if( type >= 3 )
+ if (ifr->type >= 3)
{
- /* map s => 1/s
- */
- for( j=0; j<np; j++ )
+ /* map s => 1/s */
+ for (j = 0; j < ds->n_poles; j++)
{
ir = j + j;
ii = ir + 1;
- b = zs[ir]*zs[ir] + zs[ii]*zs[ii];
+ double b = zs[ir]*zs[ir] + zs[ii]*zs[ii];
zs[ir] = zs[ir] / b;
zs[ii] = zs[ii] / b;
}
/* The zeros at infinity map to the origin.
*/
- nz = np;
- if( type == 4 )
+ ds->n_zeros = ds->n_poles;
+ if( ifr->type == 4 )
{
- nz += n/2;
+ ds->n_zeros += ifr->order / 2;
}
- for( j=0; j<nz; j++ )
+ for (j=0; j<ds->n_zeros; j++)
{
ir = ii + 1;
ii = ir + 1;
@@ -2060,39 +1291,39 @@
}
}
}
- if( kind == 2 )
+ if( ifr->kind == 2 )
{
/* For Chebyshev, find radii of two Butterworth circles
* See Gold & Rader, page 60
*/
- rho = (phi - 1.0)*(phi+1); /* rho = eps^2 = {sqrt(1+eps^2)}^2 - 1 */
- eps = sqrt(rho);
- /* sqrt( 1 + 1/eps^2 ) + 1/eps = {sqrt(1 + eps^2) + 1} / eps
+ rho = (phi - 1.0)*(phi+1); /* rho = ds->ripple_epsilon^2 = {sqrt(1+ds->ripple_epsilon^2)}^2 - 1 */
+ ds->ripple_epsilon = sqrt(rho);
+ /* sqrt( 1 + 1/ds->ripple_epsilon^2 ) + 1/ds->ripple_epsilon = {sqrt(1 + ds->ripple_epsilon^2) + 1} / ds->ripple_epsilon
*/
- phi = (phi + 1.0) / eps;
- EVERBOSE ("Chebychev: phi-before=%.20g ripple=%.20g\n", phi, eps); // BSE info
- phi = pow( phi, 1.0/rn ); /* raise to the 1/n power */
- EVERBOSE ("Chebychev: phi-raised=%.20g rn=%.20g\n", phi, rn); // BSE info
- b = 0.5 * (phi + 1.0/phi); /* y coordinates are on this circle */
- a = 0.5 * (phi - 1.0/phi); /* x coordinates are on this circle */
- if( n & 1 )
+ phi = (phi + 1.0) / ds->ripple_epsilon;
+ EVERBOSE ("Chebychev: phi-before=%.20g ripple=%.20g\n", phi, ds->ripple_epsilon); // BSE info
+ phi = pow (phi, 1.0 / ifr->order); /* raise to the 1/n power */
+ EVERBOSE ("Chebychev: phi-raised=%.20g rn=%.20g\n", phi, ifr->order*1.0); // BSE info
+ double b = 0.5 * (phi + 1.0/phi); /* y coordinates are on this circle */
+ double a = 0.5 * (phi - 1.0/phi); /* x coordinates are on this circle */
+ double m;
+ if (ifr->order & 1)
m = 0.0;
else
- m = PI / (2.0*n);
- for( i=0; i<np; i++ )
+ m = PI / (2.0 * ifr->order);
+ for (i = 0; i < ds->n_poles; i++)
{ /* poles */
- lr = i + i;
+ int lr = i + i;
zs[lr] = -a * cos(m);
zs[lr+1] = b * sin(m);
- m += PI / n;
+ m += PI / ifr->order;
}
/* high pass or band reject
*/
- if( type >= 3 )
+ if( ifr->type >= 3 )
{
- /* map s => 1/s
- */
- for( j=0; j<np; j++ )
+ /* map s => 1/s */
+ for (j = 0; j < ds->n_poles; j++)
{
ir = j + j;
ii = ir + 1;
@@ -2102,12 +1333,12 @@
}
/* The zeros at infinity map to the origin.
*/
- nz = np;
- if( type == 4 )
+ ds->n_zeros = ds->n_poles;
+ if( ifr->type == 4 )
{
- nz += n/2;
+ ds->n_zeros += ifr->order/2;
}
- for( j=0; j<nz; j++ )
+ for (j=0; j<ds->n_zeros; j++)
{
ir = ii + 1;
ii = ir + 1;
@@ -2116,51 +1347,52 @@
}
}
}
- if( kind == 3 ) /* elliptic filter -- stw */
+ if( ifr->kind == 3 ) /* elliptic filter -- stw */
{
- nz = n/2;
- ellpj( u, 1.0-m, &sn1, &cn1, &dn1, &phi1 );
- for( i=0; i<ARRSIZ; i++ )
+ double m = ds->elliptic_m;
+ ds->n_zeros = ifr->order/2;
+ ellpj (ds->elliptic_u, 1.0 - m, &sn1, &cn1, &dn1, &phi1);
+ for (i=0; i<ARRSIZ; i++)
zs[i] = 0.0;
- for( i=0; i<nz; i++ )
+ for (i = 0; i < ds->n_zeros; i++)
{ /* zeros */
- a = n - 1 - i - i;
- b = (Kk * a) / rn;
+ double a = ifr->order - 1 - i - i;
+ double b = (Kk * a) / ifr->order;
ellpj( b, m, &sn, &cn, &dn, &phi );
- lr = 2*np + 2*i;
+ int lr = 2 * ds->n_poles + 2 * i;
zs[ lr ] = 0.0;
- a = wc/(k*sn); /* k = sqrt(m) */
+ a = ds->wc / (ds->elliptic_k * sn); /* elliptic_k = sqrt(m) */
zs[ lr + 1 ] = a;
}
- for( i=0; i<np; i++ )
+ for (i = 0; i < ds->n_poles; i++)
{ /* poles */
- a = n - 1 - i - i;
- b = a * Kk / rn;
+ double a = ifr->order - 1 - i - i;
+ double b = a * Kk / ifr->order;
ellpj( b, m, &sn, &cn, &dn, &phi );
- r = k * sn * sn1;
+ double r = ds->elliptic_k * sn * sn1;
b = cn1*cn1 + r*r;
- a = -wc*cn*dn*sn1*cn1/b;
- lr = i + i;
+ a = -ds->wc*cn*dn*sn1*cn1/b;
+ int lr = i + i;
zs[lr] = a;
- b = wc*sn*dn1/b;
+ b = ds->wc*sn*dn1/b;
zs[lr+1] = b;
}
- if( type >= 3 )
+ if( ifr->type >= 3 )
{
- nt = np + nz;
- for( j=0; j<nt; j++ )
+ int nt = ds->n_poles + ds->n_zeros;
+ for (j = 0; j < nt; j++)
{
ir = j + j;
ii = ir + 1;
- b = zs[ir]*zs[ir] + zs[ii]*zs[ii];
+ double b = zs[ir]*zs[ir] + zs[ii]*zs[ii];
zs[ir] = zs[ir] / b;
zs[ii] = zs[ii] / b;
}
- while( np > nz )
+ while (ds->n_poles > ds->n_zeros)
{
ir = ii + 1;
ii = ir + 1;
- nz += 1;
+ ds->n_zeros += 1;
zs[ir] = 0.0;
zs[ii] = 0.0;
}
@@ -2168,26 +1400,20 @@
}
printf( "s plane poles:\n" );
j = 0;
- for( i=0; i<np+nz; i++ )
+ for (i = 0; i < ds->n_poles + ds->n_zeros; i++)
{
- a = zs[j];
+ double a = zs[j];
++j;
- b = zs[j];
+ double b = zs[j];
++j;
printf( "%.9E %.9E\n", a, b );
- if( i == np-1 )
+ if (i == ds->n_poles - 1)
printf( "s plane zeros:\n" );
}
return 0;
}
-
-
-
-
-
-/* cay()
- *
+/* jacobi_theta_by_nome():
* Find parameter corresponding to given nome by expansion
* in theta functions:
* AMS55 #16.38.5, 16.38.7
@@ -2208,67 +1434,54 @@
* 1/2
* Given q, this program returns m .
*/
-double cay(q)
- double q;
+static double
+jacobi_theta_by_nome (double q)
{
- double a, b, p, r;
- double t1, t2;
-
- a = 1.0;
- b = 1.0;
- r = 1.0;
- p = q;
-
+ double t1, a = 1.0, b = 1.0, r = 1.0, p = q;
do
{
r *= p;
a += 2.0 * r;
- t1 = fabs( r/a );
+ t1 = fabs (r / a);
r *= p;
b += r;
p *= q;
- t2 = fabs( r/b );
- if( t2 > t1 )
+ double t2 = fabs (r / b);
+ if (t2 > t1)
t1 = t2;
}
- while( t1 > MACHEP );
-
- a = b/a;
- a = 4.0 * sqrt(q) * a * a; /* see above formulas, solved for m */
- return(a);
+ while (t1 > MACHEP);
+ a = b / a;
+ a = 4.0 * sqrt (q) * a * a; /* see above formulas, solved for m */
+ return a;
}
-
-
-
-/* zpln.c
- * Program to convert s plane poles and zeros to the z plane.
- */
-
-extern cmplx cone;
-
-int zplna()
+/* convert s plane poles and zeros to the z plane. */
+static int
+convert_s_plane_to_z_plane (const BseIIRFilterRequirements *ifr,
+ DesignState *ds)
{
- cmplx r, cnum, cden, cwc, ca, cb, b4ac;
+ Complex r, cnum, cden, cwc, ca, cb, b4ac;
double C;
- if( kind == 3 )
- C = c;
+ if (ifr->kind == 3)
+ C = ds->tan_angle_frequency;
else
- C = wc;
+ C = ds->wc;
- for( i=0; i<ARRSIZ; i++ )
+ int i;
+ for (i = 0; i < ARRSIZ; i++)
{
z[i].r = 0.0;
z[i].i = 0.0;
}
- nc = np;
+ int nc = ds->n_poles;
jt = -1;
ii = -1;
- for( icnt=0; icnt<2; icnt++ )
+ for (icnt=0; icnt<2; icnt++)
{
/* The maps from s plane to z plane */
do
@@ -2278,7 +1491,7 @@
r.r = zs[ir];
r.i = zs[ii];
- switch( type )
+ switch( ifr->type )
{
case 1:
case 3:
@@ -2295,7 +1508,7 @@
cden.r = 1 - C * r.r;
cden.i = -C * r.i;
jt += 1;
- cdiv( &cden, &cnum, &z[jt] );
+ Cdiv( &cden, &cnum, &z[jt] );
if( r.i != 0.0 )
{
/* fill in complex conjugate root */
@@ -2319,30 +1532,30 @@
*
* and solve for the roots in the z plane.
*/
- if( kind == 2 )
- cwc.r = cbp;
+ if( ifr->kind == 2 )
+ cwc.r = ds->chebyshev_band_cbp;
else
- cwc.r = c;
+ cwc.r = ds->tan_angle_frequency;
cwc.i = 0.0;
- cmul( &r, &cwc, &cnum ); /* r wc */
- csub( &cnum, &cone, &ca ); /* a = 1 - r wc */
- cmul( &cnum, &cnum, &b4ac ); /* 1 - (r wc)^2 */
- csub( &b4ac, &cone, &b4ac );
+ Cmul( &r, &cwc, &cnum ); /* r wc */
+ Csub( &cnum, &COMPLEX_ONE, &ca ); /* a = 1 - r wc */
+ Cmul( &cnum, &cnum, &b4ac ); /* 1 - (r wc)^2 */
+ Csub( &b4ac, &COMPLEX_ONE, &b4ac );
b4ac.r *= 4.0; /* 4ac */
b4ac.i *= 4.0;
- cb.r = -2.0 * cgam; /* b */
+ cb.r = -2.0 * ds->cgam; /* b */
cb.i = 0.0;
- cmul( &cb, &cb, &cnum ); /* b^2 */
- csub( &b4ac, &cnum, &b4ac ); /* b^2 - 4 ac */
- csqrt( &b4ac, &b4ac );
+ Cmul( &cb, &cb, &cnum ); /* b^2 */
+ Csub( &b4ac, &cnum, &b4ac ); /* b^2 - 4 ac */
+ Csqrt( &b4ac, &b4ac );
cb.r = -cb.r; /* -b */
cb.i = -cb.i;
ca.r *= 2.0; /* 2a */
ca.i *= 2.0;
- cadd( &b4ac, &cb, &cnum ); /* -b + sqrt( b^2 - 4ac) */
- cdiv( &ca, &cnum, &cnum ); /* ... /2a */
+ Cadd( &b4ac, &cb, &cnum ); /* -b + sqrt( b^2 - 4ac) */
+ Cdiv( &ca, &cnum, &cnum ); /* ... /2a */
jt += 1;
- cmov( &cnum, &z[jt] );
+ Cmov( &cnum, &z[jt] );
if( cnum.i != 0.0 )
{
jt += 1;
@@ -2351,10 +1564,10 @@
}
if( (r.i != 0.0) || (cnum.i == 0) )
{
- csub( &b4ac, &cb, &cnum ); /* -b - sqrt( b^2 - 4ac) */
- cdiv( &ca, &cnum, &cnum ); /* ... /2a */
+ Csub( &b4ac, &cb, &cnum ); /* -b - sqrt( b^2 - 4ac) */
+ Cdiv( &ca, &cnum, &cnum ); /* ... /2a */
jt += 1;
- cmov( &cnum, &z[jt] );
+ Cmov( &cnum, &z[jt] );
if( cnum.i != 0.0 )
{
jt += 1;
@@ -2369,42 +1582,41 @@
if( icnt == 0 )
{
zord = jt+1;
- if( nz <= 0 )
+ if( ds->n_zeros <= 0 )
{
- if( kind != 3 )
+ if( ifr->kind != 3 )
return(0);
else
break;
}
}
- nc = nz;
+ nc = ds->n_zeros;
} /* end for() loop */
return 0;
}
-
-
-
-int zplnb()
+static int
+z_plane_zeros_poles_to_numerator_denomerator (const BseIIRFilterRequirements *ifr,
+ DesignState *ds)
{
- cmplx lin[2];
+ Complex lin[2];
lin[1].r = 1.0;
lin[1].i = 0.0;
- if( kind != 3 )
+ if( ifr->kind != 3 )
{ /* Butterworth or Chebyshev */
/* generate the remaining zeros */
while( 2*zord - 1 > jt )
{
- if( type != 3 )
+ if( ifr->type != 3 )
{
printf( "adding zero at Nyquist frequency\n" );
jt += 1;
z[jt].r = -1.0; /* zero at Nyquist frequency */
z[jt].i = 0.0;
}
- if( (type == 2) || (type == 3) )
+ if( (ifr->type == 2) || (ifr->type == 3) )
{
printf( "adding zero at 0 Hz\n" );
jt += 1;
@@ -2420,7 +1632,7 @@
jt += 1;
z[jt].r = -1.0; /* zero at Nyquist frequency */
z[jt].i = 0.0;
- if( (type == 2) || (type == 4) )
+ if( (ifr->type == 2) || (ifr->type == 4) )
{
jt += 1;
z[jt].r = 1.0; /* zero at 0 Hz */
@@ -2429,41 +1641,44 @@
}
}
printf( "order = %d\n", zord );
+
+ int j;
/* Expand the poles and zeros into numerator and
* denominator polynomials
*/
- for( icnt=0; icnt<2; icnt++ )
+ for (icnt=0; icnt<2; icnt++)
{
- for( j=0; j<ARRSIZ; j++ )
+ for (j=0; j<ARRSIZ; j++)
{
pp[j] = 0.0;
y[j] = 0.0;
}
pp[0] = 1.0;
- for( j=0; j<zord; j++ )
+ for (j=0; j<zord; j++)
{
- jj = j;
+ int jj = j;
if( icnt )
jj += zord;
- a = z[jj].r;
- b = z[jj].i;
- for( i=0; i<=j; i++ )
+ double a = z[jj].r;
+ double b = z[jj].i;
+ int i;
+ for (i = 0; i <= j; i++)
{
- jh = j - i;
+ int jh = j - i;
pp[jh+1] = pp[jh+1] - a * pp[jh] + b * y[jh];
y[jh+1] = y[jh+1] - b * pp[jh] - a * y[jh];
}
}
if( icnt == 0 )
{
- for( j=0; j<=zord; j++ )
+ for (j=0; j<=zord; j++)
aa[j] = pp[j];
}
}
/* Scale factors of the pole and zero polynomials */
- a = 1.0;
- switch( type )
+ double a = 1.0;
+ switch( ifr->type )
{
case 3:
a = -1.0;
@@ -2473,7 +1688,7 @@
pn = 1.0;
an = 1.0;
- for( j=1; j<=zord; j++ )
+ for (j=1; j<=zord; j++)
{
pn = a * pn + pp[j];
an = a * an + aa[j];
@@ -2481,23 +1696,23 @@
break;
case 2:
- gam = PI/2.0 - asin( cgam ); /* = acos( cgam ) */
- mh = zord/2;
+ gam = PI/2.0 - asin( ds->cgam ); /* = acos( cgam ) */
+ int mh = zord/2;
pn = pp[mh];
an = aa[mh];
- ai = 0.0;
+ double ai = 0.0;
if( mh > ((zord/4)*2) )
{
ai = 1.0;
pn = 0.0;
an = 0.0;
}
- for( j=1; j<=mh; j++ )
+ for (j=1; j<=mh; j++)
{
a = gam * j - ai * PI / 2.0;
- cng = cos(a);
- jh = mh + j;
- jl = mh - j;
+ double cng = cos (a);
+ int jh = mh + j;
+ int jl = mh - j;
pn = pn + cng * (pp[jh] + (1.0 - 2.0 * ai) * pp[jl]);
an = an + cng * (aa[jh] + (1.0 - 2.0 * ai) * aa[jl]);
}
@@ -2505,21 +1720,20 @@
return 0;
}
-
-
-
-int zplnc()
+static int
+gainscale_and_print_deno_nume_zeros2_poles2 (const BseIIRFilterRequirements *ifr, /* zplnc */
+ DesignState *ds)
{
-
- gain = an/(pn*scale);
- if( (kind != 3) && (pn == 0) )
+ int j;
+ gain = an/(pn*ds->gain_scale);
+ if( (ifr->kind != 3) && (pn == 0) )
gain = 1.0;
printf( "constant gain factor %23.13E\n", gain );
- for( j=0; j<=zord; j++ )
+ for (j = 0; j <= zord; j++)
pp[j] = gain * pp[j];
printf( "z plane Denominator Numerator\n" );
- for( j=0; j<=zord; j++ )
+ for (j=0; j<=zord; j++)
{
printf( "%2d %17.9E %17.9E\n", j, aa[j], pp[j] );
}
@@ -2528,35 +1742,33 @@
* so that it can be implemented without stability problems -- stw
*/
printf( "poles and zeros with corresponding quadratic factors\n" );
- for( j=0; j<zord; j++ )
+ for (j=0; j<zord; j++)
{
- a = z[j].r;
- b = z[j].i;
+ double a = z[j].r;
+ double b = z[j].i;
if( b >= 0.0 )
{
printf( "pole %23.13E %23.13E\n", a, b );
- quadf( a, b, 1 );
+ print_quadratic_factors (ifr, ds, a, b, 1);
}
- jj = j + zord;
+ int jj = j + zord;
a = z[jj].r;
b = z[jj].i;
if( b >= 0.0 )
{
printf( "zero %23.13E %23.13E\n", a, b );
- quadf( a, b, 0 );
+ print_quadratic_factors (ifr, ds, a, b, 0);
}
}
return 0;
}
-
-
-
-/* display quadratic factors
- */
-int quadf( x, y, pzflg )
- double x, y;
- int pzflg; /* 1 if poles, 0 if zeros */
+/* display quadratic factors */
+static int
+print_quadratic_factors (const BseIIRFilterRequirements *ifr,
+ DesignState *ds,
+ double x, double y,
+ int pzflg) /* 1 if poles, 0 if zeros */
{
double a, b, r, f, g, g0;
@@ -2576,7 +1788,7 @@
/* resonant frequency */
r = sqrt(b);
f = PI/2.0 - asin( -a/(2.0*r) );
- f = f * fs / (2.0 * PI );
+ f = f * ifr->sampling_frequency / (2.0 * PI );
/* gain at resonance */
g = 1.0 + r;
g = g*g - (a*a/r);
@@ -2586,9 +1798,9 @@
else
{
/* It is really a first-order network.
- * Give the gain at fnyq and D.C.
+ * Give the gain at ds->nyquist_frequency and D.C.
*/
- f = fnyq;
+ f = ds->nyquist_frequency;
g = 1.0 - a;
g0 = 1.0 + a;
}
@@ -2608,40 +1820,37 @@
return 0;
}
-
-
-/* Print table of filter frequency response
- */
-void
-print_filter_table (void)
+/* Print table of filter frequency response */
+static void
+print_filter_table (const BseIIRFilterRequirements *ifr,
+ DesignState *ds)
{
- double f, limit = 0.05 * fnyq * 21;
+ double f, limit = 0.05 * ds->nyquist_frequency * 21;
for (f=0; f < limit; f += limit / 21.)
{
- double r = response( f, gain );
+ double r = response (ifr, ds, f, gain);
if( r <= 0.0 )
r = -999.99;
else
- r = 2.0 * dbfac * log( r );
+ r = 2.0 * DECIBELL_FACTOR * log( r );
printf( "%10.1f %10.2f\n", f, r );
- // f = f + 0.05 * fnyq;
+ // f = f + 0.05 * ds->nyquist_frequency;
}
}
-
-/* Calculate frequency response at f Hz
- * mulitplied by amp
- */
-double response( f, amp )
- double f, amp;
+/* Calculate frequency response at f Hz mulitplied by amp */
+static double
+response (const BseIIRFilterRequirements *ifr,
+ DesignState *ds,
+ double f, double amp)
{
- cmplx x, num, den, w;
+ Complex x, num, den, w;
double u;
int j;
/* exp( j omega T ) */
- u = 2.0 * PI * f /fs;
+ u = 2.0 * PI * f /ifr->sampling_frequency;
x.r = cos(u);
x.i = sin(u);
@@ -2649,43 +1858,66 @@
num.i = 0.0;
den.r = 1.0;
den.i = 0.0;
- for( j=0; j<zord; j++ )
+ for (j=0; j<zord; j++)
{
- csub( &z[j], &x, &w );
- cmul( &w, &den, &den );
- csub( &z[j+zord], &x, &w );
- cmul( &w, &num, &num );
+ Csub( &z[j], &x, &w );
+ Cmul( &w, &den, &den );
+ Csub( &z[j+zord], &x, &w );
+ Cmul( &w, &num, &num );
}
- cdiv( &den, &num, &w );
+ Cdiv( &den, &num, &w );
w.r *= amp;
w.i *= amp;
- u = cabs( &w );
+ u = Cabs( &w );
return(u);
}
-
-
-
-/* Get a number from keyboard.
- * Display previous value and keep it if user just hits <CR>.
- */
-int getnum( line, val )
- char *line;
- double *val;
+static double
+my_getnum (const char *text)
{
- char s[40];
-
- printf( "%s = %.9E ? ", line, *val );
+ printf ("%s ? ", text);
+ char s[4096];
if (!fgets (s, sizeof (s), stdin))
exit (0);
- if( s[0] != '\0' )
+ double val = 0;
+ sscanf (s, "%lf", &val);
+ return val;
+}
+
+
+int
+main (int argc,
+ char *argv[])
+{
+ init_constants();
+ BseIIRFilterRequirements ifr = { 0 };
+ DesignState ds = default_design_state;
+ ifr.kind = my_getnum ("kind");
+ ifr.type = my_getnum ("type");
+ ifr.order = my_getnum ("order");
+ if (ifr.kind > BSE_IIR_FILTER_BUTTERWORTH) /* not Butterworth */
+ ifr.passband_ripple_db = my_getnum ("passband_ripple_db");
+ ifr.sampling_frequency = my_getnum ("sampling_frequency");
+ ifr.passband_edge = my_getnum ("passband_edge");
+ if (ifr.type == BSE_IIR_FILTER_BAND_PASS ||
+ ifr.type == BSE_IIR_FILTER_BAND_STOP)
+ ifr.passband_edge2 = my_getnum ("passband_edge2");
+ if (ifr.kind == BSE_IIR_FILTER_ELLIPTIC)
+ ifr.stopband_db = ifr.stopband_edge = my_getnum ("stopband_edge or stopband_db");
+ printf ("\n");
+ const char *errmsg = iir_filter_design (&ifr, &ds);
+ fflush (stdout);
+ fflush (stderr);
+ // VERBOSE ("DEBUG: %.20g %.20g %.20g %.20g %.20g\n", a, cos(a), cang, ds->cgam, 0.);
+ if (errmsg)
{
- sscanf( s, "%lf", val );
- printf( "%.9E\n", *val );
+ fprintf (stderr, "Invalid specification: %s\n", errmsg);
+ fflush (stderr);
+ return 1;
}
+
return 0;
}
-/* === ellf.c - end === */
-/* compile with: gcc -Wall -O2 -g bseiirfilter.c -lm -o ellf */
+/* compile with: gcc -Wall -O2 -g bseellipticfilter.c -lm -o bseellipticfilter */
Modified: trunk/bse/bseellipticfilter.h
===================================================================
--- trunk/bse/bseellipticfilter.h 2006-10-15 22:48:08 UTC (rev 3971)
+++ trunk/bse/bseellipticfilter.h 2006-10-16 01:05:27 UTC (rev 3972)
@@ -19,10 +19,69 @@
#ifndef BSE_IIR_FILTER_H__
#define BSE_IIR_FILTER_H__
-#include <bse/bsemath.h>
+// FIXME: #include <bse/bsemath.h>
+#include <stdbool.h> // FIXME
+typedef unsigned int uint; // FIXME
-BIRNET_EXTERN_C_BEGIN();
+// FIXME: BIRNET_EXTERN_C_BEGIN();
-BIRNET_EXTERN_C_END();
+typedef enum {
+ BSE_IIR_FILTER_BUTTERWORTH = 1,
+ BSE_IIR_FILTER_CHEBYSHEV = 2,
+ BSE_IIR_FILTER_ELLIPTIC = 3,
+} BseIIRFilterKind;
+typedef enum {
+ BSE_IIR_FILTER_LOW_PASS = 1,
+ BSE_IIR_FILTER_BAND_PASS = 2,
+ BSE_IIR_FILTER_HIGH_PASS = 3,
+ BSE_IIR_FILTER_BAND_STOP = 4,
+} BseIIRFilterType;
+typedef struct {
+ BseIIRFilterKind kind;
+ BseIIRFilterType type;
+ uint order; /* >= 1 */
+ double passband_ripple_db; /* dB, not Butterworth */
+ double sampling_frequency; /* Hz, > 0.0 */
+ double passband_edge; /* Hz, > 0.0 && < nyquist_frequency */
+ double passband_edge2; /* Hz, > 0.0 && < nyquist_frequency, for BAND filters */
+ double stopband_db; /* dB < 0, elliptic only */
+ double stopband_edge; /* Hz, > 0.0 && < nyquist_frequency, elliptic only */
+} BseIIRFilterRequirements;
+typedef struct {
+ int n_poles;
+ int n_zeros;
+ double gain_scale;
+ double ripple_epsilon;
+ double nyquist_frequency;
+ double tan_angle_frequency;
+ double wc; /* tan_angle_frequency or normalized to 1.0 for elliptic */
+ double cgam; /* angle frequency temporary */
+ double stopband_edge; /* derived from ifr->stopband_edge or ifr->stopband_db */
+ double wr;
+ double elliptic_k;
+ double elliptic_u;
+ double elliptic_m;
+ double chebyshev_band_cbp;
+} DesignState;
+
+static const DesignState default_design_state = {
+ .n_poles = 0.0,
+ .n_zeros = 0.0,
+ .gain_scale = 0.0,
+ .ripple_epsilon = 0.0,
+ .nyquist_frequency = 0.0,
+ .tan_angle_frequency = 0.0,
+ .wc = 0.0,
+ .cgam = 0.0,
+ .stopband_edge = 2400,
+ .wr = 0.0,
+ .elliptic_k = 0.0,
+ .elliptic_u = 0.0,
+ .elliptic_m = 0.0,
+ .chebyshev_band_cbp = 0.0,
+};
+
+// FIXME: BIRNET_EXTERN_C_END();
+
#endif /* BSE_IIR_FILTER_H__ */ /* vim:set ts=8 sw=2 sts=2: */
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
