src/firfit.c

   1 /* Effect: firfit filter     Copyright (c) 2009 robs@users.sourceforge.net
   2  *
   3  * This library is free software; you can redistribute it and/or modify it
   4  * under the terms of the GNU Lesser General Public License as published by
   5  * the Free Software Foundation; either version 2.1 of the License, or (at
   6  * your option) any later version.
   7  *
   8  * This library is distributed in the hope that it will be useful, but
   9  * WITHOUT ANY WARRANTY; without even the implied warranty of
  10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
  11  * General Public License for more details.
  12  *
  13  * You should have received a copy of the GNU Lesser General Public License
  14  * along with this library; if not, write to the Free Software Foundation,
  15  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  16  */
  17
  18 /* This is W.I.P. hence marked SOX_EFF_ALPHA for now.
  19  * Need to add other interpolation types e.g. linear, bspline, window types,
  20  * and filter length, maybe phase response type too.
  21  */
  22
  23 #include "sox_i.h"
  24 #include "dft_filter.h"
  25
  26 typedef struct {
  27   dft_filter_priv_t base;
  28   char const         * filename;
  29   struct {double f, gain;} * knots;
  30   int        num_knots, n;
  31 } priv_t;
  32
  33 static int create(sox_effect_t * effp, int argc, char **argv)
  34 {
  35   priv_t * p = (priv_t *)effp->priv;
  36   dft_filter_priv_t * b = &p->base;
  37   b->filter_ptr = &b->filter;
  38   --argc, ++argv;
  39   if (argc == 1)
  40     p->filename = argv[0], --argc;
  41   p->n = 2047;
  42   return argc? lsx_usage(effp) : SOX_SUCCESS;
  43 }
  44
  45 static double * make_filter(sox_effect_t * effp)
  46 {
  47   priv_t * p = (priv_t *)effp->priv;
  48   double * log_freqs, * gains, * d, * work, * h;
  49   sox_rate_t rate = effp->in_signal.rate;
  50   int i, work_len;
  51
  52   lsx_valloc(log_freqs , p->num_knots);
  53   lsx_valloc(gains, p->num_knots);
  54   lsx_valloc(d  , p->num_knots);
  55   for (i = 0; i < p->num_knots; ++i) {
  56     log_freqs[i] = log(max(p->knots[i].f, 1));
  57     gains[i] = p->knots[i].gain;
  58   }
  59   lsx_prepare_spline3(log_freqs, gains, p->num_knots, HUGE_VAL, HUGE_VAL, d);
  60
  61   for (work_len = 8192; work_len < rate / 2; work_len <<= 1);
  62   work = lsx_calloc(work_len + 2, sizeof(*work));
  63   lsx_valloc(h, p->n);
  64
  65   for (i = 0; i <= work_len; i += 2) {
  66     double f = rate * 0.5 * i / work_len;
  67     double spl1 = f < max(p->knots[0].f, 1)? gains[0] :
  68                   f > p->knots[p->num_knots - 1].f? gains[p->num_knots - 1] :
  69                   lsx_spline3(log_freqs, gains, d, p->num_knots, log(f));
  70     work[i] = dB_to_linear(spl1);
  71   }
  72   LSX_PACK(work, work_len);
  73   lsx_safe_rdft(work_len, -1, work);
  74   for (i = 0; i < p->n; ++i)
  75     h[i] = work[(work_len - p->n / 2 + i) % work_len] * 2. / work_len;
  76   lsx_apply_blackman_nutall(h, p->n);
  77
  78   free(work);
  79   return h;
  80 }
  81
  82 static sox_bool read_knots(sox_effect_t * effp)
  83 {
  84   priv_t * p = (priv_t *) effp->priv;
  85   FILE * file = lsx_open_input_file(effp, p->filename, sox_true);
  86   sox_bool result = sox_false;
  87   int num_converted = 1;
  88   char c;
  89
  90   if (file) {
  91     lsx_valloc(p->knots, 1);
  92     while (fscanf(file, " #%*[^\n]%c", &c) >= 0) {
  93       num_converted = fscanf(file, "%lf %lf",
  94           &p->knots[p->num_knots].f, &p->knots[p->num_knots].gain);
  95       if (num_converted == 2) {
  96         if (p->num_knots && p->knots[p->num_knots].f <= p->knots[p->num_knots - 1].f) {
  97           lsx_fail("knot frequencies must be strictly increasing");
  98           break;
  99         }
 100         lsx_revalloc(p->knots, ++p->num_knots + 1);
 101       } else if (num_converted != 0)
 102         break;
 103     }
 104     lsx_report("%i knots", p->num_knots);
 105     if (feof(file) && num_converted != 1)
 106       result = sox_true;
 107     else lsx_fail("error reading knot file `%s', line number %u", p->filename, 1 + p->num_knots);
 108     if (file != stdin)
 109       fclose(file);
 110   }
 111   return result;
 112 }
 113
 114 static int start(sox_effect_t * effp)
 115 {
 116   priv_t * p = (priv_t *) effp->priv;
 117   dft_filter_t * f = p->base.filter_ptr;
 118
 119   if (!f->num_taps) {
 120     double * h;
 121     if (!p->num_knots && !read_knots(effp))
 122       return SOX_EOF;
 123     h = make_filter(effp);
 124     if (effp->global_info->plot != sox_plot_off) {
 125       lsx_plot_fir(h, p->n, effp->in_signal.rate,
 126           effp->global_info->plot, "SoX effect: firfit", -30., +30.);
 127       return SOX_EOF;
 128     }
 129     lsx_set_dft_filter(f, h, p->n, p->n >> 1);
 130   }
 131   return lsx_dft_filter_effect_fn()->start(effp);
 132 }
 133
 134 sox_effect_handler_t const * lsx_firfit_effect_fn(void)
 135 {
 136   static sox_effect_handler_t handler;
 137   handler = *lsx_dft_filter_effect_fn();
 138   handler.name = "firfit";
 139   handler.usage = "[knots-file]";
 140   handler.flags |= SOX_EFF_ALPHA;
 141   handler.getopts = create;
 142   handler.start = start;
 143   handler.priv_size = sizeof(priv_t);
 144   return &handler;
 145 }