load "kaiser.5c"

real sinc(real x) = x != 0 ? sin(x)/x : 1;

real sum(real[...] x) { real s = 0; for(int i = 0; i < dim(x); i++) s += x[i]; return s; }

real[...] convolve(real[...] d, real[...] e) {
	real sample(n) = n < 0 ? d[0] : n >= dim(d) ? d[dim(d)-1] : d[n];
	real w = floor ((dim(e) - 1) / 2);
	real c(int center) {
		real	v = 0;
		real	s = 0;
		int	start;
		int	end;

		start = center - w;
		if (start < 0)
			start = 0;
		end = center + w;
		if (end >= dim(d))
			end = dim(d) - 1;
		for (int o = start; o <= end; o++) {
			real	e_o = e[o - center + w];
			real	part = d[o] * e_o;
			v += part;
			s += e_o;
		}
		v /= s;
		return v;
	}
	return (real[dim(d)]) { [n] = c(n) };
}

real[...] kaiser_filter(real[...] d, int half_width) {
#	real[half_width * 2 + 1] fir = { [n] = sinc(2 * pi * n / (2 * half_width)) };
	real M = half_width * 2 + 1;
	real[M] fir = { [n] = kaiser(n+1, M+1, 8) };
	return convolve(d, fir);
}


real[1024] square = { [n] = 2 * ((n >> 8) & 1) - 1 };

real[] cut = low_pass_filter(π*1/128, π*1/64, 1e-10);

real[dim(square)] filtered_square = convolve(square, cut);

for (int i = 0; i < dim(filtered_square); i++)
    printf ("%d %g %g\n", i, square[i], filtered_square[i]);

exit(0);

int size = 100;
real[size] raw = { [n] = n };
real[size] filtered = kaiser_filter(raw, 8);

for (int i = 0; i < size; i++)
	printf ("%d %g %g\n", i, raw[i], filtered[i]);