Bump version

In addition to C-style pointer-to-array and size arguments, allow for
C++-style std::vector<> interfaces. This also integrates nicely with
swig and allows the Python bindings to use native lists/tuples.

config

@@ -33,4 +33,5 @@ source "scripts/paths"
 
 # Project-specific variables below.
 [ -z "${CC}" ] && CC="gcc"
+[ -z "${CXX}" ] && CXX="g++"
 [ -z "${CFLAGS}" ] && CFLAGS="-g -O2 -W -Wall"

python/.gitignore

@@ -0,0 +1,5 @@
+.ipynb_checkpoints/
+__pycache__/
+build/
+iir.py
+pyiir_wrap.cxx

python/README

@@ -0,0 +1,25 @@
+Experimental Python bindings
+============================
+
+These bindings are currently experimental, as it is not yet clear how best
+to create and maintain them. They are built upon the C++ interface as the
+object oriented nature of it sits much better with Python.
+
+At present, it is necessary to compile and install the C++ library before
+working with the bindings.
+
+To generate the bindings, swig 3.0.2 has been used, but earlier versions
+may work too:
+
+  swig -c++ -python -py3 pyiir.i
+
+To build the extension, distutils is used:
+
+  python3 setup.py build_ext
+
+Note that the -L ${LIBDIR} option may be handy if libiir++ is not installed
+into the system library path.
+
+To install the extension:
+
+ python3 setup.py install

python/iir++.h

@@ -0,0 +1 @@
+../obj/iir++.h

python/pyiir.i

@@ -0,0 +1,142 @@
+%module iir
+%include <complex.i>
+
+/*
+ * List conversion
+ *
+ *  SWIG knows how to turn std::vector<T> func() into a function that returns a
+ *  Python list. We do need to tell it about all <T> types with %template,
+ *  however.
+ */
+%include <std_vector.i>
+%template(vectord) std::vector<double>;
+
+
+
+%{
+#include "iir++.h"
+%}
+
+
+
+/*
+ * Add bounds checking
+ *
+ *  The Python library is targetted at experimentation. It is therefore
+ *  reasonable to protect its object query and creation interfaces with code
+ *  that performs bounds checking on arguments.
+ */
+%exception IIR::Coeff::c {
+    if(arg2 < 0 || arg2 >= arg1->nc()) {
+        char buf[80];
+        snprintf(buf, sizeof(buf), "Index of c coefficient (%d) out of bounds "
+            "(0,%d).", arg2, arg1->nc() - 1);
+        PyErr_SetString(PyExc_RuntimeError, buf);
+        return NULL;
+    }
+    $action
+}
+
+%exception IIR::Coeff::d {
+    if(arg2 < 0 || arg2 >= arg1->nd()) {
+        char buf[80];
+        snprintf(buf, sizeof(buf), "Index of d coefficient (%d) out of bounds "
+            "(0,%d).", arg2, arg1->nd() - 1);
+        PyErr_SetString(PyExc_RuntimeError, buf);
+        return NULL;
+    }
+    $action
+}
+
+%exception IIR::Filter::getCoeffSet {
+    if(arg2 < 0 || arg2 >= arg1->numCoeffSet()) {
+        char buf[80];
+        snprintf(buf, sizeof(buf), "Index of coefficient set (%d) out of bounds "
+            "(0,%d).", arg2, arg1->numCoeffSet() - 1);
+        PyErr_SetString(PyExc_RuntimeError, buf);
+        return NULL;
+    }
+    $action
+}
+
+%exception IIR::ButterworthLowPass {
+    if(arg1 < 1) {
+        char buf[80];
+        snprintf(buf, sizeof(buf), "Invalid filter order (%d). Must be > 0.",
+            arg1);
+        PyErr_SetString(PyExc_RuntimeError, buf);
+        return NULL;
+    }
+    $action
+}
+
+%exception IIR::ButterworthHighPass {
+    if(arg1 < 1) {
+        char buf[80];
+        snprintf(buf, sizeof(buf), "Invalid filter order (%d). Must be > 0.",
+            arg1);
+        PyErr_SetString(PyExc_RuntimeError, buf);
+        return NULL;
+    }
+    $action
+}
+
+%exception IIR::ButterworthBandPass {
+    if(arg1 < 1) {
+        char buf[80];
+        snprintf(buf, sizeof(buf), "Invalid filter order (%d). Must be > 0.",
+            arg1);
+        PyErr_SetString(PyExc_RuntimeError, buf);
+        return NULL;
+    }
+    if(arg3 > arg4) {
+        PyErr_SetString(PyExc_RuntimeError, "Low corner frequency exceeds "
+            "high corner frequency.");
+        return NULL;
+    }
+    $action
+}
+
+%exception IIR::ButterworthBandStop {
+    if(arg1 < 1) {
+        char buf[80];
+        snprintf(buf, sizeof(buf), "Invalid filter order (%d). Must be > 0.",
+            arg1);
+        PyErr_SetString(PyExc_RuntimeError, buf);
+        return NULL;
+    }
+    if(arg3 > arg4) {
+        PyErr_SetString(PyExc_RuntimeError, "Low corner frequency exceeds "
+            "high corner frequency.");
+        return NULL;
+    }
+    $action
+}
+
+
+
+/*
+ * Invalid filter string handling
+ *
+ *  The C++ library uses a "valid()" member in IIR::Filter to avoid throwing an
+ *  exception if the constructor is passed an invalid filter string. Python's
+ *  exceptions are a much nicer fit, however, so use those instead.
+ *
+ *  Notes:
+ *   1. Doesn't match if we prepend it with the namespace.
+ *   2. Applies to all constructors.
+ */
+%exception Filter {
+    $action
+    if(!result->valid()) {
+        PyErr_SetString(PyExc_RuntimeError, "Invalid filter description string.");
+        delete result;
+        return NULL;
+    }
+}
+
+
+
+%include "iir++.h"
+
+// vim: ts=4:sw=4

python/setup.py

@@ -0,0 +1,16 @@
+#!/usr/bin/env python3
+
+from distutils.core import setup, Extension
+
+iir_module = Extension('_iir',
+	sources=['pyiir_wrap.cxx'],
+	libraries=['iir++'],
+	)
+
+setup(name = 'iir',
+	version = '1.0.3',
+	author = 'Laurence Withers',
+	description = 'IIR filter',
+	ext_modules = [iir_module],
+	py_modules = ['iir'],
+	)

run-test.sh

@@ -1,7 +1,7 @@
 #!/bin/bash
 # libiir/test.sh
 #
-#  Copyright: ©2010, Laurence Withers.
+#  Copyright: ©2010–2011, Laurence Withers.
 #  Author: Laurence Withers <l@lwithers.me.uk>
 #  License: GPLv3
 #
@@ -48,7 +48,7 @@ then
     exit 0
 fi
 
-run_test $*
+run_test "$@"
 
 # kate: replace-trailing-space-save true; space-indent true; tab-width 4;
 # vim: expandtab:ts=4:sw=4

scripts/.gitignore

@@ -7,6 +7,7 @@ build.docs.none
 build.files.none
 build.firmware.gpasm
 build.firmware.sdcc
+build.header.c
 build.lib.c
 build.lib.c++
 build.make.none

src/docs/MainPage.dox

@@ -1,16 +1,16 @@
 /* libiir/src/docs/MainPage.dox
  *
- *  Copyright: ©2010, Laurence Withers.
+ *  Copyright: ©2010–2014, Laurence Withers.
  *  Author: Laurence Withers <l@lwithers.me.uk>
  *  License: GPLv3
 */
 
 /*! \mainpage
 
-This library allows the implementation of arbitrary IIR filters in C. It has
+This library allows the implementation of arbitrary %IIR filters in C. It has
 functions for generating and manipulating filters in terms of coefficients, for
 chaining arbitrary filters together, and for generating coefficients for some
-common types of filter. See \ref iir_structure for a definition of the IIR
+common types of filter. See \ref iir_structure for a definition of the %IIR
 filter equation.
 
 \section creation Filter creation
@@ -23,7 +23,7 @@ Otherwise, the library user must first create a set of coefficients using
 \ref iir_coeff_new(). Any number of filters can then be instantiated using that
 set of coefficients with \ref iir_filter_new(), or the coefficients can be
 chained on to the end of an existing filter instance with
-\ref iir_filter_chain(). See \ref common_filters for functions to generate
+\ref iir_filter_chain(). See \ref libiir_common_filters for functions to generate
 coefficients.
 
 \section operation Filter operation
@@ -47,6 +47,11 @@ installed) which will generate a Bode plot for a given filter chain. Note the
 phase response can be a little rough due to the simplistic time-domain analysis
 of the output signal's phase.
 
+\section cpp C++
+
+Two C++ wrapper objects, \ref IIR::Coeff and \ref IIR::Filter, are supplied which provide
+a true object-oriented interface on top of the C implementation.
+
 */
 
 /* options for text editors

src/docs/iir_structure.dox

@@ -1,6 +1,6 @@
 /* libiir/src/docs/iir_structure.dox
  * 
- *  Copyright: ©2010, Laurence Withers.
+ *  Copyright: ©2010–2014, Laurence Withers.
  *  Author: Laurence Withers <l@lwithers.me.uk>
  *  License: GPLv3
 */
@@ -14,16 +14,17 @@ For the purposes of this library, the following notation is used:
 \li <code>c[n]</code>: array of \c x(t) coefficients
 \li <code>d[n]</code>: array of \c y(t) coefficients
 
-This leads to a general IIR filter equation:
+This leads to a general %IIR filter equation:
 
 <code>y(t) = x(t).c[0] + x(t-1).c[1] + … + x(t-N).c[N] - y(t-1).d[0] - y(t-2).d[1] - … - y(t-1-M).d[M]</code>
 
-For initial conditions, the library sets <code>y(t)</code> for t &lt; 0 to 0,
-and <code>x(t)</code> for t &lt; 0 to <code>x(0)</code>.
+For initial conditions, the library sets <code>y(t)</code> for t &lt; 0 to
+<code>x(0).G</code> (where <code>G</code> is the DC gain of the filter, or 0 if
+it looks like the filter is unstable at DC), and <code>x(t)</code> for t &lt; 0
+to <code>x(0)</code>.
 
 */
 
 /* options for text editors
-kate: replace-trailing-space-save true; space-indent true; tab-width 4;
 vim: expandtab:ts=4:sw=4:syntax=doxygen
 */

src/docs/string_desc.dox

@@ -1,17 +1,17 @@
 /* libiir/src/docs/string_desc.dox
  *
- *  Copyright: ©2010, Laurence Withers.
+ *  Copyright: ©2010–2014, Laurence Withers.
  *  Author: Laurence Withers <l@lwithers.me.uk>
  *  License: GPLv3
 */
 
 /*! \page string_desc Describing IIR filters as strings
 
-This library allows the user to describe an IIR filter chain as a string, which
+This library allows the user to describe an %IIR filter chain as a string, which
 is useful to allow configurable filtering using e.g. a configuration file. This
 page describes the format of such strings.
 
-The string is first split into individual IIR filters. Each filter is written as
+The string is first split into individual %IIR filters. Each filter is written as
 \c type(params) and separated by whitespace. The string must contain at least
 one filter but may contain an arbitrary number.
 
@@ -23,7 +23,7 @@ one filter but may contain an arbitrary number.
 
 \subsection string_desc_coeff Raw coefficients
 
-An IIR filter may be specified as raw coefficients, in which case the \c type
+An %IIR filter may be specified as raw coefficients, in which case the \c type
 is \c raw and the \c params consists of a string:
 
 <code>c[0],c[1],…,c[n]/d[0],d[1]…d[n]</code>

src/libiir++/.params

@@ -0,0 +1 @@
+lib c++ libiir++ iir++.h

src/libiir++/000_TopHeader.h

@@ -0,0 +1,32 @@
+/* libiir/src/libiir++/000_TopHeader.h
+ * 
+ *  Copyright: ©2014, Laurence Withers.
+ *  Author: Laurence Withers <l@lwithers.me.uk>
+ *  License: GPLv3
+ */
+
+#ifndef HEADER_libiirpp
+#define HEADER_libiirpp
+
+// standard includes, or includes needed for type declarations
+#include <complex>
+#include <vector>
+
+// forward declares of opaque C structures (to avoid pulling in header)
+extern "C" {
+struct iir_coeff_t;
+struct iir_filter_t;
+}
+
+/*! \defgroup libiirpp C++ library */
+
+/*! \brief All %IIR library functions and classes.
+
+\ingroup libiirpp
+
+*/
+namespace IIR {
+
+/* options for text editors
+vim: expandtab:ts=4:sw=4:syntax=cpp.doxygen
+*/

src/libiir++/000_TopSource.cc

@@ -0,0 +1,26 @@
+/* libiir/src/libiir++/000_TopSource.cc
+ *
+ *  Copyright: ©2014, Laurence Withers.
+ *  Author: Laurence Withers <l@lwithers.me.uk>
+ *  License: GPLv3
+ */
+
+// we need to include "iir.h" first, as it pulls in <complex.h>, which we need
+// to take effect before "iir++.h" pulls in <complex>
+#include "iir.h"
+
+// now remove the preprocessor definition of complex to _Complex, which is fine
+// for the C header but not good for the C++ header
+#undef complex
+
+#include "iir++.h"
+
+// Below are all the includes used throughout the library.
+#include <alloca.h>
+#include <stdlib.h>
+
+namespace IIR {
+
+/* options for text editors
+vim: expandtab:ts=4:sw=4
+*/

src/libiir++/100_Coeff.cc

@@ -0,0 +1,98 @@
+/* libiir/src/libiir++/100_Coeff.cc
+ * 
+ *  Copyright: ©2014, Laurence Withers.
+ *  Author: Laurence Withers <l@lwithers.me.uk>
+ *  License: GPLv3
+ */
+
+
+
+Coeff::Coeff(int nc, const double* c, int nd, const double* d)
+{
+    cc_ = iir_coeff_new(nc, c, nd, d);
+}
+
+
+
+Coeff::Coeff(const std::vector<double>& c, const std::vector<double>& d)
+{
+    int nc, nd;
+    double* ac, * ad;
+
+    nc = c.size();
+    ac = static_cast<double*>(alloca(nc * sizeof(double)));
+    for(int i = 0; i < nc; ++i) ac[i] = c[i];
+
+    nd = d.size();
+    ad = static_cast<double*>(alloca(nd * sizeof(double)));
+    for(int i = 0; i < nd; ++i) ad[i] = d[i];
+
+    cc_ = iir_coeff_new(nc, ac, nd, ad);
+}
+
+
+
+Coeff::Coeff(const struct iir_coeff_t* cc)
+{
+    cc_ = iir_coeff_copy(cc);
+}
+
+
+
+Coeff::Coeff(const Coeff& other)
+{
+    cc_ = iir_coeff_copy(other.cc_);
+}
+
+
+
+Coeff::~Coeff()
+{
+    iir_coeff_free(cc_);
+}
+
+
+
+int
+Coeff::nc() const
+{
+    return iir_coeff_get_nc(cc_);
+}
+
+
+
+double
+Coeff::c(int idx) const
+{
+    return iir_coeff_get_c(cc_, idx);
+}
+
+
+
+int
+Coeff::nd() const
+{
+    return iir_coeff_get_nd(cc_);
+}
+
+
+
+double
+Coeff::d(int idx) const
+{
+    return iir_coeff_get_d(cc_, idx);
+}
+
+
+
+std::complex<double>
+Coeff::response(double freq) const
+{
+    return std::complex<double>(iir_response_c(cc_, freq));
+}
+
+
+
+/* options for text editors
+vim: expandtab:ts=4:sw=4:syntax=cpp.doxygen
+*/

src/libiir++/100_Coeff.h

@@ -0,0 +1,130 @@
+/* libiir/src/libiir++/100_Coeff.h
+ * 
+ *  Copyright: ©2014, Laurence Withers.
+ *  Author: Laurence Withers <l@lwithers.me.uk>
+ *  License: GPLv3
+ */
+
+
+
+/*! \defgroup libiirpp_coeff Coefficient handling
+
+\ingroup libiirpp
+
+A set of coefficients for a single stage %IIR filter is represented by an
+instance of the \ref Coeff class. This class is often instantiated through
+one of the common filter generator functions (see \ref libiirpp_common_filters).
+
+*/
+
+
+
+/*! \brief %IIR filter coefficient set
+
+\ingroup libiirpp_coeff
+
+An instance of this class represents the set of coefficients used to implement
+a single-phase %IIR filter. See \ref iir_structure for more information.
+
+*/
+class Coeff {
+public:
+    /*! \brief Constructor with raw coefficients.
+
+    \param nc Number of \a c coefficients, ≥ 1.
+    \param c Array of \a c coefficients.
+    \param nd Number of \a d coefficients, ≥ 1.
+    \param d Array of \a d coefficients.
+
+    */
+    Coeff(int nc, const double* c, int nd, const double* d);
+
+    /*! \brief Constructor with raw coefficients.
+
+    \param c Array of \a c coefficients, size ≥ 1.
+    \param d Array of \a d coefficients, size ≥ 1.
+
+    */
+    Coeff(const std::vector<double>& c, const std::vector<double>& d);
+
+    /*! \brief Constructor from C library object.
+
+    \param cc Pointer to existing coefficient object.
+
+    */
+    Coeff(const struct iir_coeff_t* cc);
+
+    /*! \brief Copy constructor.
+
+    \param other Object to copy.
+
+    Performs a full deep copy of the set of coefficients from \a other.
+
+    */
+    Coeff(const Coeff& other);
+
+    /*! \brief Destructor. */
+    virtual ~Coeff();
+
+
+
+    /*! \brief Count number of \a c coefficients.
+
+    \returns Number of \a c coefficients, ≥ 1.
+
+    */
+    int nc() const;
+
+    /*! \brief Retrieve \a c coefficient by index.
+
+    \param idx Index of coefficient (0 ≤ \a idx < \ref nc()).
+    \returns Value of \a c coefficient at index \a idx.
+
+    */
+    double c(int idx) const;
+
+    /*! \brief Count number of \a d coefficients.
+
+    \returns Number of \a d coefficients, ≥ 1.
+
+    */
+    int nd() const;
+
+    /*! \brief Retrieve \a d coefficient by index.
+
+    \param idx Index of coefficient (0 ≤ \a idx < \ref nd()).
+    \returns Value of \a d coefficient at index \a idx.
+
+    */
+    double d(int idx) const;
+
+
+
+    /*! \brief Evaluate response of filter at given frequency.
+
+    \param freq Frequency (0 ≤ \a freq ≤ 1).
+    \returns Complex response of filter at given frequency.
+
+    Returns the steady-state response of the filter to a pure input signal of
+    the given frequency. \a freq is expressed as a fraction of the intended
+    sampling rate (with 0.5 being the Nyquist frequency).
+
+    */
+    std::complex<double> response(double freq) const;
+
+private:
+    // allow the Filter class to access cc_ directly
+    friend class Filter;
+
+    // our internal representation is the C object
+    struct iir_coeff_t* cc_;
+
+    // disallow assignment; object is immutable
+    Coeff& operator=(const Coeff& other);
+};
+
+
+
+/* options for text editors
+vim: expandtab:ts=4:sw=4:syntax=cpp.doxygen
+*/

src/libiir++/200_Filter.cc

@@ -0,0 +1,119 @@
+/* libiir/src/libiir++/200_Filter.cc
+ * 
+ *  Copyright: ©2014, Laurence Withers.
+ *  Author: Laurence Withers <l@lwithers.me.uk>
+ *  License: GPLv3
+ */
+
+
+
+Filter::Filter(const Coeff& coeff)
+{
+    coeff_.push_back(new Coeff(coeff));
+    fi_ = iir_filter_new(coeff.cc_);
+}
+
+
+
+Filter::Filter(const Filter& other)
+{
+    for(int pos = 0, max = other.coeff_.size(); pos < max; ++pos) {
+        Coeff* cc = new Coeff(*other.coeff_[pos]);
+        coeff_.push_back(cc);
+
+        if(!pos) fi_ = iir_filter_new(cc->cc_);
+        else iir_filter_chain(fi_, cc->cc_);
+    }
+}
+
+
+
+Filter::Filter(const char* desc)
+{
+    fi_ = iir_parse(desc);
+    if(!fi_) return;
+
+    for(int pos = 0, max = iir_filter_coeff_sets(fi_); pos < max; ++pos) {
+        struct iir_coeff_t* c = iir_filter_coeff_set(fi_, pos);
+        coeff_.push_back(new Coeff(c));
+        iir_coeff_free(c);
+    }
+}
+
+
+
+Filter::~Filter()
+{
+    for(int pos = 0, max = coeff_.size(); pos < max; ++pos) delete coeff_[pos];
+    iir_filter_free(fi_);
+}
+
+
+
+bool
+Filter::valid() const
+{
+    return fi_;
+}
+
+
+
+double
+Filter::f(double x)
+{
+    return iir_filter(fi_, x);
+}
+
+
+
+int
+Filter::numCoeffSet() const
+{
+    return coeff_.size();
+}
+
+
+
+const Coeff&
+Filter::getCoeffSet(int idx) const
+{
+    return *(coeff_[idx]);
+}
+
+
+
+void
+Filter::chainCoeffSet(const Coeff& coeff)
+{
+    coeff_.push_back(new Coeff(coeff));
+    iir_filter_chain(fi_, coeff.cc_);
+}
+
+
+
+std::complex<double>
+Filter::response(double freq) const
+{
+    return iir_response(fi_, freq);
+}
+
+
+
+std::vector<double>
+Filter::responseMinus3dB(double gain) const
+{
+    double* cf;
+    int num_cf, i;
+
+    iir_response_minus_3dB(fi_, gain, &cf, &num_cf);
+    std::vector<double> ret(num_cf);
+    for(i = 0; i < num_cf; ++i) ret[i] = cf[i];
+    free(cf);
+    return ret;
+}
+
+
+
+/* options for text editors
+vim: expandtab:ts=4:sw=4:syntax=cpp.doxygen
+*/

src/libiir++/200_Filter.h

@@ -0,0 +1,170 @@
+/* libiir/src/libiir++/200_Filter.h
+ * 
+ *  Copyright: ©2014, Laurence Withers.
+ *  Author: Laurence Withers <l@lwithers.me.uk>
+ *  License: GPLv3
+ */
+
+
+
+/*! \defgroup libiirpp_filter Basic IIR filtering
+
+\ingroup libiirpp
+
+An instance of the \ref Filter class represents the structure and state of an
+%IIR filter. It may be built up with \ref Coeff instances or parsed from a
+user-supplied string (see \ref string_desc).
+
+*/
+
+
+
+/*! \brief %IIR filter coefficient set
+
+\ingroup libiirpp_filter
+
+An instance of this class contains both a copy of the set of coefficients used
+in each filter stage (of which there may be multiple), and all of the state
+(i.e.  previous inputs and outputs) associated with each stage.
+
+*/
+class Filter {
+public:
+    /*! \brief Construct from coefficient set.
+
+    \param coeff Coefficient set.
+
+    This constructor builds a filter chain with one stage. Additional stages may
+    be added via \ref chainCoeffSet(). This constructor will always result in a
+    valid object.
+
+    */
+    Filter(const Coeff& coeff);
+
+    /*! \brief Copy constructor.
+
+    \param other Filter to copy.
+
+    Copies the structure (but not the state) of the filter \a other. This
+    constructor will always result in a valid object.
+
+    */
+    Filter(const Filter& other);
+
+    /*! \brief Construct from filter description string.
+
+    \param desc Description string.
+
+    This constructor will parse the description string \a desc (see
+    \ref string_desc for details on what is accepted) and build up a filter
+    as described.
+
+    If the parsing step fails, then \ref valid() will return false and the
+    object must not be used (no other functions except the destructor will be
+    safe to call).
+
+    */
+    Filter(const char* desc);
+
+    /*! \brief Destructor. */
+    virtual ~Filter();
+
+
+
+    /*! \brief Determine whether object is valid.
+
+    \retval true if object is valid.
+    \retval false if object is not valid.
+
+    A Filter object might not be valid if it was constructed from a description
+    string that contained an error. In that case this function will return
+    \a false and the object should not be used (no other functions except the
+    destructor will be safe to call).
+
+    */
+    bool valid() const;
+
+
+
+    /*! \brief %Filter a sample.
+
+    \param x Sample.
+    \returns %Filter output.
+
+    Injects a single sample into the filter and returns the output value.
+
+    */
+    double f(double x);
+
+
+
+    /*! \brief Retrieve number of coefficient sets.
+
+    \returns Number of sets in use (≥1).
+
+    */
+    int numCoeffSet() const;
+
+    /*! \brief Retrieve coefficient set.
+
+    \param idx Index of set (0 ≤ \a idx < \ref numCoeffSet()).
+    \returns Reference to coefficient set.
+
+    */
+    const Coeff& getCoeffSet(int idx) const;
+
+    /*! \brief Add another set of coefficients to the filter chain.
+
+    \param coeff Coefficient set to add.
+
+    */
+    void chainCoeffSet(const Coeff& coeff);
+
+
+
+    /*! \brief Evaluate response of filter at given frequency.
+
+    \param freq Frequency (0 ≤ \a freq ≤ 1).
+    \returns Complex response of filter at given frequency.
+
+    Returns the steady-state response of the filter to a pure input signal of
+    the given frequency. \a freq is expressed as a fraction of the intended
+    sampling rate (with 0.5 being the Nyquist frequency).
+
+    */
+    std::complex<double> response(double freq) const;
+
+
+
+    /*! \brief Search for -3dB points.
+
+    \param gain Nominal gain of filter.
+    \returns List of frequency ratios at which filter magnitude is -3dB.
+
+    This function will perform a linear scan of the frequency response of
+    the filter between DC and the Nyquist frequency. It will then search
+    for any pair of points between which the magnitude transitions the -3dB
+    level. It then performs a few levels of binary search in order to provide
+    a more accurate answer.
+
+    Returned frequencies are expressed as a fraction of the sampling rate.
+
+    */
+    std::vector<double> responseMinus3dB(double gain) const;
+
+private:
+    // one entry per set of coefficients, in order of execution
+    std::vector<const Coeff*> coeff_;
+
+    // the C object has both state and its own copy of the coefficients
+    struct iir_filter_t* fi_;
+
+    // disallow assignment
+    Filter& operator=(const Filter& other);
+};
+
+
+
+/* options for text editors
+vim: expandtab:ts=4:sw=4:syntax=cpp.doxygen
+*/

src/libiir++/300_common_filters.cc

@@ -0,0 +1,58 @@
+/* libiir/src/libiir++/300_common_filters.cc
+ * 
+ *  Copyright: ©2014, Laurence Withers.
+ *  Author: Laurence Withers <l@lwithers.me.uk>
+ *  License: GPLv3
+ */
+
+
+
+Coeff
+ButterworthLowPass(int order, double gain, double corner)
+{
+    struct iir_coeff_t* cc = iir_butterworth_lowpass(order, gain, corner);
+    Coeff ret(cc);
+    iir_coeff_free(cc);
+    return ret;
+}
+
+
+
+Coeff
+ButterworthHighPass(int order, double gain, double corner)
+{
+    struct iir_coeff_t* cc = iir_butterworth_highpass(order, gain, corner);
+    Coeff ret(cc);
+    iir_coeff_free(cc);
+    return ret;
+}
+
+
+
+Coeff
+ButterworthBandPass(int order, double gain, double corner1, double corner2)
+{
+    struct iir_coeff_t* cc = iir_butterworth_bandpass(order, gain, corner1,
+        corner2);
+    Coeff ret(cc);
+    iir_coeff_free(cc);
+    return ret;
+}
+
+
+
+Coeff
+ButterworthBandStop(int order, double gain, double corner1, double corner2)
+{
+    struct iir_coeff_t* cc = iir_butterworth_bandstop(order, gain, corner1,
+        corner2);
+    Coeff ret(cc);
+    iir_coeff_free(cc);
+    return ret;
+}
+
+
+
+/* options for text editors
+vim: expandtab:ts=4:sw=4:syntax=cpp.doxygen
+*/

src/libiir++/300_common_filters.h

@@ -0,0 +1,144 @@
+/* libiir/src/libiir++/300_common_filters.h
+ * 
+ *  Copyright: ©2014, Laurence Withers.
+ *  Author: Laurence Withers <l@lwithers.me.uk>
+ *  License: GPLv3
+ */
+
+
+
+/*! \defgroup libiirpp_common_filters Common types of IIR filter
+
+\ingroup libiirpp
+
+Routines in this section can be used to build a coefficient set that matches a
+common %IIR filter type. These can be passed to the \ref Filter class to build
+up an actual instance of a filter.
+
+\note It is common to build higher-order filters by chaining together several
+    lower-order filters. The functions in this section do not do that directly,
+    but are intended to build the individual components of a chain. See
+    \ref Filter::Filter(const char*) for an interface which performs this
+    automatically.
+
+*/
+/*!@{*/
+
+
+
+/*! \brief nth-order Butterworth low-pass
+
+\param order Order of filter (≥1).
+\param gain Linear gain of filter.
+\param corner Corner frequency expressed as a fraction of Nyquist
+    (0 ≤ \a corner ≤ 1).
+
+\returns %IIR filter coefficient set.
+
+Uses the Exstrom labs code to compute the coefficients of an nth-order (param
+\a order) Butterworth-type low pass filter with gain \a gain and corner
+frequency \a corner.
+
+Note it is recommended to chain multiple filters together to build anything
+greater than a 4th-order filter. This function won't do that directly for you.
+\a gain will usually be set to be 1.0.
+
+The corner frequency corner is expressed as a fraction of the sampling
+frequency (which is of course not known by the %IIR code). It should lie between
+0 (0Hz) and 1 (the Nyquist frequency, or ½ the sampling frequency). 
+
+*/
+Coeff ButterworthLowPass(int order, double gain, double corner);
+
+
+
+/*! \brief nth-order Butterworth high-pass
+
+\param order Order of filter (≥1).
+\param gain Linear gain of filter.
+\param corner Corner frequency expressed as a fraction of Nyquist
+    (0 ≤ \a corner ≤ 1).
+
+\returns %IIR filter coefficient set.
+
+Uses the Exstrom labs code to compute the coefficients of an nth-order (param
+\a order) Butterworth-type high pass filter with gain \a gain and corner
+frequency \a corner.
+
+Note it is recommended to chain multiple filters together to build anything
+greater than a 4th-order filter. This function won't do that directly for you.
+\a gain will usually be set to be 1.0.
+
+The corner frequency corner is expressed as a fraction of the sampling
+frequency (which is of course not known by the %IIR code). It should lie between
+0 (0Hz) and 1 (the Nyquist frequency, or ½ the sampling frequency). 
+
+*/
+Coeff ButterworthHighPass(int order, double gain, double corner);
+
+
+
+/*! \brief nth-order Butterworth band-pass
+
+\param order Order of filter (≥1).
+\param gain Linear gain of filter.
+\param corner1 Low corner frequency expressed as a fraction of Nyquist
+    (0 ≤ \a corner ≤ 1).
+\param corner2 High corner frequency expressed as a fraction of Nyquist
+    (0 ≤ \a corner ≤ 1).
+
+\returns %IIR filter coefficient set.
+
+Uses the Exstrom labs code to compute the coefficients of an nth-order (param
+\a order) Butterworth-type band pass filter with gain \a gain and corner
+frequencies \a corner1 and \a corner2.
+
+Note it is recommended to chain multiple filters together to build anything
+greater than a 4th-order filter. This function won't do that directly for you.
+\a gain will usually be set to be 1.0.
+
+The corner frequencies \a corner1 and \a corner2 are expressed as a fraction of
+the sampling frequency (which is of course not known by the %IIR code). They
+should lie between 0 (0Hz) and 1 (the Nyquist frequency, or ½ the sampling
+frequency), and \a corner2 should be greater than \a corner1.
+
+*/
+Coeff ButterworthBandPass(int order, double gain, double corner1,
+    double corner2);
+
+
+
+/*! \brief nth-order Butterworth band-stop
+
+\param order Order of filter (≥1).
+\param gain Linear gain of filter.
+\param corner1 Low corner frequency expressed as a fraction of Nyquist
+    (0 ≤ \a corner ≤ 1).
+\param corner2 High corner frequency expressed as a fraction of Nyquist
+    (0 ≤ \a corner ≤ 1).
+
+\returns %IIR filter coefficient set.
+
+Uses the Exstrom labs code to compute the coefficients of an nth-order (param
+\a order) Butterworth-type band stop filter with gain \a gain and corner
+frequencies \a corner1 and \a corner2.
+
+Note it is recommended to chain multiple filters together to build anything
+greater than a 4th-order filter. This function won't do that directly for you.
+\a gain will usually be set to be 1.0.
+
+The corner frequencies \a corner1 and \a corner2 are expressed as a fraction of
+the sampling frequency (which is of course not known by the %IIR code). They
+should lie between 0 (0Hz) and 1 (the Nyquist frequency, or ½ the sampling
+frequency), and \a corner2 should be greater than \a corner1.
+
+*/
+Coeff ButterworthBandStop(int order, double gain, double corner1,
+    double corner2);
+
+
+
+/*!@}*/
+/* options for text editors
+vim: expandtab:ts=4:sw=4:syntax=cpp.doxygen
+*/

src/libiir++/999_BottomHeader.h

@@ -0,0 +1,14 @@
+/* libiir/src/libiir++/999_BottomHeader.h
+ *
+ *  Copyright: ©2014, Laurence Withers.
+ *  Author: Laurence Withers <l@lwithers.me.uk>
+ *  License: GPLv3
+ */
+
+} // end namespace IIR
+
+#endif
+
+/* options for text editors
+vim: expandtab:ts=4:sw=4:syntax=cpp.doxygen
+*/

src/libiir++/999_BottomSource.cc

@@ -0,0 +1,12 @@
+/* libiir/src/libiir++/999_BottomSource.cpp
+ *
+ *  Copyright: ©2014, Laurence Withers.
+ *  Author: Laurence Withers <l@lwithers.me.uk>
+ *  License: GPLv3
+ */
+
+} // end namspace IIR
+
+/* options for text editors
+vim: expandtab:ts=4:sw=4
+*/

src/libiir++/build.default

@@ -0,0 +1 @@
+source src/libiir++/build.lib

src/libiir++/build.install

@@ -0,0 +1 @@
+source src/libiir++/build.install-lib

src/libiir++/build.install-lib

@@ -0,0 +1,35 @@
+build_target libiir++
+
+# make paths (this is for Gentoo in particular)
+build_dir_tree "${LIBDIR}" || return 1
+build_dir_tree "${BINDIR}" || return 1
+build_dir_tree "${INCLUDEDIR}" || return 1
+
+# install library
+echo "Installing libraries into '${LIBDIR}'"
+install_file ${libiirpp} ${LIBDIR} 0755 || return 1
+BASE="${libiirpp_BASE}.so"
+MAJOR="${BASE}.${SOMAJOR}"
+MICRO="${MAJOR}.${SOMICRO}"
+install_symlink "${BASE}" "${MICRO}" "${LIBDIR}"
+
+# install header
+echo "Installing header file '${libiirpp_HEADER}' into ${INCLUDEDIR}"
+install_header ${libiirpp_HEADER} ${INCLUDEDIR} 0644 || return 1
+
+# install config script
+echo "Installing config script into ${PKGCONFDIR}"
+CONFFILE="${INSTALL_PREFIX}${BINDIR}/libiir++-config"
+
+do_cmd rm -f "${CONFFILE}"
+do_cmd_redir "${CONFFILE}" sed \
+    -e "s,@VERSION@,${VERSION}," \
+    -e "s,@LIB_DIR@,${LIBDIR}," \
+    -e "s,@INCLUDE_DIR@,${INCLUDEDIR}," \
+    -e "s,@DEP_CFLAGS@,${libiirpp_DEP_CFLAGS}," \
+    -e "s,@DEP_LIBS@,${libiirpp_DEP_LIBS}," \
+    src/libiir++/config-script
+do_cmd chmod 0644 "${CONFFILE}"
+print_success "Done"
+
+# vim: syntax=sh:expandtab:ts=4:sw=4

src/libiir++/build.lib

@@ -0,0 +1,58 @@
+# These are external variables, and shouldn't clash with anything else
+#  libiirpp
+#  libiirpp_BUILT
+#  libiirpp_HEADER
+#  libiirpp_BASE
+#  libiirpp_DEP_CFLAGS
+#  libiirpp_DEP_LIBS
+
+build_target libiir
+
+if [ -z ${libiirpp_BUILT} ]
+then
+    libiirpp_BASE=libiir++
+    source src/libiir++/soversion
+
+    libiirpp="obj/${libiirpp_BASE}.so.${SOMAJOR}.${SOMICRO}"
+    libiirpp_DEP_CFLAGS=""
+    libiirpp_DEP_LIBS="${libiir}"
+    SO_EXTRA="${libiirpp_DEP_CFLAGS} ${libiirpp_DEP_LIBS} -lstdc++ -lc \
+        -D_GNU_SOURCE"
+
+    echo "Building library ${libiirpp}..."
+
+    do_cmd source src/libiir++/build.monolithic || return 1
+
+    MODIFIED=0
+    for test in ${MONOLITHIC_TESTS} ${HDR} ${SRC}
+    do
+        if [ ${test} -nt ${libiirpp} ]
+        then
+            MODIFIED=1
+            break
+        fi
+    done
+
+    if [ ${MODIFIED} -ne 0 ]
+    then
+        echo " Compiling"
+
+        SONAME="${libiirpp_BASE}.so.${SOMAJOR}"
+        do_cmd ${CXX} ${CFLAGS} -Iobj -shared -fpic -o "${libiirpp}" \
+            -Wl,-soname,${SONAME} \
+            ${SRC} ${SO_EXTRA} || return 1
+
+        # make tests work
+        do_cmd ln -sf $(basename ${libiirpp}) obj/${SONAME} || return 1
+
+        print_success "Library built"
+    else
+        print_success "Library up to date"
+    fi
+
+    libiirpp_BUILT=1
+    libiirpp_HEADER=${HDR}
+
+fi
+
+# vim: syntax=sh:expandtab:ts=4:sw=4

src/libiir++/build.monolithic

@@ -0,0 +1,21 @@
+# These are external variables, and shouldn't clash with anything else
+#  libiirpp_MONOLITHIC
+
+SRC="obj/libiir++.cpp"
+HDR="obj/iir++.h"
+
+MONOLITHIC_TESTS="src/libiir++/build.lib src/libiir++/build.monolithic"
+
+if [ -z "${libiirpp_MONOLITHIC}" ]
+then
+    MONOLITHIC_SOURCE="$(find src/libiir++/ -name '*.h' | sort)"
+    make_monolithic ${HDR} Ch || return 1
+
+    MONOLITHIC_SOURCE="$(find src/libiir++/ -name '*.cc' | sort)"
+    make_monolithic ${SRC} C || return 1
+
+    libiirpp_MONOLITHIC=1
+    MONOLITHIC_DOC="${MONOLITHIC_DOC} ${HDR}"
+fi
+
+# vim: syntax=sh:expandtab:ts=4:sw=4

src/libiir++/config-script

@@ -0,0 +1,97 @@
+#!/bin/bash
+# libiir/src/libiir++/config-script
+#
+# libiir++-config template. Variables are finalised at install time.
+#
+dep_cflags="@DEP_CFLAGS@"
+dep_libs="@DEP_LIBS@"
+include_dir="@INCLUDE_DIR@"
+include_dir_set="no"
+lib_dir="@LIB_DIR@"
+lib_dir_set="no"
+
+
+
+usage() {
+    cat <<EOF
+Usage: libiir++-config [options]
+Options:
+        [--version]
+        [--libs]
+        [--libdir[=DIR]]
+        [--cflags]
+        [--includedir[=DIR]]
+EOF
+    exit $1
+}
+
+
+
+[ $# -eq 0 ] && usage 1 1>&2
+
+
+
+while [ $# -gt 0 ]
+do
+    case "$1" in
+        -*=*) 
+            optarg="$(echo "$1" | sed 's/[-_a-zA-Z0-9]*=//')"
+        ;;
+
+        *)
+            optarg=""
+        ;;
+    esac
+
+    case "$1" in
+        --libdir=*)
+            lib_dir="${optarg}"
+            lib_dir_set="yes"
+        ;;
+
+        --libdir)
+            echo_lib_dir="yes"
+        ;;
+
+        --includedir=*)
+            include_dir="${optarg}"
+            include_dir_set="yes"
+        ;;
+
+        --includedir)
+            echo_include_dir="yes"
+        ;;
+
+        --version)
+            echo "@VERSION@"
+            exit 0
+        ;;
+
+        --cflags)
+            [ "${include_dir}" != "/usr/include" ] && includes="-I${include_dir}"
+            echo_cflags="yes"
+        ;;
+
+        --libs)
+            echo_libs="yes"
+        ;;
+
+        *)
+            usage 1 1>&2
+        ;;
+    esac
+
+    shift
+done
+
+
+
+[ "${echo_prefix}" == "yes" ] && echo "${prefix}"
+[ "${echo_exec_prefix}" == "yes" ] && echo "${exec_prefix}"
+[ "${echo_cflags}" == "yes" ] && echo "${dep_cflags} ${includes}"
+[ "${echo_libs}" == "yes" ] && echo "${dep_libs} -L${lib_dir} -liir++"
+true
+
+
+
+# kate: 

src/libiir++/soversion

@@ -0,0 +1,15 @@
+# libiir/src/libiir++/soversion
+#
+#  Copyright: ©2014, Laurence Withers.
+#  Author: Laurence Withers <l@lwithers.me.uk>
+#  License: GPLv3
+#
+
+
+
+# SOMAJOR is included in the library's soname, and needs to be bumped
+# after a binary-incompatible release. It is a single integer.
+SOMAJOR=0
+
+# SOMICRO is bumped every time there is a binary-compatible release.
+SOMICRO=0

src/libiir/000_TopHeader.h

@@ -1,6 +1,6 @@
 /* libiir/src/libiir/000_TopHeader.h
  * 
- *  Copyright: ©2010, Laurence Withers.
+ *  Copyright: ©2010–2011, Laurence Withers.
  *  Author: Laurence Withers <l@lwithers.me.uk>
  *  License: GPLv3
 */
@@ -9,8 +9,14 @@
 #define HEADER_libiir
 
 /* standard includes, or includes needed for type declarations */
+#include <complex.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/*! \defgroup libiir C library */
 
 /* options for text editors
-kate: replace-trailing-space-save true; space-indent true; tab-width 4;
 vim: expandtab:ts=4:sw=4:syntax=c.doxygen
 */

src/libiir/000_TopSource.c

@@ -1,6 +1,6 @@
 /* libiir/src/libiir/000_TopSource.c
  * 
- *  Copyright: ©2010, Laurence Withers.
+ *  Copyright: ©2010–2014, Laurence Withers.
  *  Author: Laurence Withers <l@lwithers.me.uk>
  *  License: GPLv3
 */
@@ -15,6 +15,5 @@
 #include <string.h>
 
 /* options for text editors
-kate: replace-trailing-space-save true; space-indent true; tab-width 4;
 vim: expandtab:ts=4:sw=4
 */

src/libiir/100_coeff.c

@@ -0,0 +1,113 @@
+/* libiir/src/libiir/100_coeff.c
+ * 
+ *  Copyright: ©2010–2014, Laurence Withers.
+ *  Author: Laurence Withers <l@lwithers.me.uk>
+ *  License: GPLv3
+*/
+
+
+
+/* struct iir_coeff_t
+ *  Holds a general IIR filter (i.e. the set of coefficients that define it).
+ *  nc >= 1 and nd >= 1.
+ */
+struct iir_coeff_t {
+    int nc, nd;
+    double* c, * d;
+};
+
+
+
+/* iir_coeff_new()
+ *  Allocates a new set of coefficient objects.
+ */
+struct iir_coeff_t*
+iir_coeff_new(int nc, const double* c, int nd, const double* d)
+{
+    struct iir_coeff_t* coeff;
+
+    if(nc < 1 || nd < 1) {
+        errno = EINVAL;
+        return 0;
+    }
+
+    coeff = malloc(sizeof(struct iir_coeff_t));
+    coeff->nc = nc;
+    coeff->nd = nd;
+    coeff->c = malloc(sizeof(double) * nc);
+    coeff->d = malloc(sizeof(double) * nd);
+    memcpy(coeff->c, c, sizeof(double) * nc);
+    memcpy(coeff->d, d, sizeof(double) * nd);
+
+    return coeff;
+}
+
+
+
+/* iir_coeff_copy()
+ *  Deep copy of coefficient object.
+ */
+struct iir_coeff_t*
+iir_coeff_copy(const struct iir_coeff_t* other)
+{
+    struct iir_coeff_t* coeff;
+
+    coeff = malloc(sizeof(struct iir_coeff_t));
+    coeff->nc = other->nc;
+    coeff->nd = other->nd;
+    coeff->c = malloc(sizeof(double) * coeff->nc);
+    coeff->d = malloc(sizeof(double) * coeff->nd);
+    memcpy(coeff->c, other->c, sizeof(double) * coeff->nc);
+    memcpy(coeff->d, other->d, sizeof(double) * coeff->nd);
+
+    return coeff;
+}
+
+
+
+/* iir_coeff_free()
+ *  Frees memory associated with ‘coeff’.
+ */
+void
+iir_coeff_free(struct iir_coeff_t* coeff)
+{
+    if(!coeff) return;
+    free(coeff->c);
+    free(coeff->d);
+    free(coeff);
+}
+
+
+
+/* iir_coeff_get_*()
+ *  Functions for querying the coefficients that make up an IIR filter.
+ */
+int
+iir_coeff_get_nc(const struct iir_coeff_t* coeff)
+{
+    return coeff->nc;
+}
+
+int
+iir_coeff_get_nd(const struct iir_coeff_t* coeff)
+{
+    return coeff->nd;
+}
+
+double
+iir_coeff_get_c(const struct iir_coeff_t* coeff, int idx)
+{
+    return coeff->c[idx];
+}
+
+double
+iir_coeff_get_d(const struct iir_coeff_t* coeff, int idx)
+{
+    return coeff->d[idx];
+}
+
+
+
+/* options for text editors
+vim: expandtab:ts=4:sw=4
+*/

src/libiir/100_coeff.h

@@ -0,0 +1,126 @@
+/* libiir/src/libiir/100_coeff.h
+ * 
+ *  Copyright: ©2010–2014, Laurence Withers.
+ *  Author: Laurence Withers <l@lwithers.me.uk>
+ *  License: GPLv3
+*/
+
+
+
+/*! \defgroup libiir_coeff Coefficient handling
+
+\ingroup libiir
+
+The functions in this module present an interface for managing the set of
+coefficients of a single %IIR filter instance. See \ref iir_structure for
+the definition of a coefficient set. The coefficient sets form the
+basic building block of an %IIR filter instance (see \ref libiir_filter),
+which may chain several %IIR filters in sequence.
+
+*/
+/*!@{*/
+
+
+
+/* opaque structure */
+struct iir_coeff_t;
+
+
+
+/*! \brief Create %IIR coefficient set
+
+\param nc Number of \a c coefficients.
+\param c Array of \a c coefficients.
+\param nd Number of \a d coefficients.
+\param d Array of \a d coefficients.
+\returns Pointer to new general %IIR filter object.
+
+This function creates a new set of %IIR filter coefficients which may be used to
+create filter instances through \ref iir_filter_new() or chained on to existing
+instances through \ref iir_filter().
+
+See \ref iir_structure for a full explanation of the parameters.
+
+*/
+struct iir_coeff_t* iir_coeff_new(int nc, const double* c, int nd,
+    const double* d)
+#ifndef DOXYGEN
+    __attribute__((malloc,nonnull))
+#endif
+;
+
+
+
+/*! \brief Create copy of a set of %IIR coefficients
+
+\param other Set of %IIR filter coefficients to copy.
+\returns Pointer to new general %IIR filter object.
+
+Performs a deep copy of the set of %IIR coefficients contained within \a other.
+
+*/
+struct iir_coeff_t* iir_coeff_copy(const struct iir_coeff_t* other);
+
+
+
+/*! \brief Free %IIR coefficient set
+
+\param coeff Pointer to %IIR coefficient object. May be 0.
+
+Frees a set of %IIR filter coefficients previously allocated through
+\ref iir_coeff_new(). Can be called on a null pointer without consequences.
+Note that \ref iir_filter_new() and \ref iir_filter_chain() actually store a
+copy of the coefficients, so it is possible to free \a coeff even if existing
+filters are still using its coefficient values.
+
+*/
+void iir_coeff_free(struct iir_coeff_t* coeff);
+
+
+
+/*! \brief Query number of C coefficients.
+
+\param coeff Pointer to %IIR coefficient object.
+\returns Number of C coefficients (≥1).
+
+*/
+int iir_coeff_get_nc(const struct iir_coeff_t* coeff);
+
+
+
+/*! \brief Query number of D coefficients.
+
+\param coeff Pointer to %IIR coefficient object.
+\returns Number of D coefficients (≥1).
+
+*/
+int iir_coeff_get_nd(const struct iir_coeff_t* coeff);
+
+
+
+/*! \brief Get value of C coefficient.
+
+\param coeff Pointer to %IIR coefficient object.
+\param idx Index of coefficient (≥0, < \ref iir_coeff_get_nc()).
+\returns C coefficient value.
+
+*/
+double iir_coeff_get_c(const struct iir_coeff_t* coeff, int idx);
+
+
+
+/*! \brief Get value of D coefficient.
+
+\param coeff Pointer to %IIR coefficient object.
+\param idx Index of coefficient (≥0, < \ref iir_coeff_get_nd()).
+\returns D coefficient value.
+
+*/
+double iir_coeff_get_d(const struct iir_coeff_t* coeff, int idx);
+
+
+
+/*!@}*/
+/* options for text editors
+vim: expandtab:ts=4:sw=4:syntax=c.doxygen
+*/

src/libiir/200_filter.c

@@ -1,63 +1,12 @@
-/* libiir/src/libiir/200_iir.c
+/* libiir/src/libiir/200_filter.c
  * 
- *  Copyright: ©2010, Laurence Withers.
+ *  Copyright: ©2010–2014, Laurence Withers.
  *  Author: Laurence Withers <l@lwithers.me.uk>
  *  License: GPLv3
 */
 
 
 
-/* struct iir_coeff_t
- *  Holds a general IIR filter (i.e. the set of coefficients that define it).
- *  nc >= 1 and nd >= 1.
- */
-struct iir_coeff_t {
-    int nc, nd;
-    double* c, * d;
-};
-
-
-
-/* iir_coeff_new()
- *  Allocates a new set of coefficient objects.
- */
-struct iir_coeff_t*
-iir_coeff_new(int nc, double* c, int nd, double* d)
-{
-    struct iir_coeff_t* coeff;
-
-    if(nc < 1 || nd < 1) {
-        errno = EINVAL;
-        return 0;
-    }
-
-    coeff = malloc(sizeof(struct iir_coeff_t));
-    coeff->nc = nc;
-    coeff->nd = nd;
-    coeff->c = malloc(sizeof(double) * nc);
-    coeff->d = malloc(sizeof(double) * nd);
-    memcpy(coeff->c, c, sizeof(double) * nc);
-    memcpy(coeff->d, d, sizeof(double) * nd);
-
-    return coeff;
-}
-
-
-
-/* iir_coeff_free()
- *  Frees memory associated with ‘coeff’.
- */
-void
-iir_coeff_free(struct iir_coeff_t* coeff)
-{
-    if(!coeff) return;
-    free(coeff->c);
-    free(coeff->d);
-    free(coeff);
-}
-
-
-
 /* struct iir_filter_t
  *  An instantiated IIR filter. This is actually a linked list node, so that we
  *  can create chains of filters. It also has a copy of the coefficients so that
@@ -104,7 +53,6 @@ iir_filter_new(const struct iir_coeff_t* coeff)
     /* allocate space for state */
     fi->x = malloc(sizeof(double) * fi->nc);
     fi->y = malloc(sizeof(double) * fi->nd);
-    memset(fi->y, 0, sizeof(double) * fi->nd);
 
     return fi;
 }
@@ -150,7 +98,8 @@ iir_filter_chain(struct iir_filter_t* fi, const struct iir_coeff_t* coeff)
 /* iir_filter_copy()
  *  Performs a deep copy of a filter instance chain.
  */
-struct iir_filter_t* iir_filter_copy(const struct iir_filter_t* fi, int state)
+struct iir_filter_t*
+iir_filter_copy(const struct iir_filter_t* fi, int state)
 {
     struct iir_filter_t* head = 0, * tail = 0, * copy;
 
@@ -166,14 +115,13 @@ struct iir_filter_t* iir_filter_copy(const struct iir_filter_t* fi, int state)
         memcpy(copy->c, fi->c, sizeof(double) * fi->nc);
         memcpy(copy->d, fi->d, sizeof(double) * fi->nd);
 
-        if(state) {
+        if(state && fi->ready) {
             copy->ready = 1;
             memcpy(copy->x, fi->x, sizeof(double) * fi->nc);
             memcpy(copy->y, fi->y, sizeof(double) * fi->nd);
             copy->xpos = fi->xpos;
             copy->ypos = fi->ypos;
         } else {
-            memset(copy->y, 0, sizeof(double) * fi->nd);
             copy->ready = copy->xpos = copy->ypos = 0;
         }
 
@@ -203,15 +151,42 @@ iir_get_xy(const double* xy, int pos, int max, int step)
     return xy[pos];
 }
 
+static double
+iir_dc_gain(const struct iir_filter_t* fi)
+{
+    int i;
+    double sum, gain;
+
+    /* If the forward path DC gain is less than a threshold then we assume ‘fi’
+     * has high pass characteristics.
+     */
+    sum = 0;
+    for(i = 0; i < fi->nc; ++i) sum += fi->c[i];
+    gain = sum;
+
+    sum = 0;
+    for(i = 0; i < fi->nd; ++i) sum += fi->d[i];
+    if(fabs(sum - 1) < 0.001) {
+        /* ergh, looks unstable at DC; pretend DC gain is 0 */
+        return 0;
+    }
+
+    gain /= (1 + sum);
+    return gain;
+}
+
 double
 iir_filter(struct iir_filter_t* fi, double samp)
 {
     int i;
+    double gain;
 
     while(fi) {
         if(!fi->ready) {
             /* initial conditions */
+            gain = iir_dc_gain(fi);
             for(i = 0; i < fi->nc; ++i) fi->x[i] = samp;
+            for(i = 0; i < fi->nd; ++i) fi->y[i] = samp * gain;
             fi->ready = 1;
         }
 
@@ -242,7 +217,25 @@ iir_filter(struct iir_filter_t* fi, double samp)
 
 
 
+int
+iir_filter_coeff_sets(const struct iir_filter_t* fi)
+{
+    int count;
+    for(count = 0; fi; ++count) fi = fi->next;
+    return count;
+}
+
+
+
+struct iir_coeff_t*
+iir_filter_coeff_set(const struct iir_filter_t* fi, int idx)
+{
+    while(idx--) fi = fi->next;
+    return iir_coeff_new(fi->nc, fi->c, fi->nd, fi->d);
+}
+
+
+
 /* options for text editors
-kate: replace-trailing-space-save true; space-indent true; tab-width 4;
 vim: expandtab:ts=4:sw=4
 */

src/libiir/200_filter.h

@@ -1,23 +1,25 @@
-/* libiir/src/libiir/200_iir.h
+/* libiir/src/libiir/200_filter.h
  * 
- *  Copyright: ©2010, Laurence Withers.
+ *  Copyright: ©2010–2014, Laurence Withers.
  *  Author: Laurence Withers <l@lwithers.me.uk>
  *  License: GPLv3
 */
 
 
 
-/*! \defgroup iir Basic IIR filtering
+/*! \defgroup libiir_filter Basic IIR filtering
 
-The functions in this module present a basic interface for representing IIR
-filters, creating instances of (possibly chained) IIR filters, and filtering an
+\ingroup libiir
+
+The functions in this module present a basic interface for representing %IIR
+filters, creating instances of (possibly chained) %IIR filters, and filtering an
 input sample.
 
-A general IIR filter consists of a set of coefficients, and may be created
+A general %IIR filter consists of a set of coefficients, and may be created
 through \ref iir_coeff_new(). The filter object (the opaque <code>struct
 iir_coeff_t</code>) is then used to instantiate specific filters (the opaque
 <code>struct iir_filter_t</code>) through \ref iir_filter_new(). Each filter
-instance may have an arbitrary further number of IIR filters chained on to it
+instance may have an arbitrary further number of %IIR filters chained on to it
 through \ref iir_filter_chain(). The filter processes one sample at a time
 through \ref iir_filter().
 
@@ -26,65 +28,22 @@ through \ref iir_filter().
 
 
 
-/* opaque structure */
-struct iir_coeff_t;
-
-
-
-/*! \brief Create general IIR filter
-
-\param nc Number of \a c coefficients.
-\param c Array of \a c coefficients.
-\param nd Number of \a d coefficients.
-\param d Array of \a d coefficients.
-\returns Pointer to new general IIR filter object.
-
-This function creates a new general IIR filter object which may be used to
-create filter instances through \ref iir_filter_new() or chained on to existing
-instances through \ref iir_filter().
-
-See \ref iir_structure for a full explanation of the parameters.
-
-*/
-struct iir_coeff_t* iir_coeff_new(int nc, double* c, int nd, double* d)
-#ifndef DOXYGEN
-    __attribute__((malloc,nonnull))
-#endif
-;
-
-
-
-/*! \brief Free general IIR filter
-
-\param coeff Pointer to IIR filter object. May be 0.
-
-Frees a set of IIR filter coefficients previously allocated through
-\ref iir_coeff_new(). Can be called on a null pointer without consequences.
-Note that \ref iir_filter_new() and \ref iir_filter_chain() actually store a
-copy of the coefficients, so it is possible to free \a coeff even if existing
-filters are still using its coefficient values.
-
-*/
-void iir_coeff_free(struct iir_coeff_t* coeff);
-
-
-
 /* opaque structure */
 struct iir_filter_t;
 
 
 
-/*! \brief Create IIR filter instance
+/*! \brief Create %IIR filter instance
 
 \param coeff Filter coefficients to use.
-\returns Pointer to new instance of IIR filter.
+\returns Pointer to new instance of %IIR filter.
 
-Creates a new instance of a general IIR filter. The set of coefficients \a coeff
+Creates a new instance of a general %IIR filter. The set of coefficients \a coeff
 is copied into the returned structure, meaning the coefficients can be freed
 after this function returns if they will not be needed again.
 
 The first sample passed through \ref iir_filter() will be used to set initial
-conditions.
+conditions (see \ref iir_structure). 
 
 An arbitrary number of further filters may be chained on to the end of this
 instance through \ref iir_filter_chain().
@@ -98,11 +57,11 @@ struct iir_filter_t* iir_filter_new(const struct iir_coeff_t* coeff)
 
 
 
-/*! \brief Free IIR filter instance
+/*! \brief Free %IIR filter instance
 
 \param fi Filter object to free. May be 0.
 
-Frees a previously-allocated IIR filter instance. Can be called on a null
+Frees a previously-allocated %IIR filter instance. Can be called on a null
 pointer without consequences.
 
 */
@@ -110,12 +69,12 @@ void iir_filter_free(struct iir_filter_t* fi);
 
 
 
-/*! \brief Add a further IIR filter to a filter instance
+/*! \brief Add a further %IIR filter to a filter instance
 
 \param fi Filter instance to chain onto.
-\param coeff New IIR filter coefficients to add to chain.
+\param coeff New %IIR filter coefficients to add to chain.
 
-Extends an existing IIR filter by chaining a new set of coefficients onto the
+Extends an existing %IIR filter by chaining a new set of coefficients onto the
 end. This can be used for &gt;4th order Butterworth filters, for example. This
 copies the set of coefficients from \a coeff so the coefficients can be freed
 after this function returns if they are no longer required.
@@ -129,7 +88,7 @@ void iir_filter_chain(struct iir_filter_t* fi, const struct iir_coeff_t* coeff)
 
 
 
-/*! \brief Create a deep copy of an IIR filter instance
+/*! \brief Create a deep copy of an %IIR filter instance
 
 \param fi Filter instance to copy.
 \param state Non-zero to copy state as well.
@@ -162,8 +121,36 @@ double iir_filter(struct iir_filter_t* fi, double samp);
 
 
 
+/*! \brief Count number of coefficient sets in %IIR filter chain
+
+\param fi Filter object.
+\returns Number of coefficient sets in chain (≥1).
+
+Returns the number of discrete %IIR filter coefficient sets used in the filter
+object \a fi. There will always be at least one set.
+
+*/
+int iir_filter_coeff_sets(const struct iir_filter_t* fi);
+
+
+
+/*! \brief Get %IIR coefficient set from filter chain
+
+\param fi Filter object.
+\param idx Index of coefficient set (0 ≤ \a idx &lt; \ref iir_filter_coeff_sets()).
+\returns Newly-allocated coefficient set object.
+
+Extracts the %IIR filter coefficient set for the given step (\a idx) in the chain
+of filters in \a fi. Returns a newly-allocated filter coefficient set object
+which can be freed with \ref iir_coeff_free().
+
+*/
+struct iir_coeff_t* iir_filter_coeff_set(const struct iir_filter_t* fi,
+    int idx);
+
+
+
 /*!@}*/
 /* options for text editors
-kate: replace-trailing-space-save true; space-indent true; tab-width 4;
 vim: expandtab:ts=4:sw=4:syntax=c.doxygen
 */

src/libiir/300_common_filters.h

@@ -1,16 +1,18 @@
 /* libiir/src/libiir/300_common_filters.h
  * 
- *  Copyright: ©2010, Laurence Withers.
+ *  Copyright: ©2010–2014, Laurence Withers.
  *  Author: Laurence Withers <l@lwithers.me.uk>
  *  License: GPLv3
 */
 
 
 
-/*! \defgroup common_filters Common types of IIR filter
+/*! \defgroup libiir_common_filters Common types of IIR filter
 
-Functions to create coefficients for various common types of IIR filter. The
-coefficient structures which are returned may be used to instantiate IIR
+\ingroup libiir
+
+Functions to create coefficients for various common types of %IIR filter. The
+coefficient structures which are returned may be used to instantiate %IIR
 filters using \ref iir_filter_new().
 
 The Butterworth filter code comes from the Exstrom Labs LLC code available under
@@ -28,7 +30,7 @@ is a copy of the original code available in the top level of this project.
 \param gain Linear gain of filter.
 \param corner Corner frequency expressed as a fraction of Nyquist
     (0 ≤ \a corner ≤ 1)
-\returns Newly-allocated IIR filter coefficients.
+\returns Newly-allocated %IIR filter coefficients.
 
 Uses the Exstrom labs code to compute the coefficients of an nth-order (param
 \a order) Butterworth-type low pass filter with gain \a gain and corner
@@ -39,7 +41,7 @@ greater than a 4th-order filter. This function won't do that directly for you.
 \a gain will usually be set to be 1.0.
 
 The corner frequency \a corner is expressed as a fraction of the sampling
-frequency (which is of course not known by the IIR code). It should lie between
+frequency (which is of course not known by the %IIR code). It should lie between
 0 (0Hz) and 1 (the Nyquist frequency, or ½ the sampling frequency).
 
 */
@@ -59,7 +61,7 @@ struct iir_coeff_t* iir_butterworth_lowpass(int order,
 \param gain Linear gain of filter.
 \param corner Corner frequency expressed as a fraction of Nyquist
     (0 ≤ \a corner ≤ 1)
-\returns Newly-allocated IIR filter coefficients.
+\returns Newly-allocated %IIR filter coefficients.
 
 Uses the Exstrom labs code to compute the coefficients of an nth-order (param
 \a order) Butterworth-type high pass filter with gain \a gain and corner
@@ -70,7 +72,7 @@ greater than a 4th-order filter. This function won't do that directly for you.
 \a gain will usually be set to be 1.0.
 
 The corner frequency \a corner is expressed as a fraction of the sampling
-frequency (which is of course not known by the IIR code). It should lie between
+frequency (which is of course not known by the %IIR code). It should lie between
 0 (0Hz) and 1 (the Nyquist frequency, or ½ the sampling frequency).
 
 */
@@ -92,7 +94,7 @@ struct iir_coeff_t* iir_butterworth_highpass(int order,
     (0 ≤ \a c1 ≤ 1)
 \param c2 High corner frequency expressed as a fraction of Nyquist
     (0 ≤ \a c2 ≤ 1, and \a c1 &lt; \a c2)
-\returns Newly-allocated IIR filter coefficients.
+\returns Newly-allocated %IIR filter coefficients.
 
 Uses the Exstrom labs code to compute the coefficients of an nth-order (param
 \a order) Butterworth-type band pass filter with gain \a gain and corner
@@ -103,7 +105,7 @@ greater than a 4th-order filter. This function won't do that directly for you.
 \a gain will usually be set to be 1.0.
 
 The corner frequencies \a c1 and \a c2 are expressed as a fraction of the
-sampling frequency (which is of course not known by the IIR code). They should
+sampling frequency (which is of course not known by the %IIR code). They should
 lie between 0 (0Hz) and 1 (the Nyquist frequency, or ½ the sampling frequency),
 and \a c2 should be greater than \a c1.
 
@@ -127,7 +129,7 @@ struct iir_coeff_t* iir_butterworth_bandpass(int order,
     (0 ≤ \a c1 ≤ 1)
 \param c2 High corner frequency expressed as a fraction of Nyquist
     (0 ≤ \a c2 ≤ 1, and \a c1 &lt; \a c2)
-\returns Newly-allocated IIR filter coefficients.
+\returns Newly-allocated %IIR filter coefficients.
 
 Uses the Exstrom labs code to compute the coefficients of an nth-order (param
 \a order) Butterworth-type band stop filter with gain \a gain and corner
@@ -138,7 +140,7 @@ greater than a 4th-order filter. This function won't do that directly for you.
 \a gain will usually be set to be 1.0.
 
 The corner frequencies \a c1 and \a c2 are expressed as a fraction of the
-sampling frequency (which is of course not known by the IIR code). They should
+sampling frequency (which is of course not known by the %IIR code). They should
 lie between 0 (0Hz) and 1 (the Nyquist frequency, or ½ the sampling frequency),
 and \a c2 should be greater than \a c1.
 
@@ -156,6 +158,5 @@ struct iir_coeff_t* iir_butterworth_bandstop(int order,
 
 /*!@}*/
 /* options for text editors
-kate: replace-trailing-space-save true; space-indent true; tab-width 4;
 vim: expandtab:ts=4:sw=4:syntax=c.doxygen
 */

src/libiir/400_parser.c

@@ -1,6 +1,6 @@
 /* libiir/src/libiir/400_parser.c
  * 
- *  Copyright: ©2010, Laurence Withers.
+ *  Copyright: ©2010–2011, Laurence Withers.
  *  Author: Laurence Withers <l@lwithers.me.uk>
  *  License: GPLv3
 */
@@ -103,8 +103,13 @@ iir_parser_raw(struct iir_filter_t** fi, const char* desc)
  *  The maximum order in any single Butterworth-type filter. If the given order
  *  exceeds this, we split the resulting filter up into multiple sets of
  *  coefficients of this order or less.
+ *
+ *  This is required since IIR filters, even with double-precision floating
+ *  point numbers, are particularly susceptible to quantization errors. In
+ *  practice, anything with a steeper roll-off than 2nd order can go unstable
+ *  when presented with real input data.
  */
-#define IIR_PARSER_BW_MAX_ORDER         (4)
+#define IIR_PARSER_BW_MAX_ORDER         (2)
 
 
 
@@ -176,7 +181,7 @@ iir_parser_bw_aux2(struct iir_filter_t** fi,
         }
         iir_coeff_free(coeff);
         if(!order) return 0;
-        /* add a <4th order segment */
+        /* add a <2nd order segment */
     }
 
     coeff = bw(order, gain, corner);
@@ -214,7 +219,7 @@ iir_parser_bw_aux3(struct iir_filter_t** fi,
         }
         iir_coeff_free(coeff);
         if(!order) return 0;
-        /* add a <4th order segment */
+        /* add a <2nd order segment */
     }
 
     coeff = bw(order, gain, c1, c2);
@@ -382,7 +387,7 @@ iir_parse_n(const char* desc, int n)
 
     a = malloc(sizeof(struct iir_filter_t*) * n);
     a[0] = fi;
-    for(i = 1; i < n; ++i) a[n] = iir_filter_copy(fi, 0);
+    for(i = 1; i < n; ++i) a[i] = iir_filter_copy(fi, 0);
 
     return a;
 }
@@ -390,6 +395,5 @@ iir_parse_n(const char* desc, int n)
 
 
 /* options for text editors
-kate: replace-trailing-space-save true; space-indent true; tab-width 4;
 vim: expandtab:ts=4:sw=4:syntax=c.doxygen
 */

src/libiir/400_parser.h

@@ -1,17 +1,19 @@
 /* libiir/src/libiir/400_parser.h
  * 
- *  Copyright: ©2010, Laurence Withers.
+ *  Copyright: ©2010–2014, Laurence Withers.
  *  Author: Laurence Withers <l@lwithers.me.uk>
  *  License: GPLv3
 */
 
 
 
-/*! \defgroup parser Parser for user-specified IIR filters
+/*! \defgroup iir_parser Parser for user-specified IIR filters
 
-This is a high-level interface that can instantiate a set of IIR filters based
+\ingroup libiir
+
+This is a high-level interface that can instantiate a set of %IIR filters based
 on a user-specified, human-readable string. The intention of this interface is
-to allow IIR filters to be specified in configuration files so that they can be
+to allow %IIR filters to be specified in configuration files so that they can be
 easily modified by the user and easily understood/parsed by the system.
 
 See \ref string_desc for details on the string description format.
@@ -21,14 +23,14 @@ See \ref string_desc for details on the string description format.
 
 
 
-/*! \brief Instantiate an IIR filter based on a string description
+/*! \brief Instantiate an %IIR filter based on a string description
 
-\param desc IIR filter description.
-\returns Pointer to newly-allocated IIR filter instance.
+\param desc %IIR filter description.
+\returns Pointer to newly-allocated %IIR filter instance.
 \retval 0 on error.
 
-Parses the human-readable description of an IIR filter chain in \a desc,
-instantiating an IIR filter object to match. Returns the new filter. If \a desc
+Parses the human-readable description of an %IIR filter chain in \a desc,
+instantiating an %IIR filter object to match. Returns the new filter. If \a desc
 cannot be parsed correctly, returns 0 and sets \a errno to \c EINVAL.
 
 */
@@ -40,15 +42,15 @@ struct iir_filter_t* iir_parse(const char* desc)
 
 
 
-/*! \brief Instantiate a set of IIR filters based on a string description
+/*! \brief Instantiate a set of %IIR filters based on a string description
 
-\param desc IIR filter description.
+\param desc %IIR filter description.
 \param n Number of instances to allocate.
-\returns Pointer to array of \a n newly-allocated IIR filter instances.
+\returns Pointer to array of \a n newly-allocated %IIR filter instances.
 \retval 0 on error.
 
-Parses the human-readable description of an IIR filter chain in \a desc,
-instantiating a set of \a n identical IIR filter objects to match. Returns a
+Parses the human-readable description of an %IIR filter chain in \a desc,
+instantiating a set of \a n identical %IIR filter objects to match. Returns a
 pointer to an array of new filters. If \a desc cannot be parsed correctly,
 returns 0 and sets \a errno to \c EINVAL.
 
@@ -66,6 +68,5 @@ struct iir_filter_t** iir_parse_n(const char* desc, int n)
 
 /*!@}*/
 /* options for text editors
-kate: replace-trailing-space-save true; space-indent true; tab-width 4;
 vim: expandtab:ts=4:sw=4:syntax=c.doxygen
 */

src/libiir/500_response.c

@@ -0,0 +1,127 @@
+/* libiir/src/libiir/500_response.c
+ * 
+ *  Copyright: ©2011, Laurence Withers.
+ *  Author: Laurence Withers <l@lwithers.me.uk>
+ *  License: GPLv3
+*/
+
+
+
+static double complex
+iir_response_aux(int nc, const double* c, int nd, const double* d, double freq)
+{
+    int i;
+    double complex sum, resp, N;
+
+    /* Z transform; substitutions:
+     *
+     * z = exp(jω)
+     *   where j is the imaginary unit (written as I in the C code below)
+     *   and where ω = 2πf
+     *
+     * z^-n = exp(jωn)
+     *
+     * exp(jx) = cos(x) + j.sin(x)
+     *
+     * z^-n = cos(-ωn) + j.sin(-ωn)
+     *
+     * N = -ω = -2πf
+     */
+    N = -2 * M_PI * freq;
+
+    sum = 0;
+    for(i = 0; i < nc; ++i) {
+        sum += c[i] * (cos(i*N) + I * sin(i*N));
+    }
+
+    resp = sum;
+    sum = 1;
+    for(i = 1; i <= nd; ++i) {
+        sum += d[i - 1] * (cos(i*N) + I * sin(i*N));
+    }
+
+    return resp / sum;
+}
+
+
+
+double complex
+iir_response_c(const struct iir_coeff_t* coeff, double freq)
+{
+    return iir_response_aux(coeff->nc, coeff->c, coeff->nd, coeff->d, freq);
+}
+
+
+
+double complex
+iir_response(const struct iir_filter_t* fi, double freq)
+{
+    double complex resp;
+
+    for(resp = 1.0; fi; fi = fi->next) {
+        resp *= iir_response_aux(fi->nc, fi->c, fi->nd, fi->d, freq);
+    }
+    return resp;
+}
+
+
+
+void
+iir_response_minus_3dB_aux(const struct iir_filter_t* fi, double gain,
+    double** freqs, int* num_freqs, int* sz_freqs, double f1, double f2)
+{
+    double f3;
+    double mag;
+    int i;
+
+    for(i = 0; i < 20; ++i) {
+        f3 = (f1 + f2) / 2;
+        mag = cabs(iir_response(fi, f3)) / gain;
+        if(mag < 0.5) {
+            f1 = f3;
+        } else {
+            f2 = f3;
+        }
+    }
+
+    if(*num_freqs == *sz_freqs) {
+        *sz_freqs = (*sz_freqs) ? (*sz_freqs << 1) : 4;
+        *freqs = realloc(*freqs, sizeof(double) * (*sz_freqs));
+    }
+    (*freqs)[*num_freqs] = f3;
+    ++(*num_freqs);
+}
+
+void
+iir_response_minus_3dB(const struct iir_filter_t* fi, double gain,
+    double** freqs, int* num_freqs)
+{
+    int sz_freqs = 0, pos;
+    double mag, last_mag;
+
+    *freqs = 0;
+    *num_freqs = 0;
+
+    last_mag = cabs(iir_response(fi, 0)) / gain;
+    for(pos = 1; pos < 1000; ++pos) {
+        mag = cabs(iir_response(fi, 0.0005 * pos)) / gain;
+        if(mag < 0.5) {
+            if(last_mag >= 0.5) {
+                iir_response_minus_3dB_aux(fi, gain, freqs, num_freqs,
+                    &sz_freqs, 0.0005 * pos, 0.0005 * (pos - 1));
+            }
+        } else {
+            if(last_mag < 0.5) {
+                iir_response_minus_3dB_aux(fi, gain, freqs, num_freqs,
+                    &sz_freqs, 0.0005 * (pos - 1), 0.0005 * pos);
+            }
+        }
+        last_mag = mag;
+    }
+}
+
+
+
+/* options for text editors
+vim: expandtab:ts=4:sw=4
+*/

src/libiir/500_response.h

@@ -0,0 +1,77 @@
+/* libiir/src/libiir/500_response.h
+ * 
+ *  Copyright: ©2011–2014, Laurence Withers.
+ *  Author: Laurence Withers <l@lwithers.me.uk>
+ *  License: GPLv3
+*/
+
+
+
+/*! \defgroup iir_response Measuring IIR response
+
+\ingroup libiir
+
+These functions allow measuring the %IIR response (i.e. frequency and phase of
+the transfer function) at a given frequency. The response is a complex number;
+see \c complex(5) for details. The magnitude \c cabs(z) is the gain of the
+filter and the phase \c carg(z) its phase delay.
+
+*/
+/*!@{*/
+
+
+
+/*! \brief Find response of coefficients at a single frequency.
+
+\param coeff Set of coefficients.
+\param freq Frequency of interest, as a fraction of sample rate.
+\returns Response as complex number.
+
+*/
+double complex iir_response_c(const struct iir_coeff_t* coeff, double freq);
+
+
+
+/*! \brief Find response of filter at a single frequency.
+
+\param fi Filter to evaluate.
+\param freq Frequency of interest, as a fraction of sample rate.
+\returns Response as a complex number.
+
+Similar to \ref iir_response_c(), except that it works on an instantiated
+filter (and takes all parts of a chained filter into account). Does not alter
+the state of the filter in any way and calls may be interleaved with filter
+evaluation.
+
+*/
+double complex iir_response(const struct iir_filter_t* fi, double freq);
+
+
+
+/*! \brief Search for -3dB points.
+
+\param fi Filter to evaluate.
+\param gain Nominal gain of the filter.
+\param[out] freqs Set to point to array of frequencies
+\param[out] num_freqs Number of frequencies in \a freq.
+
+This function will perform a linear scan of the frequency response of the
+filter between DC and the Nyquist frequency. It will then search for any pair
+of points between which the magnitude transitions the -3dB level. It then
+performs a few levels of binary search in order to provide a more accurate
+answer.
+
+Returned frequencies are expressed as a fraction of the sampling rate. The
+function will allocate an array of \c double with \c malloc(3), so the result
+should later be passed to \c free(3).
+
+*/
+void iir_response_minus_3dB(const struct iir_filter_t* fi, double gain,
+    double** freqs, int* num_freqs);
+
+
+
+/*!@}*/
+/* options for text editors
+vim: expandtab:ts=4:sw=4:syntax=c.doxygen
+*/

src/libiir/999_BottomHeader.h

@@ -1,13 +1,16 @@
 /* libiir/src/libiir/999_BottomHeader.h
  * 
- *  Copyright: ©2010, Laurence Withers.
+ *  Copyright: ©2010–2014, Laurence Withers.
  *  Author: Laurence Withers <l@lwithers.me.uk>
  *  License: GPLv3
 */
 
+#ifdef __cplusplus
+}
+#endif
+
 #endif
 
 /* options for text editors
-kate: replace-trailing-space-save true; space-indent true; tab-width 4;
 vim: expandtab:ts=4:sw=4:syntax=c.doxygen
 */

src/libiir/soversion

@@ -12,4 +12,4 @@
 SOMAJOR=0
 
 # SOMICRO is bumped every time there is a binary-compatible release.
-SOMICRO=0
+SOMICRO=2

src/tests/plot_filter.c

@@ -1,6 +1,6 @@
 /* libiir/src/tests/plot_filter.c
  * 
- *  Copyright: ©2010, Laurence Withers.
+ *  Copyright: ©2010–2011, Laurence Withers.
  *  Author: Laurence Withers <l@lwithers.me.uk>
  *  License: GPLv3
 */
@@ -50,18 +50,7 @@ do_plot(const char* filter_desc, double samp_rat, const char* png_filename)
         "set output '%s'\n"
         "set multiplot layout 2,1 title \"Bode plot for filter '",
         png_filename);
-    ret = 0;
-    while(*filter_desc) {
-        if(isspace(*filter_desc)) {
-            if(!ret) {
-                ret = 1;
-                putc('\n', fp);
-            }
-        } else {
-            putc(*filter_desc, fp);
-        }
-        ++filter_desc;
-    }
+    fputs(strstr(filter_desc, "raw") ? "raw" : filter_desc, fp);
     fprintf(fp, "' at %fHz\"\n"
         "set grid\n"
         "set logscale\n"
@@ -209,7 +198,7 @@ main(int argc, char* argv[])
 
     /* process commandline arguments */
     if(argc == 2 && !strcmp(argv[1], "--print-summary")) {
-        fputs("Generates Bode plot for a filter.\n", stdout);
+        fputs("Generates Bode plot for a filter (simulation method).\n", stdout);
         return 0;
     }
 
@@ -266,6 +255,5 @@ main(int argc, char* argv[])
 
 
 /* options for text editors
-kate: replace-trailing-space-save true; space-indent true; tab-width 4;
 vim: expandtab:ts=4:sw=4
 */

src/tests/zplot_filter.c

@@ -0,0 +1,194 @@
+/* libiir/src/tests/zplot_filter.c
+ * 
+ *  Copyright: ©2010–2011, Laurence Withers.
+ *  Author: Laurence Withers <l@lwithers.me.uk>
+ *  License: GPLv3
+*/
+
+#include "iir.h"
+
+#include <math.h>
+#include <ctype.h>
+#include <errno.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <unistd.h>
+#include <complex.h>
+
+
+
+#define NPOINTS                         (1000)
+#define STEADY_STATE_CYCLES             (20)
+
+
+
+char* tmp_fname;
+void
+unlink_tmpfile(void)
+{
+    unlink(tmp_fname);
+}
+
+
+
+int
+do_plot(const char* filter_desc, double samp_rat, const char* png_filename)
+{
+    int fd, ret;
+    FILE* fp;
+    char cmd_file[] = "/tmp/libiir-plot_filter.cmd.XXXXXX",
+        cmd[200];
+
+    fd = mkstemp(cmd_file);
+    if(fd == -1) {
+        perror("mkstemp");
+        return -1;
+    }
+
+    fp = fdopen(fd, "w");
+    fprintf(fp, "set terminal png size 1000,1000\n"
+        "set output '%s'\n"
+        "set multiplot layout 2,1 title \"Bode plot for filter '",
+        png_filename);
+    fputs(strstr(filter_desc, "raw") ? "raw" : filter_desc, fp);
+    fprintf(fp, "' at %fHz\"\n"
+        "set grid\n"
+        "set logscale\n"
+        "set ytics add ('-3dB' %f)\n"
+        "set xlabel 'Frequency (Hz)'\n"
+        "set ylabel 'Gain'\n"
+        "plot '%s' using 1:2 notitle\n"
+        "unset logscale y\n"
+        "set yrange [-180:180]\n"
+        "set ytics -180,45,180\n"
+        "set ylabel 'Phase (degrees)'\n"
+        "plot '%s' using 1:3 notitle\n",
+        samp_rat, 
+        pow(10, -3.0/20),
+        tmp_fname,
+        tmp_fname);
+    fclose(fp);
+
+    snprintf(cmd, sizeof(cmd), "gnuplot %s", cmd_file);
+    ret = system(cmd);
+    unlink(cmd_file);
+
+    return ret;
+}
+
+
+
+double
+interp(int step, int max, double start, double end)
+{
+    return start + step * ((end - start) / (max - 1));
+}
+
+void
+calc_response(FILE* fp,
+    const struct iir_filter_t* fi,
+    double samp_rat,
+    double start_freq,
+    double end_freq)
+{
+    int step;
+    double freq;
+    double complex resp;
+
+    for(step = 0; step < NPOINTS; ++step) {
+        freq = exp(interp(step, NPOINTS, log(start_freq), log(end_freq)));
+        resp = iir_response(fi, freq / samp_rat);
+
+        fprintf(fp, "%e\t%e\t% 6.2f\n",
+            freq,
+            cabs(resp),
+            carg(resp) * 180 / M_PI);
+    }
+}
+
+
+
+int
+safe_strtod(const char* str, double* d)
+{
+    char* endp = 0;
+    errno = 0;
+    *d = strtod(str, &endp);
+    if(errno || !endp || *endp) return -1;
+    return 0;
+}
+
+
+
+int
+main(int argc, char* argv[])
+{
+    int fd;
+    double samp_rat, start_freq, end_freq;
+    FILE* fp;
+    struct iir_filter_t* fi;
+
+    /* process commandline arguments */
+    if(argc == 2 && !strcmp(argv[1], "--print-summary")) {
+        fputs("Generates Bode plot for a filter (Z transform method).\n", stdout);
+        return 0;
+    }
+
+    if(argc != 6) {
+        fputs("Usage: plot_filter 'filter_desc' samp_rat start_freq end_freq out.png\n",
+            stderr);
+        return 1;
+    }
+
+    fi = iir_parse(argv[1]);
+    if(!fi) {
+        fputs("Invalid filter description string.\n", stderr);
+        return 1;
+    }
+
+    if(safe_strtod(argv[2], &samp_rat) || samp_rat < 1e-6) {
+        fputs("Invalid sample rate. Positive float in Hz.\n", stderr);
+        return 1;
+    }
+
+    if(safe_strtod(argv[3], &start_freq) || start_freq < 1e-6) {
+        fputs("Invalid start frequency. Positive float in Hz.\n", stderr);
+        return 1;
+    }
+
+    if(safe_strtod(argv[4], &end_freq) || end_freq < 1e-6
+        || end_freq > samp_rat || end_freq < start_freq)
+    {
+        fputs("Invalid end frequency. Positive float in Hz, less than sample\n"
+            "rate, but greater than start frequency.\n", stderr);
+        return 1;
+    }
+
+    printf("DC response: %f\n", cabs(iir_response(fi, 0)));
+
+    /* create temporary file for results; gnuplot will use this */
+    tmp_fname = strdup("/tmp/libiir-plot_filter.data.XXXXXX");
+    fd = mkstemp(tmp_fname);
+    if(fd == -1) {
+        perror("mkstemp");
+        return 1;
+    }
+    atexit(unlink_tmpfile);
+
+    fp = fdopen(fd, "w");
+    calc_response(fp, fi, samp_rat, start_freq, end_freq);
+    fclose(fp);
+
+    /* clean up (for valgrind) */
+    iir_filter_free(fi);
+
+    /* draw the plot */
+    return do_plot(argv[1], samp_rat, argv[5]);
+}
+
+
+
+/* options for text editors
+vim: expandtab:ts=4:sw=4
+*/

version

@@ -10,7 +10,7 @@
 # VERSION contains the full version number of the library, which is 
 # expected to be in 'major.minor.micro' format.
 VERMAJOR=1
-VERMINOR=0
+VERMINOR=1
 VERMICRO=0
 
 # kate: replace-trailing-space-save true; space-indent true; tab-width 4;