From 724fcb6a0b03fd047758529ee212d59d13631f13 Mon Sep 17 00:00:00 2001
From: Laurence Withers <lwithers@guralp.com>
Date: Tue, 8 Jul 2014 08:27:50 +0000
Subject: [PATCH] Add -3dB search function

This commit adds a heuristic search function which finds -3dB points of
the specified filter. It first performs a 1000-point linear scan and
looks for transitions; it then performs a 20-point binary chop within
the transition band to better approximate the crossing point.
---
 src/libiir++/000_TopSource.cc |  1 +
 src/libiir++/200_Filter.cc    | 15 ++++++++++
 src/libiir++/200_Filter.h     | 18 +++++++++++
 src/libiir/500_response.c     | 56 +++++++++++++++++++++++++++++++++++
 src/libiir/500_response.h     | 23 ++++++++++++++
 5 files changed, 113 insertions(+)

diff --git a/src/libiir++/000_TopSource.cc b/src/libiir++/000_TopSource.cc
index 6daa3bf..f6b2040 100644
--- a/src/libiir++/000_TopSource.cc
+++ b/src/libiir++/000_TopSource.cc
@@ -16,6 +16,7 @@
 #include "iir++.h"
 
 // Below are all the includes used throughout the library.
+#include <stdlib.h>
 
 namespace IIR {
 
diff --git a/src/libiir++/200_Filter.cc b/src/libiir++/200_Filter.cc
index 027b0de..d0fd159 100644
--- a/src/libiir++/200_Filter.cc
+++ b/src/libiir++/200_Filter.cc
@@ -99,6 +99,21 @@ Filter::response(double freq) const
 
 
 
+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
 */
diff --git a/src/libiir++/200_Filter.h b/src/libiir++/200_Filter.h
index c073c4f..5a4d79f 100644
--- a/src/libiir++/200_Filter.h
+++ b/src/libiir++/200_Filter.h
@@ -134,6 +134,24 @@ public:
     */
     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_;
diff --git a/src/libiir/500_response.c b/src/libiir/500_response.c
index b19ee9b..8c6ad0d 100644
--- a/src/libiir/500_response.c
+++ b/src/libiir/500_response.c
@@ -66,6 +66,62 @@ iir_response(const struct iir_filter_t* fi, double freq)
 
 
 
+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
 */
diff --git a/src/libiir/500_response.h b/src/libiir/500_response.h
index 9572a28..885e1bf 100644
--- a/src/libiir/500_response.h
+++ b/src/libiir/500_response.h
@@ -48,6 +48,29 @@ 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