X-Git-Url: https://git.gag.com/?a=blobdiff_plain;f=Misc%2Fkaiser-test.5c;fp=Misc%2Fkaiser-test.5c;h=41fd5e5c9eab56e323162ae93850858f49507895;hb=b7adf88eaf34060be679b93d0761d345543e3e97;hp=0000000000000000000000000000000000000000;hpb=77b1bf3a9c62d0ffe73be18905e26e5f8087ec4a;p=fw%2Ftmflights

diff --git a/Misc/kaiser-test.5c b/Misc/kaiser-test.5c
new file mode 100644
index 0000000..41fd5e5
--- /dev/null
+++ b/Misc/kaiser-test.5c
@@ -0,0 +1,58 @@
+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]);