git.gag.com Git - fw/tmflights/blob - kaiser-test.5c

   1 load "kaiser.5c"
   2
   3 real sinc(real x) = x != 0 ? sin(x)/x : 1;
   4
   5 real sum(real[...] x) { real s = 0; for(int i = 0; i < dim(x); i++) s += x[i]; return s; }
   6
   7 real[...] convolve(real[...] d, real[...] e) {
   8         real sample(n) = n < 0 ? d[0] : n >= dim(d) ? d[dim(d)-1] : d[n];
   9         real w = floor ((dim(e) - 1) / 2);
  10         real c(int center) {
  11                 real    v = 0;
  12                 real    s = 0;
  13                 int     start;
  14                 int     end;
  15
  16                 start = center - w;
  17                 if (start < 0)
  18                         start = 0;
  19                 end = center + w;
  20                 if (end >= dim(d))
  21                         end = dim(d) - 1;
  22                 for (int o = start; o <= end; o++) {
  23                         real    e_o = e[o - center + w];
  24                         real    part = d[o] * e_o;
  25                         v += part;
  26                         s += e_o;
  27                 }
  28                 v /= s;
  29                 return v;
  30         }
  31         return (real[dim(d)]) { [n] = c(n) };
  32 }
  33
  34 real[...] kaiser_filter(real[...] d, int half_width) {
  35 #       real[half_width * 2 + 1] fir = { [n] = sinc(2 * pi * n / (2 * half_width)) };
  36         real M = half_width * 2 + 1;
  37         real[M] fir = { [n] = kaiser(n+1, M+1, 8) };
  38         return convolve(d, fir);
  39 }
  40
  41
  42 real[1024] square = { [n] = 2 * ((n >> 8) & 1) - 1 };
  43
  44 real[] cut = low_pass_filter(π*1/128, π*1/64, 1e-10);
  45
  46 real[dim(square)] filtered_square = convolve(square, cut);
  47
  48 for (int i = 0; i < dim(filtered_square); i++)
  49     printf ("%d %g %g\n", i, square[i], filtered_square[i]);
  50
  51 exit(0);
  52
  53 int size = 100;
  54 real[size] raw = { [n] = n };
  55 real[size] filtered = kaiser_filter(raw, 8);
  56
  57 for (int i = 0; i < size; i++)
  58         printf ("%d %g %g\n", i, raw[i], filtered[i]);