src/sinc.c

   1 /* Effect: sinc filters     Copyright (c) 2008-9 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 #include "sox_i.h"
  19 #include "dft_filter.h"
  20 #include <string.h>
  21
  22 typedef struct {
  23   dft_filter_priv_t  base;
  24   double             att, beta, phase, Fc0, Fc1, tbw0, tbw1;
  25   int                num_taps[2];
  26   sox_bool           round;
  27 } priv_t;
  28
  29 static int create(sox_effect_t * effp, int argc, char * * argv)
  30 {
  31   priv_t * p = (priv_t *)effp->priv;
  32   dft_filter_priv_t * b = &p->base;
  33   char * parse_ptr = argv[0];
  34   int i = 0;
  35   lsx_getopt_t optstate;
  36   lsx_getopt_init(argc, argv, "+ra:b:p:MILt:n:", NULL, lsx_getopt_flag_none, 1, &optstate);
  37
  38   b->filter_ptr = &b->filter;
  39   p->phase = 50;
  40   p->beta = -1;
  41   while (i < 2) {
  42     int c = 1;
  43     while (c && (c = lsx_getopt(&optstate)) != -1) switch (c) {
  44       char * parse_ptr2;
  45       case 'r': p->round = sox_true; break;
  46       GETOPT_NUMERIC(optstate, 'a', att,  40 , 180)
  47       GETOPT_NUMERIC(optstate, 'b', beta,  0 , 256)
  48       GETOPT_NUMERIC(optstate, 'p', phase, 0, 100)
  49       case 'M': p->phase =  0; break;
  50       case 'I': p->phase = 25; break;
  51       case 'L': p->phase = 50; break;
  52       GETOPT_NUMERIC(optstate, 'n', num_taps[1], 11, 32767)
  53       case 't': p->tbw1 = lsx_parse_frequency(optstate.arg, &parse_ptr2);
  54         if (p->tbw1 < 1 || *parse_ptr2) return lsx_usage(effp);
  55         break;
  56       default: c = 0;
  57     }
  58     if ((p->att && p->beta >= 0) || (p->tbw1 && p->num_taps[1]))
  59       return lsx_usage(effp);
  60     if (!i || !p->Fc1)
  61       p->tbw0 = p->tbw1, p->num_taps[0] = p->num_taps[1];
  62     if (!i++ && optstate.ind < argc) {
  63       if (*(parse_ptr = argv[optstate.ind++]) != '-')
  64         p->Fc0 = lsx_parse_frequency(parse_ptr, &parse_ptr);
  65       if (*parse_ptr == '-')
  66         p->Fc1 = lsx_parse_frequency(parse_ptr + 1, &parse_ptr);
  67     }
  68   }
  69   return optstate.ind != argc || p->Fc0 < 0 || p->Fc1 < 0 || *parse_ptr ?
  70       lsx_usage(effp) : SOX_SUCCESS;
  71 }
  72
  73 static void invert(double * h, int n)
  74 {
  75   int i;
  76   for (i = 0; i < n; ++i)
  77     h[i] = -h[i];
  78   h[(n - 1) / 2] += 1;
  79 }
  80
  81 static double * lpf(double Fn, double Fc, double tbw, int * num_taps, double att, double * beta, sox_bool round)
  82 {
  83   int n = *num_taps;
  84   if ((Fc /= Fn) <= 0 || Fc >= 1) {
  85     *num_taps = 0;
  86     return NULL;
  87   }
  88   att = att? att : 120;
  89   lsx_kaiser_params(att, Fc, (tbw? tbw / Fn : .05) * .5, beta, num_taps);
  90   if (!n) {
  91     n = *num_taps;
  92     *num_taps = range_limit(n, 11, 32767);
  93     if (round)
  94       *num_taps = 1 + 2 * (int)((int)((*num_taps / 2) * Fc + .5) / Fc + .5);
  95     lsx_report("num taps = %i (from %i)", *num_taps, n);
  96   }
  97   return lsx_make_lpf(*num_taps |= 1, Fc, *beta, 0., 1., sox_false);
  98 }
  99
 100 static int start(sox_effect_t * effp)
 101 {
 102   priv_t * p = (priv_t *)effp->priv;
 103   dft_filter_t * f = p->base.filter_ptr;
 104
 105   if (!f->num_taps) {
 106     double Fn = effp->in_signal.rate * .5;
 107     double * h[2];
 108     int i, n, post_peak, longer;
 109
 110     if (p->Fc0 >= Fn || p->Fc1 >= Fn) {
 111       lsx_fail("filter frequency must be less than sample-rate / 2");
 112       return SOX_EOF;
 113     }
 114     h[0] = lpf(Fn, p->Fc0, p->tbw0, &p->num_taps[0], p->att, &p->beta,p->round);
 115     h[1] = lpf(Fn, p->Fc1, p->tbw1, &p->num_taps[1], p->att, &p->beta,p->round);
 116     if (h[0])
 117       invert(h[0], p->num_taps[0]);
 118
 119     longer = p->num_taps[1] > p->num_taps[0];
 120     n = p->num_taps[longer];
 121     if (h[0] && h[1]) {
 122       for (i = 0; i < p->num_taps[!longer]; ++i)
 123         h[longer][i + (n - p->num_taps[!longer])/2] += h[!longer][i];
 124
 125       if (p->Fc0 < p->Fc1)
 126         invert(h[longer], n);
 127
 128       free(h[!longer]);
 129     }
 130     if (p->phase != 50)
 131       lsx_fir_to_phase(&h[longer], &n, &post_peak, p->phase);
 132     else post_peak = n >> 1;
 133
 134     if (effp->global_info->plot != sox_plot_off) {
 135       char title[100];
 136       sprintf(title, "SoX effect: sinc filter freq=%g-%g",
 137           p->Fc0, p->Fc1? p->Fc1 : Fn);
 138       lsx_plot_fir(h[longer], n, effp->in_signal.rate,
 139           effp->global_info->plot, title, -p->beta * 10 - 25, 5.);
 140       return SOX_EOF;
 141     }
 142     lsx_set_dft_filter(f, h[longer], n, post_peak);
 143   }
 144   return lsx_dft_filter_effect_fn()->start(effp);
 145 }
 146
 147 sox_effect_handler_t const * lsx_sinc_effect_fn(void)
 148 {
 149   static sox_effect_handler_t handler;
 150   handler = *lsx_dft_filter_effect_fn();
 151   handler.name = "sinc";
 152   handler.usage = "[-a att|-b beta] [-p phase|-M|-I|-L] [-t tbw|-n taps] [freqHP][-freqLP [-t tbw|-n taps]]";
 153   handler.getopts = create;
 154   handler.start = start;
 155   handler.priv_size = sizeof(priv_t);
 156   return &handler;
 157 }