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, 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 i, j, l, order = f->order;
  52   lsmp_t *x;
  53   lsmp_t *y;
  54   lsmp_t *a = f->a->data[0];
  55   lsmp_t *b = f->b->data[0];
  56
  57   for (i = 0; i < in->channels; i++) {
  58     x = f->x->data[i];
  59     y = f->y->data[i];
  60     for (j = 0; j < in->length; j++) {
  61       /* new input */
  62       x[0] = KILL_DENORMAL (in->data[i][j]);
  63       y[0] = b[0] * x[0];
  64       for (l = 1; l < order; l++) {
  65         y[0] += b[l] * x[l];
  66         y[0] -= a[l] * y[l];
  67       }
  68       /* new output */
  69       in->data[i][j] = y[0];
  70       /* store for next sample */
  71       for (l = order - 1; l > 0; l--) {
  72         x[l] = x[l - 1];
  73         y[l] = y[l - 1];
  74       }
  75     }
  76     /* store for next run */
  77     f->x->data[i] = x;
  78     f->y->data[i] = y;
  79   }
  80 }
  81
  82 /*
  83  *
  84  * despite mirroring, end effects destroy both phse and amplitude. the longer
  85  * the buffer, the less affected they are.
  86  *
  87  * replacing with zeros clicks.
  88  *
  89  * seems broken for order > 4 (see biquad_do_filtfilt for audible one)
  90  */
  91 void aubio_filter_do_filtfilt(aubio_filter_t * f, fvec_t * in, fvec_t * tmp) {
  92   uint_t j,i=0;
  93   uint_t length = in->length;
  94   //uint_t order = f->order;
  95   //lsmp_t mir;
  96   /* mirroring */
  97   //mir = 2*in->data[i][0];
  98   //for (j=1;j<order;j++)
  99   //f->x[j] = 0.;//mir - in->data[i][order-j];
 100   /* apply filtering */
 101   aubio_filter_do(f,in);
 102   /* invert */
 103   for (j = 0; j < length; j++)
 104     tmp->data[i][length-j-1] = in->data[i][j];
 105   /* mirror inverted */
 106   //mir = 2*tmp->data[i][0];
 107   //for (j=1;j<order;j++)
 108   //f->x[j] = 0.;//mir - tmp->data[i][order-j];
 109   /* apply filtering on inverted */
 110   aubio_filter_do(f,tmp);
 111   /* invert back */
 112   for (j = 0; j < length; j++)
 113     in->data[i][j] = tmp->data[i][length-j-1];
 114 }
 115
 116 lvec_t *
 117 aubio_filter_get_feedback (aubio_filter_t * f)
 118 {
 119   return f->a;
 120 }
 121
 122 lvec_t *
 123 aubio_filter_get_feedforward (aubio_filter_t * f)
 124 {
 125   return f->b;
 126 }
 127
 128 uint_t
 129 aubio_filter_get_order (aubio_filter_t * f)
 130 {
 131   return f->order;
 132 }
 133
 134 uint_t
 135 aubio_filter_get_samplerate (aubio_filter_t * f)
 136 {
 137   return f->samplerate;
 138 }
 139
 140 uint_t
 141 aubio_filter_set_samplerate (aubio_filter_t * f, uint_t samplerate)
 142 {
 143   f->samplerate = samplerate;
 144   return AUBIO_OK;
 145 }
 146
 147 aubio_filter_t *
 148 new_aubio_filter (uint_t order, uint_t channels)
 149 {
 150   aubio_filter_t *f = AUBIO_NEW (aubio_filter_t);
 151   f->x = new_lvec (order, channels);
 152   f->y = new_lvec (order, channels);
 153   f->a = new_lvec (order, 1);
 154   f->b = new_lvec (order, 1);
 155   /* by default, samplerate is not set */
 156   f->samplerate = 0;
 157   f->order = order;
 158   /* set default to identity */
 159   f->a->data[0][1] = 1.;
 160   return f;
 161 }
 162
 163 void
 164 del_aubio_filter (aubio_filter_t * f)
 165 {
 166   del_lvec (f->a);
 167   del_lvec (f->b);
 168   del_lvec (f->x);
 169   del_lvec (f->y);
 170   AUBIO_FREE (f);
 171   return;
 172 }