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.
This commit is contained in:
Laurence Withers 2014-07-08 08:27:50 +00:00
parent e5e391411f
commit 724fcb6a0b
5 changed files with 113 additions and 0 deletions

View File

@ -16,6 +16,7 @@
#include "iir++.h" #include "iir++.h"
// Below are all the includes used throughout the library. // Below are all the includes used throughout the library.
#include <stdlib.h>
namespace IIR { namespace IIR {

View File

@ -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 /* options for text editors
vim: expandtab:ts=4:sw=4:syntax=cpp.doxygen vim: expandtab:ts=4:sw=4:syntax=cpp.doxygen
*/ */

View File

@ -134,6 +134,24 @@ public:
*/ */
std::complex<double> response(double freq) const; 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: private:
// one entry per set of coefficients, in order of execution // one entry per set of coefficients, in order of execution
std::vector<const Coeff*> coeff_; std::vector<const Coeff*> coeff_;

View File

@ -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 /* options for text editors
vim: expandtab:ts=4:sw=4 vim: expandtab:ts=4:sw=4
*/ */

View File

@ -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 /* options for text editors
vim: expandtab:ts=4:sw=4:syntax=c.doxygen vim: expandtab:ts=4:sw=4:syntax=c.doxygen