src/temporal/filter.c

   1 /*
   2    Copyright (C) 2003 Paul Brossier
   3
   4    This program is free software; you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation; either version 2 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program; if not, write to the Free Software
  16    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  17
  18 */
  19
  20
  21 /* Requires lsmp_t to be long or double. float will NOT give reliable
  22  * results */
  23
  24 #include "aubio_priv.h"
  25 #include "fvec.h"
  26 #include "lvec.h"
  27 #include "mathutils.h"
  28 #include "temporal/filter.h"
  29 #include "temporal/filter_priv.h"
  30
  31 void aubio_filter_do(aubio_filter_t * f, fvec_t * in) {
  32   aubio_filter_do_outplace(f, in, in);
  33 }
  34
  35 void aubio_filter_do_outplace(aubio_filter_t * f, fvec_t * in, fvec_t * out) {
  36   uint_t i,j,l, order = f->order;
  37   lsmp_t *x;
  38   lsmp_t *y;
  39   lsmp_t *a = f->a->data[0];
  40   lsmp_t *b = f->b->data[0];
  41
  42   for (i = 0; i < in->channels; i++) {
  43     x = f->x->data[i];
  44     y = f->y->data[i];
  45     for (j = 0; j < in->length; j++) {
  46       /* new input */
  47       if (ISDENORMAL(in->data[i][j])) {
  48         x[0] = y[0] = 0.;
  49       } else {
  50         x[0] = in->data[i][j];
  51         y[0] = b[0] * x[0];
  52         for (l=1;l<order; l++) {
  53           y[0] += b[l] * x[l];
  54           y[0] -= a[l] * y[l];
  55         }
  56       }
  57       /* new output */
  58       out->data[i][j] = y[0];
  59       /* store for next sample */
  60       for (l=order-1; l>0; l--){
  61         x[l] = x[l-1];
  62         y[l] = y[l-1];
  63       }
  64     }
  65     /* store for next run */
  66     f->x->data[i] = x;
  67     f->y->data[i] = y;
  68   }
  69 }
  70
  71 /*
  72  *
  73  * despite mirroring, end effects destroy both phse and amplitude. the longer
  74  * the buffer, the less affected they are.
  75  *
  76  * replacing with zeros clicks.
  77  *
  78  * seems broken for order > 4 (see biquad_do_filtfilt for audible one)
  79  */
  80 void aubio_filter_do_filtfilt(aubio_filter_t * f, fvec_t * in, fvec_t * tmp) {
  81   uint_t j,i=0;
  82   uint_t length = in->length;
  83   //uint_t order = f->order;
  84   //lsmp_t mir;
  85   /* mirroring */
  86   //mir = 2*in->data[i][0];
  87   //for (j=1;j<order;j++)
  88   //f->x[j] = 0.;//mir - in->data[i][order-j];
  89   /* apply filtering */
  90   aubio_filter_do(f,in);
  91   /* invert */
  92   for (j = 0; j < length; j++)
  93     tmp->data[i][length-j-1] = in->data[i][j];
  94   /* mirror inverted */
  95   //mir = 2*tmp->data[i][0];
  96   //for (j=1;j<order;j++)
  97   //f->x[j] = 0.;//mir - tmp->data[i][order-j];
  98   /* apply filtering on inverted */
  99   aubio_filter_do(f,tmp);
 100   /* invert back */
 101   for (j = 0; j < length; j++)
 102     in->data[i][j] = tmp->data[i][length-j-1];
 103 }
 104
 105 aubio_filter_t * new_aubio_filter(uint_t samplerate UNUSED, uint_t order, uint_t channels) {
 106   aubio_filter_t * f = AUBIO_NEW(aubio_filter_t);
 107   f->x = new_lvec(order, channels);
 108   f->y = new_lvec(order, channels);
 109   f->a = new_lvec(order, 1);
 110   f->b = new_lvec(order, 1);
 111   f->a->data[0][1] = 1.;
 112   f->order = order;
 113   return f;
 114 }
 115
 116 void del_aubio_filter(aubio_filter_t * f) {
 117   AUBIO_FREE(f->a);
 118   AUBIO_FREE(f->b);
 119   AUBIO_FREE(f->x);
 120   AUBIO_FREE(f->y);
 121   AUBIO_FREE(f);
 122   return;
 123 }