<commit_before>// Licensed GNU LGPL v2.1 or later: http://www.gnu.org/licenses/lgpl.html
#include <bse/bsedefs.hh>
// #define TEST_VERBOSE
#include <sfi/sfitests.hh>
#include <bse/bsefilter.hh>
#include <bse/bsemain.hh>
#include <bse/gslfilter.hh> // FIXME
#include <bse/bseglobals.hh> // FIXME
#include "topconfig.h"
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <complex>
using namespace Rapicorn;
using std::max;
using std::min;

typedef std::complex<double> Complex;

static inline double
sqr (register double a)
{
  return a * a;
}

static inline uint
complex_find_nearest (const BseComplex *zp,
                      uint              n_zps,
                      const Complex    *zps)
{
  const Complex z = Complex (zp->re, zp->im);
  uint j = 0;
  double last = sqr (real (z) - real (zps[j])) + sqr (imag (z) - imag (zps[j]));
  if (last == 0.0)
    return j;
  for (uint i = 1; i < n_zps; i++)
    {
      double dist2 = sqr (real (z) - real (zps[i])) + sqr (imag (z) - imag (zps[i]));
      if (dist2 < last)
        {
          last = dist2;
          j = i;
          if (last == 0.0)
            return j;
        }
    }
  return j;
}

static double
compare_zeros (uint              n_zeros,
               const BseComplex *czeros,
               const double     *rizeros)
{
  Complex *z2 = (Complex*) alloca (sizeof (Complex) * n_zeros);
  for (uint i = 0; i < n_zeros; i++)
    z2[i] = Complex (rizeros[i * 2], rizeros[i * 2 + 1]);
  double max_eps = 0;
  for (uint i = 0; i < n_zeros; i++)
    {
      uint j = complex_find_nearest (&czeros[i], n_zeros - i, z2);
      double reps = fabs (real (z2[j]) - czeros[i].re);
      double ieps = fabs (imag (z2[j]) - czeros[i].im);
      z2[j] = z2[n_zeros - i - 1]; // optimization of: swap (z2[last], z2[j]);
      max_eps = max (max_eps, max (reps, ieps));
    }
  return max_eps;
}

static double
to_db (double response)
{
  return response <= 0.0 ? -999.99 : max (BSE_DECIBEL20_FACTOR * log (response), -999.99);
}

static double
filter_zp_response (const BseIIRFilterDesign *fdes,
                    double                    freq)
{
  /* map freq onto unit circle, eval filter zeros and poles in z-plane */
  double omega = 2.0 * PI * freq / fdes->sampling_frequency;    // map freq onto unit circle circumference
  Complex z (cos (omega), sin (omega));                         // map freq onto unit circle (into z-plane): exp (j * omega)
  Complex num (1, 0), den (1, 0);
  for (uint i = 0; i < fdes->n_zeros; i++)
    {
      num = num * (z - Complex (fdes->zz[i].re, fdes->zz[i].im));
      if (fdes->zz[i].im != 0.0)
        num = num * (z - std::conj (Complex (fdes->zz[i].re, fdes->zz[i].im)));
    }
  for (uint i = 0; i < fdes->n_poles; i++)
    {
      den = den * (z - Complex (fdes->zp[i].re, fdes->zp[i].im));
      if (fdes->zp[i].im != 0.0)
        den = den * (z - std::conj (Complex (fdes->zp[i].re, fdes->zp[i].im)));
    }
  Complex r = fdes->gain * num / den;
  double freq_magnitude = abs (r);
  double freq_phase = arg (r);
  (void) freq_phase;
  return freq_magnitude;
}

bool
bse_iir_filter_dump_gnuplot (const BseIIRFilterDesign *fdes,
                             const char               *fname,
                             uint                      n_consts,
                             const double             *consts,<commit_msg>BSE: tests: minor compilation fix<commit_after>// Licensed GNU LGPL v2.1 or later: http://www.gnu.org/licenses/lgpl.html
#include <bse/bsedefs.hh>
// #define TEST_VERBOSE
#include <sfi/sfitests.hh>
#include <bse/bsefilter.hh>
#include <bse/bsemain.hh>
#include <bse/gslfilter.hh> // FIXME
#include <bse/bseglobals.hh> // FIXME
#include "topconfig.h"
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <complex>
using namespace Rapicorn;
using std::max;
using std::min;

typedef std::complex<double> Complex;

static inline double
sqr (register double a)
{
  return a * a;
}

static inline uint
complex_find_nearest (const BseComplex *zp,
                      uint              n_zps,
                      const Complex    *zps)
{
  const Complex z = Complex (zp->re, zp->im);
  uint j = 0;
  double last = sqr (real (z) - real (zps[j])) + sqr (imag (z) - imag (zps[j]));
  if (last == 0.0)
    return j;
  for (uint i = 1; i < n_zps; i++)
    {
      double dist2 = sqr (real (z) - real (zps[i])) + sqr (imag (z) - imag (zps[i]));
      if (dist2 < last)
        {
          last = dist2;
          j = i;
          if (last == 0.0)
            return j;
        }
    }
  return j;
}

static double
compare_zeros (uint              n_zeros,
               const BseComplex *czeros,
               const double     *rizeros)
{
  Complex *z2 = (Complex*) alloca (sizeof (Complex) * n_zeros);
  for (uint i = 0; i < n_zeros; i++)
    z2[i] = Complex (rizeros[i * 2], rizeros[i * 2 + 1]);
  double max_eps = 0;
  for (uint i = 0; i < n_zeros; i++)
    {
      uint j = complex_find_nearest (&czeros[i], n_zeros - i, z2);
      double reps = fabs (real (z2[j]) - czeros[i].re);
      double ieps = fabs (imag (z2[j]) - czeros[i].im);
      z2[j] = z2[n_zeros - i - 1]; // optimization of: swap (z2[last], z2[j]);
      max_eps = max (max_eps, max (reps, ieps));
    }
  return max_eps;
}

static double
to_db (double response)
{
  return response <= 0.0 ? -999.99 : max (BSE_DECIBEL20_FACTOR * log (response), -999.99);
}

static double
filter_zp_response (const BseIIRFilterDesign *fdes,
                    double                    freq)
{
  /* map freq onto unit circle, eval filter zeros and poles in z-plane */
  double omega = 2.0 * PI * freq / fdes->sampling_frequency;    // map freq onto unit circle circumference
  Complex z (cos (omega), sin (omega));                         // map freq onto unit circle (into z-plane): exp (j * omega)
  Complex num (1, 0), den (1, 0);
  for (uint i = 0; i < fdes->n_zeros; i++)
    {
      num = num * (z - Complex (fdes->zz[i].re, fdes->zz[i].im));
      if (fdes->zz[i].im != 0.0)
        num = num * (z - std::conj (Complex (fdes->zz[i].re, fdes->zz[i].im)));
    }
  for (uint i = 0; i < fdes->n_poles; i++)
    {
      den = den * (z - Complex (fdes->zp[i].re, fdes->zp[i].im));
      if (fdes->zp[i].im != 0.0)
        den = den * (z - std::conj (Complex (fdes->zp[i].re, fdes->zp[i].im)));
    }
  Complex r = fdes->gain * num / den;
  double freq_magnitude = abs (r);
  double freq_phase = arg (r);
  (void) freq_phase;
  return freq_magnitude;
}

static bool
bse_iir_filter_dump_gnuplot (const BseIIRFilterDesign *fdes,
                             const char               *fname,
                             uint                      n_consts,
                             const double             *consts,
{'fraction_non_alphanumeric': 0.15421185559352252, 'fraction_numerical': 0.015748031496062992, 'mean_word_length': 2.0994475138121547, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 57, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3913306', 'n_tokens_mistral': 2603, 'n_tokens_neox': 2508, 'n_words': 700}