src/temporal/filter.c

   1 /*
   2   Copyright (C) 2003-2009 Paul Brossier <piem@aubio.org>
   3
   4   This file is part of aubio.
   5
   6   aubio is free software: you can redistribute it and/or modify
   7   it under the terms of the GNU General Public License as published by
   8   the Free Software Foundation, either version 3 of the License, or
   9   (at your option) any later version.
  10
  11   aubio is distributed in the hope that it will be useful,
  12   but WITHOUT ANY WARRANTY; without even the implied warranty of
  13   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14   GNU General Public License for more details.
  15
  16   You should have received a copy of the GNU General Public License
  17   along with aubio.  If not, see <http://www.gnu.org/licenses/>.
  18
  19 */
  20
  21
  22 /* Requires lsmp_t to be long or double. float will NOT give reliable
  23  * results */
  24
  25 #include "aubio_priv.h"
  26 #include "fvec.h"
  27 #include "lvec.h"
  28 #include "mathutils.h"
  29 #include "temporal/filter.h"
  30
  31 struct _aubio_filter_t
  32 {
  33   uint_t order;
  34   uint_t samplerate;
  35   lvec_t *a;
  36   lvec_t *b;
  37   lvec_t *y;
  38   lvec_t *x;
  39 };
  40
  41 void
  42 aubio_filter_do_outplace (aubio_filter_t * f, const fvec_t * in, fvec_t * out)
  43 {
  44   fvec_copy (in, out);
  45   aubio_filter_do (f, out);
  46 }
  47
  48 void
  49 aubio_filter_do (aubio_filter_t * f, fvec_t * in)
  50 {
  51   uint_t j, l, order = f->order;
  52   lsmp_t *x = f->x->data;
  53   lsmp_t *y = f->y->data;
  54   lsmp_t *a = f->a->data;
  55   lsmp_t *b = f->b->data;
  56
  57   for (j = 0; j < in->length; j++) {
  58     /* new input */
  59     x[0] = KILL_DENORMAL (in->data[j]);
  60     y[0] = b[0] * x[0];
  61     for (l = 1; l < order; l++) {
  62       y[0] += b[l] * x[l];
  63       y[0] -= a[l] * y[l];
  64     }
  65     /* new output */
  66     in->data[j] = y[0];
  67     /* store for next sample */
  68     for (l = order - 1; l > 0; l--) {
  69       x[l] = x[l - 1];
  70       y[l] = y[l - 1];
  71     }
  72   }
  73 }
  74
  75 /* The rough way: reset memory of filter between each run to avoid end effects. */
  76 void
  77 aubio_filter_do_filtfilt (aubio_filter_t * f, fvec_t * in, fvec_t * tmp)
  78 {
  79   uint_t j;
  80   uint_t length = in->length;
  81   /* apply filtering */
  82   aubio_filter_do (f, in);
  83   aubio_filter_do_reset (f);
  84   /* mirror */
  85   for (j = 0; j < length; j++)
  86     tmp->data[length - j - 1] = in->data[j];
  87   /* apply filtering on mirrored */
  88   aubio_filter_do (f, tmp);
  89   aubio_filter_do_reset (f);
  90   /* invert back */
  91   for (j = 0; j < length; j++)
  92     in->data[j] = tmp->data[length - j - 1];
  93 }
  94
  95 lvec_t *
  96 aubio_filter_get_feedback (const aubio_filter_t * f)
  97 {
  98   return f->a;
  99 }
 100
 101 lvec_t *
 102 aubio_filter_get_feedforward (const aubio_filter_t * f)
 103 {
 104   return f->b;
 105 }
 106
 107 uint_t
 108 aubio_filter_get_order (const aubio_filter_t * f)
 109 {
 110   return f->order;
 111 }
 112
 113 uint_t
 114 aubio_filter_get_samplerate (const aubio_filter_t * f)
 115 {
 116   return f->samplerate;
 117 }
 118
 119 uint_t
 120 aubio_filter_set_samplerate (aubio_filter_t * f, uint_t samplerate)
 121 {
 122   f->samplerate = samplerate;
 123   return AUBIO_OK;
 124 }
 125
 126 void
 127 aubio_filter_do_reset (aubio_filter_t * f)
 128 {
 129   lvec_zeros (f->x);
 130   lvec_zeros (f->y);
 131 }
 132
 133 aubio_filter_t *
 134 new_aubio_filter (uint_t order)
 135 {
 136   aubio_filter_t *f = AUBIO_NEW (aubio_filter_t);
 137   if ((sint_t)order < 1) {
 138     AUBIO_FREE(f);
 139     return NULL;
 140   }
 141   f->x = new_lvec (order);
 142   f->y = new_lvec (order);
 143   f->a = new_lvec (order);
 144   f->b = new_lvec (order);
 145   /* by default, samplerate is not set */
 146   f->samplerate = 0;
 147   f->order = order;
 148   /* set default to identity */
 149   f->a->data[0] = 1.;
 150   f->b->data[0] = 1.;
 151   return f;
 152 }
 153
 154 void
 155 del_aubio_filter (aubio_filter_t * f)
 156 {
 157   del_lvec (f->a);
 158   del_lvec (f->b);
 159   del_lvec (f->x);
 160   del_lvec (f->y);
 161   AUBIO_FREE (f);
 162   return;
 163 }