3 real sinc(real x) = x != 0 ? sin(x)/x : 1;
5 real sum(real[...] x) { real s = 0; for(int i = 0; i < dim(x); i++) s += x[i]; return s; }
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);
22 for (int o = start; o <= end; o++) {
23 real e_o = e[o - center + w];
24 real part = d[o] * e_o;
31 return (real[dim(d)]) { [n] = c(n) };
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);
42 real[1024] square = { [n] = 2 * ((n >> 8) & 1) - 1 };
44 real[] cut = low_pass_filter(π*1/128, π*1/64, 1e-10);
46 real[dim(square)] filtered_square = convolve(square, cut);
48 for (int i = 0; i < dim(filtered_square); i++)
49 printf ("%d %g %g\n", i, square[i], filtered_square[i]);
54 real[size] raw = { [n] = n };
55 real[size] filtered = kaiser_filter(raw, 8);
57 for (int i = 0; i < size; i++)
58 printf ("%d %g %g\n", i, raw[i], filtered[i]);