From e154c135fc52966e70f205958b8d91300e855b61 Mon Sep 17 00:00:00 2001 From: Jean Pierre Cimalando Date: Mon, 15 Feb 2021 04:56:10 +0100 Subject: [PATCH] Add the filter designer from hiir --- cmake/SfizzDeps.cmake | 6 + .../hiir/hiir/PolyphaseIir2Designer.cpp | 444 ++++++++++++++++++ .../hiir/hiir/PolyphaseIir2Designer.h | 143 ++++++ 3 files changed, 593 insertions(+) create mode 100644 src/external/hiir/hiir/PolyphaseIir2Designer.cpp create mode 100644 src/external/hiir/hiir/PolyphaseIir2Designer.h diff --git a/cmake/SfizzDeps.cmake b/cmake/SfizzDeps.cmake index 864f1225..464d6186 100644 --- a/cmake/SfizzDeps.cmake +++ b/cmake/SfizzDeps.cmake @@ -105,6 +105,12 @@ add_library(sfizz_hiir INTERFACE) add_library(sfizz::hiir ALIAS sfizz_hiir) target_include_directories(sfizz_hiir INTERFACE "src/external/hiir") +# The hiir filter designer +add_library(sfizz_hiir_polyphase_iir2designer STATIC + "src/external/hiir/hiir/PolyphaseIir2Designer.cpp") +add_library(sfizz::hiir_polyphase_iir2designer ALIAS sfizz_hiir_polyphase_iir2designer) +target_link_libraries(sfizz_hiir_polyphase_iir2designer PUBLIC sfizz::hiir) + # The kissfft library add_library(sfizz_kissfft STATIC "src/external/kiss_fft/kiss_fft.c" diff --git a/src/external/hiir/hiir/PolyphaseIir2Designer.cpp b/src/external/hiir/hiir/PolyphaseIir2Designer.cpp new file mode 100644 index 00000000..d36883d6 --- /dev/null +++ b/src/external/hiir/hiir/PolyphaseIir2Designer.cpp @@ -0,0 +1,444 @@ +/***************************************************************************** + + PolyphaseIir2Designer.cpp + Author: Laurent de Soras, 2005 + +--- Legal stuff --- + +This program is free software. It comes without any warranty, to +the extent permitted by applicable law. You can redistribute it +and/or modify it under the terms of the Do What The Fuck You Want +To Public License, Version 2, as published by Sam Hocevar. See +http://sam.zoy.org/wtfpl/COPYING for more details. + +*Tab=3***********************************************************************/ + + + +#if defined (_MSC_VER) + #pragma warning (1 : 4130) // "'operator' : logical operation on address of string constant" + #pragma warning (1 : 4223) // "nonstandard extension used : non-lvalue array converted to pointer" + #pragma warning (1 : 4705) // "statement has no effect" + #pragma warning (1 : 4706) // "assignment within conditional expression" + #pragma warning (4 : 4786) // "identifier was truncated to '255' characters in the debug information" + #pragma warning (4 : 4800) // "forcing value to bool 'true' or 'false' (performance warning)" + #pragma warning (4 : 4355) // "'this' : used in base member initializer list" +#endif + + + +/*\\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ + +#include "hiir/def.h" +#include "hiir/fnc.h" +#include "hiir/PolyphaseIir2Designer.h" + +#include +#include + + + +namespace hiir +{ + + + +/*\\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ + + + +/* +============================================================================== +Name: compute_nbr_coefs_from_proto +Description: + Finds the minimum number of coefficients for a given filter specification +Input parameters: + - attenuation: stopband attenuation, dB. > 0. + - transition: normalized transition bandwith. Range ]0 ; 1/2[ +Returns: Number of coefficients, > 0 +Throws: Nothing +============================================================================== +*/ + +int PolyphaseIir2Designer::compute_nbr_coefs_from_proto (double attenuation, double transition) +{ + assert (attenuation > 0); + assert (transition > 0); + assert (transition < 0.5); + + double k; + double q; + compute_transition_param (k, q, transition); + const int order = compute_order (attenuation, q); + const int nbr_coefs = (order - 1) / 2; + + return nbr_coefs; +} + + + +/* +============================================================================== +Name: compute_atten_from_order_tbw +Description: + Compute the attenuation correspounding to a given number of coefficients + and the transition bandwith. +Input parameters: + - nbr_coefs: Number of desired coefficients. > 0. + - transition: normalized transition bandwith. Range ]0 ; 1/2[ +Returns: stopband attenuation, dB. > 0. +Throws: Nothing +============================================================================== +*/ + +double PolyphaseIir2Designer::compute_atten_from_order_tbw (int nbr_coefs, double transition) +{ + assert (nbr_coefs > 0); + assert (transition > 0); + assert (transition < 0.5); + + double k; + double q; + compute_transition_param (k, q, transition); + const int order = nbr_coefs * 2 + 1; + const double attenuation = compute_atten (q, order); + + return attenuation; +} + + + +/* +============================================================================== +Name: compute_coefs +Description: + Computes coefficients for a half-band polyphase IIR filter, function of a + given stopband gain / transition bandwidth specification. + Order is automatically calculated. +Input parameters: + - attenuation: stopband attenuation, dB. > 0. + - transition: normalized transition bandwith. Range ]0 ; 1/2[ +Output parameters: + - coef_arr: Coefficient list, must be large enough to store all the + coefficients. Filter order = nbr_coefs * 2 + 1 +Returns: number of coefficients +Throws: Nothing +============================================================================== +*/ + +int PolyphaseIir2Designer::compute_coefs (double coef_arr [], double attenuation, double transition) +{ + assert (attenuation > 0); + assert (transition > 0); + assert (transition < 0.5); + + double k; + double q; + compute_transition_param (k, q, transition); + + // Computes number of required coefficients + const int order = compute_order (attenuation, q); + const int nbr_coefs = (order - 1) / 2; + + // Coefficient calculation + for (int index = 0; index < nbr_coefs; ++index) + { + coef_arr [index] = compute_coef (index, k, q, order); + } + + return nbr_coefs; +} + + + +/* +============================================================================== +Name: compute_coefs_spec_order_tbw +Description: + Computes coefficients for a half-band polyphase IIR filter, function of a + given transition bandwidth and desired filter order. Bandstop attenuation + is set to the maximum value for these constraints. +Input parameters: + - nbr_coefs: Number of desired coefficients. > 0. + - transition: normalized transition bandwith. Range ]0 ; 1/2[ +Output parameters: + - coef_arr: Coefficient list, must be large enough to store all the + coefficients. +Throws: Nothing +============================================================================== +*/ + +void PolyphaseIir2Designer::compute_coefs_spec_order_tbw (double coef_arr [], int nbr_coefs, double transition) +{ + assert (nbr_coefs > 0); + assert (transition > 0); + assert (transition < 0.5); + + double k; + double q; + compute_transition_param (k, q, transition); + const int order = nbr_coefs * 2 + 1; + + // Coefficient calculation + for (int index = 0; index < nbr_coefs; ++index) + { + coef_arr [index] = compute_coef (index, k, q, order); + } +} + + + +/* +============================================================================== +Name: compute_phase_delay +Description: + Computes the phase delay introduced by a single filtering unit at a + specified frequency. + The delay is given for a constant sampling rate between input and output. +Input parameters: + - a: coefficient for the cell, [0 ; 1] + - f_fs: frequency relative to the sampling rate, [0 ; 0.5]. +Returns: + The phase delay in samples, >= 0. +Throws: Nothing +============================================================================== +*/ + +double PolyphaseIir2Designer::compute_phase_delay (double a, double f_fs) +{ + assert (a >= 0); + assert (a <= 1); + assert (f_fs >= 0); + assert (f_fs < 0.5); + + const double w = 2 * hiir::PI * f_fs; + const double c = cos (w); + const double s = sin (w); + const double x = a + c + a * (c * (a + c) + s * s); + const double y = a * a * s - s; + double ph = atan2 (y, x); + if (ph < 0) + { + ph += 2 * hiir::PI; + } + const double dly = ph / w; + + return dly; +} + + + +/* +============================================================================== +Name: compute_group_delay +Description: + Computes the group delay introduced by a single filtering unit at a + specified frequency. + The delay is given for a constant sampling rate between input and output. + To compute the group delay of a complete filter, add the group delays + of all the units in A0 (z). +Input parameters: + - a: coefficient for the cell, [0 ; 1] + - f_fs: frequency relative to the sampling rate, [0 ; 0.5]. + - ph_flag: set if filtering unit is used in pi/2-phaser mode, in the form + (a - z^-2) / (1 - az^-2) +Returns: + The group delay in samples, >= 0. +Throws: Nothing +============================================================================== +*/ + +double PolyphaseIir2Designer::compute_group_delay (double a, double f_fs, bool ph_flag) +{ + assert (a >= 0); + assert (a <= 1); + assert (f_fs >= 0); + assert (f_fs < 0.5); + + const double w = 2 * hiir::PI * f_fs; + const double a2 = a * a; + const double sig = (ph_flag) ? -2 : 2; + const double dly = 2 * (1 - a2) / (a2 + sig * a * cos (2 * w) + 1); + + return dly; +} + + + +/* +============================================================================== +Name: compute_group_delay +Description: + Computes the group delay introduced by a complete filter at a specified + frequency. + The delay is given for a constant sampling rate between input and output. +Input parameters: + - coef_arr: filter coefficient, as given by the designing functions + - nbr_coefs: Number of filter coefficients. > 0. + - f_fs: frequency relative to the sampling rate, [0 ; 0.5]. + - ph_flag: set if filter is used in pi/2-phaser mode, in the form + (a - z^-2) / (1 - az^-2) +Returns: + The group delay in samples, >= 0. +Throws: Nothing +============================================================================== +*/ + +double PolyphaseIir2Designer::compute_group_delay (const double coef_arr [], int nbr_coefs, double f_fs, bool ph_flag) +{ + assert (nbr_coefs > 0); + assert (f_fs >= 0); + assert (f_fs < 0.5); + + double dly_total = 0; + for (int k = 0; k < nbr_coefs; ++k) + { + const double dly = compute_group_delay (coef_arr [k], f_fs, ph_flag); + dly_total += dly; + } + + return dly_total; +} + + + +/*\\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ + + + +/*\\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ + + + +void PolyphaseIir2Designer::compute_transition_param (double &k, double &q, double transition) +{ + assert (transition > 0); + assert (transition < 0.5); + + k = tan ((1 - transition * 2) * hiir::PI / 4); + k *= k; + assert (k < 1); + assert (k > 0); + double kksqrt = pow (1 - k * k, 0.25); + const double e = 0.5 * (1 - kksqrt) / (1 + kksqrt); + const double e2 = e * e; + const double e4 = e2 * e2; + q = e * (1 + e4 * (2 + e4 * (15 + 150 * e4))); + assert (q > 0); +} + + + +int PolyphaseIir2Designer::compute_order (double attenuation, double q) +{ + assert (attenuation > 0); + assert (q > 0); + + const double attn_p2 = pow (10.0, -attenuation / 10); + const double a = attn_p2 / (1 - attn_p2); + int order = hiir::ceil_int (log (a * a / 16) / log (q)); + if ((order & 1) == 0) + { + ++ order; + } + if (order == 1) + { + order = 3; + } + + return order; +} + + + +double PolyphaseIir2Designer::compute_atten (double q, int order) +{ + assert (q > 0); + assert (order > 0); + assert ((order & 1) == 1); + + const double a = 4 * exp (order * 0.5 * log (q)); + assert (a != -1.0); + const double attn_p2 = a / (1 + a); + const double attenuation = -10 * log10 (attn_p2); + assert (attenuation > 0); + + return attenuation; +} + + + +double PolyphaseIir2Designer::compute_coef (int index, double k, double q, int order) +{ + assert (index >= 0); + assert (index * 2 < order); + + const int c = index + 1; + const double num = compute_acc_num (q, order, c) * pow (q, 0.25); + const double den = compute_acc_den (q, order, c) + 0.5; + const double ww = num / den; + const double wwsq = ww * ww; + + const double x = sqrt ((1 - wwsq * k) * (1 - wwsq / k)) / (1 + wwsq); + const double coef = (1 - x) / (1 + x); + + return coef; +} + + + +double PolyphaseIir2Designer::compute_acc_num (double q, int order, int c) +{ + assert (c >= 1); + assert (c < order * 2); + + int i = 0; + int j = 1; + double acc = 0; + double q_ii1; + do + { + q_ii1 = hiir::ipowp (q, i * (i + 1)); + q_ii1 *= sin ((i * 2 + 1) * c * hiir::PI / order) * j; + acc += q_ii1; + + j = -j; + ++i; + } + while (fabs (q_ii1) > 1e-100); + + return acc; +} + + + +double PolyphaseIir2Designer::compute_acc_den (double q, int order, int c) +{ + assert (c >= 1); + assert (c < order * 2); + + int i = 1; + int j = -1; + double acc = 0; + double q_i2; + do + { + q_i2 = hiir::ipowp (q, i * i); + q_i2 *= cos (i * 2 * c * hiir::PI / order) * j; + acc += q_i2; + + j = -j; + ++i; + } + while (fabs (q_i2) > 1e-100); + + return acc; +} + + + +} // namespace hiir + + + +/*\\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ diff --git a/src/external/hiir/hiir/PolyphaseIir2Designer.h b/src/external/hiir/hiir/PolyphaseIir2Designer.h new file mode 100644 index 00000000..db3b4d09 --- /dev/null +++ b/src/external/hiir/hiir/PolyphaseIir2Designer.h @@ -0,0 +1,143 @@ +/***************************************************************************** + + PolyphaseIir2Designer.h + Author: Laurent de Soras, 2005 + +Compute coefficients for 2-path polyphase IIR filter, half-band filter or +Pi/2 phaser. + + -2 + a + z + N/2-1 2k +A0 (z) = Prod ---------- + k = 0 -2 + 1 + a z + 2k + + -2 + a + z + -1 (N-1)/2 2k+1 +A1 (z) = z . Prod ------------ + k = 0 -2 + 1 + a z + 2k+1 + + 1 +H (z) = - (A0 (z) + A1 (z)) + 2 + +Sum of A0 and A1 gives a low-pass filter. +Difference of A0 and A1 gives the complementary high-pass filter. + +For the Pi/2 phaser, product form is (a - z^-2) / (1 - az^-2) +Sum and difference of A0 and A1 have a Pi/2 phase difference. + +References: + +* Artur Krukowski + Polyphase Two-Path Filter Designer in Java + http://www.cmsa.wmin.ac.uk/~artur/Poly.html + +* R.A. Valenzuela, A.G. Constantinides + Digital Signal Processing Schemes for Efficient Interpolation and Decimation + IEE Proceedings, Dec 1983 + +* Scott Wardle + A Hilbert-Transformer Frequency Shifter for Audio + International Conference on Digital Audio Effects (DAFx) 1998 + http://www.iua.upf.es/dafx98/papers/WAR19.PS + +--- Legal stuff --- + +This program is free software. It comes without any warranty, to +the extent permitted by applicable law. You can redistribute it +and/or modify it under the terms of the Do What The Fuck You Want +To Public License, Version 2, as published by Sam Hocevar. See +http://sam.zoy.org/wtfpl/COPYING for more details. + +*Tab=3***********************************************************************/ + + + +#if ! defined (hiir_PolyphaseIir2Designer_HEADER_INCLUDED) +#define hiir_PolyphaseIir2Designer_HEADER_INCLUDED + +#if defined (_MSC_VER) + #pragma once + #pragma warning (4 : 4250) // "Inherits via dominance." +#endif + + + +/*\\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ + + + +namespace hiir +{ + + + +class PolyphaseIir2Designer +{ + +/*\\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ + +public: + + static int compute_nbr_coefs_from_proto (double attenuation, double transition); + static double compute_atten_from_order_tbw (int nbr_coefs, double transition); + + static int compute_coefs (double coef_arr [], double attenuation, double transition); + static void compute_coefs_spec_order_tbw (double coef_arr [], int nbr_coefs, double transition); + + static double compute_phase_delay (double a, double f_fs); + static double compute_group_delay (double a, double f_fs, bool ph_flag); + static double compute_group_delay (const double coef_arr [], int nbr_coefs, double f_fs, bool ph_flag); + + + +/*\\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ + +protected: + + + +/*\\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ + +private: + + static void compute_transition_param (double &k, double &q, double transition); + static int compute_order (double attenuation, double q); + static double compute_atten (double q, int order); + static double compute_coef (int index, double k, double q, int order); + static double compute_acc_num (double q, int order, int c); + static double compute_acc_den (double q, int order, int c); + + + +/*\\\ FORBIDDEN MEMBER FUNCTIONS \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ + +private: + + PolyphaseIir2Designer (); + ~PolyphaseIir2Designer (); + PolyphaseIir2Designer (const PolyphaseIir2Designer &other); + PolyphaseIir2Designer & + operator = (const PolyphaseIir2Designer &other); + bool operator == (const PolyphaseIir2Designer &other); + bool operator != (const PolyphaseIir2Designer &other); + +}; // class PolyphaseIir2Designer + + + +} // namespace hiir + + + +#endif // hiir_PolyphaseIir2Designer_HEADER_INCLUDED + + + +/*\\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/