src/out_sfilt.c

   1 /* copyright 2016 Apache 2 sddekit authors */
   2
   3 #include "sddekit.h"
   4
   5 /* spatial filter (bank)
   6  *
   7  * xfilts laid out as filter per row, i.e. shape [nfilt nx] & [nfilt nc]
   8  */
   9
  10 typedef struct sfiltd {
  11         sd_out_sfilt *sfilt;
  12         sd_out *out;
  13         uint32_t nfilt, filtlen;
  14         double * restrict xfilts, * restrict cfilts, * restrict x, * restrict c;
  15         sd_out *next_out;
  16 } sfiltd;
  17
  18 uint32_t get_nfilt(sd_out_sfilt *out) { return ((sfiltd*)out->ptr)->nfilt; }
  19 uint32_t get_filtlen(sd_out_sfilt *out) { return ((sfiltd*)out->ptr)->filtlen; }
  20 double *get_xfilts(sd_out_sfilt *out) { return ((sfiltd*)out->ptr)->xfilts; }
  21 double *get_cfilts(sd_out_sfilt *out) { return ((sfiltd*)out->ptr)->cfilts; }
  22 double *get_x(sd_out_sfilt *out) { return ((sfiltd*)out->ptr)->x; }
  23 double *get_c(sd_out_sfilt *out) { return ((sfiltd*)out->ptr)->c; }
  24 sd_out *get_out(sd_out_sfilt *out) { return ((sfiltd*)out->ptr)->out; }
  25 sd_out *get_next_out(sd_out_sfilt *out) { return ((sfiltd*)out->ptr)->next_out; }
  26
  27 static sd_stat apply(sd_out *out, double t,
  28                    uint32_t nx, double * restrict x,
  29                    uint32_t nc, double * restrict c)
  30 {
  31         (void) nx; (void) nc;
  32         uint32_t i, j;
  33         sfiltd *d = out->ptr;
  34         for (i=0; i<d->nfilt; i++) {
  35                 d->x[i] = d->c[i] = 0.0;
  36                 for (j=0; j<d->filtlen; j++) {
  37                         d->x[i] += d->xfilts[i*d->filtlen+j] * x[j];
  38                         d->c[i] += d->cfilts[i*d->filtlen+j] * c[j];
  39                 }
  40         }
  41         return d->next_out->apply(d->next_out, t, d->nfilt, d->x, d->nfilt, d->c);
  42 }
  43
  44 void sfilt_free(sd_out *out) {
  45         sfiltd *d = out->ptr;
  46         sd_free(d->xfilts);
  47         sd_free(d->cfilts);
  48         sd_free(d->x);
  49         sd_free(d->c);
  50         sd_free(d->out);
  51         sd_free(d->sfilt);
  52         sd_free(d);
  53 }
  54
  55 static sd_out sfilt_out_defaults = {
  56         .ptr = NULL,
  57         .apply = &apply,
  58         .free = &sfilt_free
  59 };
  60
  61 static sd_out_sfilt sfilt_defaults = {
  62         .ptr = NULL,
  63         .get_nfilt = &get_nfilt,
  64         .get_filtlen = &get_filtlen,
  65         .get_xfilts = &get_xfilts,
  66         .get_cfilts = &get_cfilts,
  67         .get_x = &get_x,
  68         .get_c = &get_c,
  69         .out = &get_out,
  70         .get_out = &get_next_out,
  71 };
  72
  73 sd_out_sfilt *
  74 sd_out_sfilt_new(uint32_t nfilt, uint32_t filtlen,
  75                                  double * restrict xfilts, double * restrict cfilts,
  76                              sd_out *next_out)
  77 {
  78         uint32_t n_filt_bytes;
  79         sfiltd *d, zd = { 0 };
  80         n_filt_bytes = nfilt * filtlen * sizeof(double);
  81         if ((d = sd_malloc(sizeof(sfiltd))) == NULL
  82                 || (*d = zd, nfilt < 1 || filtlen < 1)
  83                 || (d->out = sd_malloc(sizeof(sd_out))) == NULL
  84                 || (d->sfilt = sd_malloc(sizeof(sd_out_sfilt))) == NULL
  85                 || (d->xfilts = sd_malloc(n_filt_bytes)) == NULL
  86                 || (d->cfilts = sd_malloc(n_filt_bytes)) == NULL
  87                 || (d->x = sd_malloc(sizeof(double) * nfilt)) == NULL
  88                 || (d->c = sd_malloc(sizeof(double) * nfilt)) == NULL
  89                 )
  90         {
  91                 if (d->xfilts!=NULL) sd_free(d->xfilts);
  92                 if (d->cfilts!=NULL) sd_free(d->cfilts);
  93                 if (d->x!=NULL) sd_free(d->x);
  94                 if (d->sfilt != NULL) sd_free(d->sfilt);
  95                 if (d->out != NULL) sd_free(d->out);
  96                 sd_free(d);
  97                 sd_err("nfilt < 1, filtlen < 1 or memory alloc failed.");
  98                 return NULL;
  99         }
 100         d->nfilt = nfilt;
 101         d->filtlen = filtlen;
 102         memcpy(d->xfilts, xfilts, n_filt_bytes);
 103         memcpy(d->cfilts, cfilts, n_filt_bytes);
 104         d->next_out = next_out;
 105         *(d->out) = sfilt_out_defaults;
 106         *(d->sfilt) = sfilt_defaults;
 107         d->out->ptr = d->sfilt->ptr = d;
 108         return d->sfilt;
 109 }